Как реализовать Алгоритм Решето Эратосфена от определенного числа, до данного?

Question

[Idwln]

У меня есть код Алгоритма Решето Эратосфена:

def primes(n):
        sieve = [True] * n
        for i in range(3, int(n ** 0.5) + 1, 2):
            if sieve[i]:
                sieve[i * i::2 * i] = [False] * ((n - i * i - 1) // (2 * i) + 1)
        return [2] + [i for i in range(3, n, 2) if sieve[i]]

Как мне его исправить так, чтобы если я вызываю эту функцию и передаю туда число 1000000000 поиск простых чисел начинался с числа 900000000. Сам не могу понять, ибо код Решето брал с интернета и для моих задач он не требовал изменений

Comments

[mayton2019]

Скорее всего это невозможно. Решето - это континуум. И чтобы продолжить с 900 миллинов до мильярда - тебе нужно предыдущее состояние этой системы а именно сведения о всех невычеркнутых числах от 2 до 9000.

Поищи другие алгоритмы не на основе биткарты.

[Wataru]

mayton2019, Есть алгоритм сегментированное решето.

Достаточно найти все простые до корня из 1000000000, а дальше их надо вручную вычеркнуть из текущего сегмента.

[mayton2019]

Wataru, в этом случае - зачем вам решето? У вас есть коллекция простых. Проверяйте каждого кандидата на делимость. Тем болеее что найденых уже 90% и осталось "добить чутка".

[Wataru]

mayton2019, Для оставшихся 100000000 чисел решето будет быстрее, чем каждое число проверять на делимость всех известных.
В решете будет N/3 +N/5+N/7+... < N (log log N) операций. (N - длина отрезка)

В наивной проверке всех чисел на делители у вас будет как минимум N/logN * (sqrt(N)/log(sqrt(N)) операций - ведь для каждого простого числа в ответе придется, таки, проверить все найденные простые из списка. Это ассимптотически растет сильно быстрее оценки для решета.

[mayton2019]

Wataru, в диапазоне от 900 млн до мильярда находится примерно 4.5 миллиона простых.
Ну не знаю. Я-бы сделал бенчмарк. Тут по сути идет трейдофф. Так или иначе нам нужен
список простых делителей. А когда он ну нас есть - упрощается проверка последовательности.
Четные. И тройки можно учесть в цикле. Тоесть мы можем двигаться более длинными прыжками
чем просто последовательность. В этом смысле наивный брутфорс - более прост т.к. ничего в алгоритме
менять не надо. Мы просто его продолжаем. Для решета - нужно аллоцировать память.

[Wataru]

mayton2019, Решето тоже элементарно учитывает четные и делящиеся на 3.
Но да, наивный брутфорс более прост. Но он медленнее. Хотя в этом диапазоне за десяток секунд отработает.
А решето - за пол секунды.

Позже сегодня сделаю бенчмарк, но можно и на пальцах посчитать.

Там 3400 простых. Поверю вам наслово, что в интересующем промежутке 4.5млн простых. Тогда получается 4.5млн*3400 = 15.3 млрд делений. И это не оптимизируется никакими пропусками четных. Это только вычисления, чтобы доказать, что простые - простые.

В решете же будет 3400 делений/умножений, чтобы найти первое делящееся в промежутке, и ~3* 100 млн сложений (log log 32000 * длина промежутка), чтобы проскакать по таблице помечая делящееся. Уже получается в 40 раз меньше более простых операций.
Но решето будет бегать по 12.5mb памяти. Все равно - на порядок быстрее, ибо x40 - хороший запас производительности.

Ну и последнее, автор явно про решето спрашивает же.

[mayton2019]

Wataru, хорошо. Бенчмарк - так бенчмарк.

Стадион - так стадион.

[Wataru]

mayton2019,

Found 3401 small primes in 0.5ms 
Sieve: Found 4838319 primes in 0.296255s
Naive: Found 4838319 primes in 25.4289s

Решето примерно в 100 раз быстрее.

Никакие пропуски делящихся на 2,3 и 5 не помогут, потому что большую часть времени наивный метод доказывает, что простые числа ни на что не делятся.

Заметьте, там уже оптимизация - проверяются только простые до корня.

Вот код

#include <iostream>
#include <vector>
#include <chrono>
#include <cmath>

const int BE = 900000000;
const int N = 1000000000;


std::vector<int> small_primes;


bool DivisibleBySmallPrimes(int x) {
	for (int i = 0; i < (int)small_primes.size() && small_primes[i]*small_primes[i] <= x; ++i) {
		if (x % small_primes[i] == 0) return true;
	}
	return false;
}

void ComputeSmallPrimes() {
	auto start = std::chrono::steady_clock::now();
	int n = sqrt(N);
	for (int i = 2; i <= n; ++i) {
		if (!DivisibleBySmallPrimes(i)) {
			small_primes.push_back(i);
		}
	}
	auto end = std::chrono::steady_clock::now();
    std::chrono::duration<double> elapsed_seconds = end-start;
	std::cout << "Found " << small_primes.size() << " small primes in " << elapsed_seconds.count()*1000 << "ms\n";
}

void ComputeSieve() {
	auto start = std::chrono::steady_clock::now();
	std::vector<bool> sieve(N-BE+1);
	std::vector<int> primes;
	for (auto &x : small_primes) {
		int i = (BE-1) - (BE-1) % x + x;
		while (i <= N) {
			sieve[i-BE] = true;
			i += x;
		}
	}
	for (int i = BE; i <= N; ++i) {
		if(!sieve[i-BE]) {
			primes.push_back(i);
		}
	}
	auto end = std::chrono::steady_clock::now();
    std::chrono::duration<double> elapsed_seconds = end-start;
	std::cout << "Sieve: Found " << primes.size() << " primes in " << elapsed_seconds.count() << "s\n";
}

void ComputeNaive() {
	auto start = std::chrono::steady_clock::now();
	std::vector<int> primes;
	for (int i = BE; i <= N; ++i) {
		if (!DivisibleBySmallPrimes(i)) {
			primes.push_back(i);
		}
	}
	auto end = std::chrono::steady_clock::now();
    std::chrono::duration<double> elapsed_seconds = end-start;
	std::cout << "Naive: Found " << primes.size() << " primes in " << elapsed_seconds.count() << "s\n";
}


int main() {
	ComputeSmallPrimes();
	ComputeSieve();
	ComputeNaive();
	return 0;
}

[mayton2019]

Wataru, Шикарно. Я возьму ваш исходник с вашего позволения.

Вот интересно наступит-ли момент когда "решето" станет невыгодным?

[Wataru]

mayton2019, По ассимптотике решето всегда будет более выгодно. Я выше уже приводил наброски. Поэтому наивный подход может быть быстрее только для мелких чисел. Скажем, до сотни. Но там вообще без разницы, как считать.

Я вообще не могу придумать ни одного примера, когда решето было бы хуже наивного подхода для больших чисел. Казалось бы, из-за большего потребления памяти можно было бы его обломать, но это только если нужен очень большой промежуток чисел. И его в любом случае можно сегментировать - самое большое потребление памяти будет на простые числа до корня. Но это нужно и в наивном методе.

Решето может быть хуже других продвинутых алгоритмов. Например вероятностных тестов на простоту. Вот ими можно проверить, а простое ли вот это вот очень большое число. Решетом до него никак не добраться, потому что нужно очень много простых чисел до корня. Но и наивный метод там по той же причине работать не будет.

[mayton2019]

Wataru, А вот кстати интересно. Если убрать из Миллера-Рабина все случайные генераторы и просто сделать его детерминистичным - то тогда получается быстрый проверятель и который дает известные ошибки. А эти известные ошибки задать таблично как часть функици и получистя вполне себе годный проверятель для известных диапазонов. И протестировать его там хотя-бы до нескольких миллиардов.

[Wataru]

mayton2019, Задолбаетесь табличку считать и хранить. Ведь он интересен как раз для очень больших чисел (десятки знаков). Там этих исключений будет куча. И для их нахождения потребуется время до тепловой смерти вселенной. А до нескольких миллиардов вам и решето за секунду все подсчитает.

Answer 1

[Rsa97]

Решето Эратосфена всегда начинается с 2. Его принцип - на каждом шаге находим следующее невычеркнутое число, оно является простым. Для построения массива от 900000000 до 1000000000 необходимо вычеркнуть из этого массива все кратные простым числам от 2 до 500000000, а для этого надо эти простые найти.

Comments

[Wataru]

Только нужны простые числа не от 2 до 500000000, а от 3 до sqrt(1000000000). Во-первых, тут алгоритм уже рассматривает только нечетные числа. Во-вторых, любое состовное число имеет простой делитель не больше его корня, поэтому не надо рассматривать простые больше sqrt(1000000000).

Их можно найти отдельным решетом или другим алгоритмом.

Answer 2

[Wataru]

Раз уж я в коментариях к вопросу и так код привел (ради доказательства, что кто-то в интернете неправ), то приведу его и в качестве ответа. Код был на коленке ради быстрой проверки написан, поэтому в продакшн его так не пихайте. Глобальные переменные, константы - все это надо бы причесать.

const int BE = 900000000;
const int N = 1000000000;

std::vector<int> small_primes;

bool DivisibleBySmallPrimes(int x) {
	for (int i = 0; i < (int)small_primes.size() && small_primes[i]*small_primes[i] <= x; ++i) {
		if (x % small_primes[i] == 0) return true;
	}
	return false;
}

void ComputeSmallPrimes() {
	int n = sqrt(N);
	for (int i = 2; i <= n; ++i) {
		if (!DivisibleBySmallPrimes(i)) {
			small_primes.push_back(i);
		}
	}
}

void ComputeSieve() {
	ComputeSmallPrimes()
	std::vector<bool> sieve(N-BE+1);
	std::vector<int> primes;
	for (auto &x : small_primes) {
		int i = (BE-1) - (BE-1) % x + x;
		while (i <= N) {
			sieve[i-BE] = true;
			i += x;
		}
	}
	for (int i = BE; i <= N; ++i) {
		if(!sieve[i-BE]) {
			primes.push_back(i);
		}
	}
}

upd:
ComputeSmallPrimes можно переписать простым решетом - будет еще быстрее.