Введение
Одним из алгоритмов для поиска простых чисел является Решето Эратосфена предложенное еще древнегреческим математиком.
Картинка из википедии:
Смысл в вычеркивании чисел кратных уже найденным простым. Остающиеся невычеркнутыми в свою очередь являются простыми. Более подробно расписано тут.
Одна из проблем при поиске решетом это объем памяти который надо выделить под фильтруемые числа. Вычеркнутые непростые удаляются уменьшая память, но изначально объем требуется большой.
Для решения используется сегментация (когда память выделяется по кускам) и другие ухищрения (см. тут).
Реализация алгоритма
Алгоритм внизу (написан на java) предполагает минимальный объем памяти — по сути для каждого найденного простого числа мы храним еще одно число — последнее зачеркнутое (наибольшее). Если я правильно оцениваю объем памяти ln(n) — число найденных простых.
Суть алгоритма:
Допустим мы имеем несколько уже найденных простых чисел отсортированных по возрастанию. (Алгоритм стартует с [2,3]). Для каждого из них храним последнее (наибольшее) зачеркнутое число. Инициализируем самими простыми числами.
Выбираем кандидата в простые например наибольшее найденное простое число +1 (в алгоритме внизу перескакиваем четные как заведомо не простые).
Кандидат сравнивается с последним зачеркнутым текущего простого. Пока текущее зачеркнутое меньше кандидата увеличиваем его на текущее простое. Если текущее зачеркнутое равно кандидату, то кандидат не простое число. Выбираем следующего кандидата.
В случае если текущее зачеркнутое больше кандидата, проверяем кандидата на следующем простом числе.
Если добрались до конца списка простых чисел (то есть все зачеркнутые больше кандидата) мы нашли очередное простое число.
Добавляем его в список и инициализируем последнее зачеркнутое найденным простым.
Код на java
import java.util.ArrayList;
import java.util.List;
public class SieveEratosthenes {
static class PrimePair {
Integer prime;
Integer lastCrossed;
PrimePair(Integer prime, Integer lastCrossed) {
this.prime = prime;
this.lastCrossed = lastCrossed;
}
}
private List<PrimePair> primes;
private SieveEratosthenes() {
primes = new ArrayList<>();
primes.add(new PrimePair(2, 2));
primes.add(new PrimePair(3, 3));
}
private void fillNPrimes(int n) {
while (primes.size()<n) {
addNextPrime();
}
}
private void addNextPrime() {
int candidate = primes.get(primes.size()-1).prime + 2;
for (int i = 1; i < primes.size(); i++) {
PrimePair p = primes.get(i);
while (p.lastCrossed < candidate) {
p.lastCrossed += p.prime;
}
if (p.lastCrossed == candidate) {
//restart
candidate+=2;
i=-1;
}
}
System.out.println(candidate);
primes.add(new PrimePair(candidate, candidate));
}
public static void main(String[] args) {
SieveEratosthenes test = new SieveEratosthenes();
test.fillNPrimes(1000);
}
}
Тот же код на питоне:
primes = [2, 3]
last_crossed = [2, 3]
def add_next_prime():
candidate = primes[-1] + 2
i = 0
while i < len(primes):
while last_crossed[i] < candidate:
last_crossed[i] += primes[i]
if last_crossed[i] == candidate:
candidate += 2
i = 0
i += 1
primes.append(candidate)
last_crossed.append(candidate)
def fill_primes(n):
while len(primes) < n:
add_next_prime()
fill_primes(1000)
print(primes)
Вся логика традиционных алгоритмов с колесами может быть применена при выборе следующего кандидата.
> GitHub
Возможно это очередное изобретение велосипеда. Готов выслушать замечания.
Комментарии (30)
dom1n1k
14.07.2017 17:52+3Писал своего Эратосфена году эдак в 2002.
На машине было то ли 128, то ли 256 мб памяти (не помню точно, до апгрейда это всё происходило или после), а хотелось проверять 4 миллиарда чисел — ну то есть максимум, что в int32 влезает.
Я поступил просто: флажок для каждого нечетного числа соответствовал одному биту, четные пропускались вообще. 4 млрд чисел / 16 = 256 мб памяти.
Внутри цикла индекс преобразовывался в битовую маску, проставлялись флажки тоже битовыми операциями.
Писалось всё на C++.
Время работы, если склероз не подводит, было что-то в районе 5 минут на Celeron-667.
Код за давностью лет, к сожалению, утерян.
StanislavL
17.07.2017 08:59Есть Например тут http://primesieve.org/ очень эффективное решето. Цель не написать чтобы найти сколько то чисел. Цель скорее показать идею.
AnROm
14.07.2017 19:02Как по мне, в оптимальной реализации решето Эратосфена есть только одна хитрость, связанная с уменьшением количества проходов по последовательности чисел:
пусть N — максимальное число, то для вычеркивания чисел из последовательности нужно рассматривать все числа от 2 до sqrt(N). Для всех оставшихся чисел, больших sqrt(N), все непростые числа к этому моменту будут вычеркнуты. Естественно четные числа нужно пропускать.
N = 12
sqrt(12) = 3
После вычеркивания из последовательности всех чисел кратных двум и трем останутся только простые.Останется только по битовым маскам пройти, чтоб вывести невычеркнутые числа.
daiver19
15.07.2017 00:24+1Если я правильно оцениваю объем памяти ln(n) — число найденных простых.
Неправильно. Количество простых чисел асимптотически растет как n/ln(n). Вы б сами себя проверили, ln (1000) всего лишь 7.
Ну и да, выходит что ваш алгоритм работает за prime_count^2 или я что-то упускаю? Какое это тогда решето Эратосфена?StanislavL
17.07.2017 10:35Да. Вы правы. По памяти n/ln(n),
Но насчет prime_count^2 нет.
для 2 мы сделаем не более n/2 сложений, для 3 — n/3 и т.д.
получаем n/2 + n/3 + n/5… то есть n * log log(n)
а если 2 выкинуть четные незачем проверять. И вот чуть ниже есть коммент что стартовать можно с candidate^2 и перепрыгивать при суммировании
if (p.lastCrossed < candidate) {
p.lastCrossed += p.prime * Math.round((candidate — p.lastCrossed) / p.prime());
}
то сложность поменьше будет. Хотя так сходу не могу прикинуть число. Буду думать.
mikhaelkh
15.07.2017 01:54+1StanislavL
Мда… Попытка минимизировать память и алгоритм запускается для 1000…
Ваше "решето" использует N^(1+o(1)) памяти и N^(2+o(1)) времени.
В Википедии написано гораздо лучше: Решето Эратосфена, Решето АткинаStanislavL
17.07.2017 09:03Цель же не найти первых 1000 или 100 млн. Цель показать идею алгоритма.
mikhaelkh
18.07.2017 13:21Приведите пример задачи, где ваша идея будет лучше обычного решета с битовым массивом. Ваша идея даже памяти больше расходует — ~N / ln(N) * 2 * log2(N) = 2.885 * N бит против ~N бит.
StanislavL
18.07.2017 14:42Памяти используется n/ln(n) — число простых. очевидно лучше n. Откуда вы взяли * 2 * log2(N)?
По операциям n log(log(n)) (n/2 + n/3 + n/5+… n/Pk… +)mikhaelkh
23.07.2017 03:57Длина простого числа порядка n равна log2(n) бит и у Вас 2 массива. На каждом простом Вы просматриваете весь список простых — сложность (n/ln(n))^2. Запустите свой код хотя бы для миллиона.
realbtr
15.07.2017 03:00Память можно экономить, выбирая только числа вида 6n+-1, далее 30n+-1,7,11,13, далее 210n+-1,11,13,17,19,23,29,… потом 2310n+-1,13,17,… и так постепено от 1/2 (нечетные) мы переходим к 1/3, 4/15 и т.д.
StanislavL
17.07.2017 09:04Wheel factorization и т.п. по прежнему применимо при выборе следующего кандидата.
xdemonx
17.07.2017 08:46Для ещё доп оптимизации можно заменить код
primes.add(new PrimePair(candidate, candidate));
на
primes.add(new PrimePair(candidate, candidate * candidate));
т.к. меньшие непростые числа будут отсеяны предыдущими простыми числамиStanislavL
17.07.2017 08:56Наверное candidate*2 в квадрате будет многовато.
vesper-bot
17.07.2017 09:58Вообще, теорема есть, что если из решета Эратосфена выкинуты все числа, кратные первым N простым, следующее оставшееся составное число будет равно квадрату N+1-го простого числа. Доказывается довольно просто, поэтому подход правильный.
speshuric
17.07.2017 08:461. Это не решето Эратосфена, а перебор делителей (trial division)
2. В решете Эратосфена без каких-либо ухищрений, используя банальнейший BitSet можно положить все положительные int 0..2^31-1 в 256 мегабайт. При этом в этом объёме примерно 105 млн простых (спасибо www.wolframalpha.com за подсказку). Т.е. при хранении в массиве это будет примерно 400 мегабайт (если хранить компактно).
3. Если из BitSet выкинуть даже биты, соответствующие составным числам в интервале 1..30 т. е. считать, что остались только 1, 7, 11, 13, 17, 19, 23, 29 по модулю 30, то это усложнит арифметику, но позволит все простые intы уместить в 70 МБ. При переходе в long разница становится еще более существенной.StanislavL
17.07.2017 09:061. Где же там делители? попыток деления не происходит.
2. Если надо искать тысячебитовые числа куда их складывать?speshuric
18.07.2017 00:59Хм… извиняюсь, да, невнимателен был. Это не trial division, но и решетом Эратосфена это сложно назвать.
Но к памяти всё равно много вопросов.
В этом варианте на каждое простое используется пара простого и границы интервала. Размер такой структуры, конечно, зависит от реализации, но это никак не может быть меньше двух целых чисел. Если у вас целые всего 32 бита, то это 64 бита. Для того, чтобы это было выгодным по сравнению с BitSet (или подобной структурой), надо чтобы в среднем между простыми было 64 составных. А в целых этого диапазона примерно каждое 20-е простое.
С 1000-битными на самом деле для вашей модели не лучше. Количество требуемых бит для числа растет быстрее, чем средний интервал между простыми примерно в 1/ln(2) раз. Так что идея интересная, но не работает.
Ну и да, строить решето Эратосфена для 1000-битных чисел — нашей вселенной не хватит (ни места ни времени).
Если надо искать 1000-битные простые числа не в теории, а на практике, то в Java всё уже до нас придумано и почти всё сделано. У класса BigInteger есть метод probablePrime, а его уже можно проверить каким-либо из подходящих под задачу primality tests (в BigInteger только вероятностный, так что, увы — придётся писать).StanislavL
18.07.2017 09:06По памяти это, как уже правильно заметили в комментах, n / ln(n) — по факту мы держим количество простых и для каждого из них последнее зачеркнутое.
BigInteger может быть. java выбрана потому как для меня более знакомый язык. Корректнее бы на С писать выжимая все возможное побитово.
Вероятностный мало. Надо доказывать что это простое, а проверок будет слишком много.
В общем компьютер должен работать, а человек думать. Буду думать дальше. Спасибо за коммент.
ijusti
19.07.2017 09:00Чтобы уменьшить память надо использовать примитивы в PrimePair, как следствие не будет ещё и автобоксинга. Цикл лучше использовать через итератор и делать next и previous, вместо get(i)
StanislavL
19.07.2017 09:01Это не принципиально. Для реальной скорости надо использовать массивы, а не inner class.
Я знаю как ускорять java. Так написано чтобы было понятнее как работает алгоритм.
realbtr
24.07.2017 20:30Я рассматривал как вариант хранить только исключения вида i*j*k*..., рекурсивно вычисляя сами числа грейдами и не сохраняя простые числа как-то иначе.
Как я уже писал (и мне ответили про wheel factirization, обязательно изучу, о чем это), можно развивать простые числа по принципу
2n+1 (дает числа 3 и 5),
6n+{1,5} (дает числа 7, 11, 13, 17, 19, 23, 25, 29} и нужно сохранить лишь одно число пока 5*5, как исключение.
30n+{1,7,11,13,17,19,23,29}, которое дает все простые числа до 209, плюс произведения i*j, где 7<=i,j<=sqrt(210)
Далее будет база 2*3*5*7=210 простых чисел до 2309 (2*3*5*7*11-1) и нужно исключать произведения i*j*k*… где 11<=i,j,k,..<=sqrt(2310)
Данная рекурсия может быть вообще чисто вычислительной, не используя память для хранения всех чисел (например при проверке числа на простоту), а для вычисления последовательных чисел можно заметить просто правило, 6n+{1,5} = 6m+-1, 30n+{1,7,11,13,17,19,23,29}= 30m+-{1,7,11,13} и т.д.то есть — зеркальная симметрия каждого такого блока. И блоки большего размера фрактально повторяются из блоков меньшего размера.
Я писал как-то статью на эту тему здесь https://habrahabr.ru/post/255527/ но видимо открывал велосипед и получил немало критических замечаний.
ibessonov
Смущает данный код:
Чем он лучше такой банальной проверки?
vesper-bot
Например, для простого p=65537 верхний код будет выполнен один раз на 65537/2 итераций, а проверка условия в нижнем — каждую итерацию, при отсутствии делимости на остальной список простых. Код написан, чтобы ещё и минимизировать количество делений, которые для больших p будут выполняться ненужное число раз.
profesor08
StanislavL
Да, это улучшение. Там еще много чего есть поулучшить. Я пытался просто идею показать…