Найти следующий штрих с учетом всех предыдущих
Я пишу рекурсивный генератор бесконечных простых чисел, и я почти уверен, что смогу оптимизировать его лучше.
Прямо сейчас, кроме таблицы поиска первых дюжин простых чисел, каждый вызов рекурсивной функции получает список всех предыдущих простых чисел.
Поскольку это ленивый генератор, сейчас я просто отфильтровываю любое число по модулю 0 для любого из предыдущих простых чисел и получаю первый нефильтрованный результат. (При проверке я использую короткие замыкания, поэтому первый раз, когда предыдущий штрих делит текущее число, равномерно прерывается с этой информацией.)
Прямо сейчас моя производительность ухудшается при поиске 400-го простого числа (37 813). Я ищу способы использовать уникальный факт, что у меня есть список всех предыдущих простых чисел, и я только ищу следующие, чтобы улучшить мой алгоритм фильтрации. (Большая часть информации, которую я могу найти, предлагает не ленивые сита для поиска простых чисел при ограничении или способы поиска p n- го простого числа при заданном p n-1, а не оптимизации для поиска p n при заданных 2... p n-1 простых числах.)
Например, я знаю, что главное число должно находиться в диапазоне (p n-1 + 1)... (p n-1 + p n-2). Прямо сейчас я начинаю свою фильтрацию целых чисел с p n-1 + 2 (поскольку p n-1 + 1 может быть простым только для p n-1 = 2, который предварительно вычисляется). Но так как это ленивый генератор, знание терминальных границ диапазона (p n-1 + p n-2) не помогает мне ничего фильтровать.
Что я могу сделать для более эффективной фильтрации с учетом всех предыдущих простых чисел?
Пример кода
@doc """
Creates an infinite stream of prime numbers.
iex> Enum.take(primes, 5)
[2, 3, 5, 7, 11]
iex> Enum.take_while(primes, fn(n) -> n < 25 end)
[2, 3, 5, 7, 11, 13, 17, 19, 23]
"""
@spec primes :: Stream.t
def primes do
Stream.unfold( [], fn primes ->
next = next_prime(primes)
{ next, [next | primes] }
end )
end
defp next_prime([]), do: 2
defp next_prime([2 | _]), do: 3
defp next_prime([3 | _]), do: 5
defp next_prime([5 | _]), do: 7
# ... etc
defp next_prime(primes) do
start = Enum.first(primes) + 2
Enum.first(
Stream.drop_while(
Integer.stream(from: start, step: 2),
fn number ->
Enum.any?(primes, fn prime ->
rem(number, prime) == 0
end )
end
)
)
end
primes функция начинается с пустого массива, получает для него следующее простое число (2 вначале), а затем 1) испускает его из потока и 2) добавляет его в начало стека простых чисел, используемого в следующем вызове. (Я уверен, что этот стек является источником некоторого замедления.)
next_primes функция берет в этом стеке. Начиная с последнего известного простого числа +2, он создает бесконечный поток целых чисел и отбрасывает каждое целое число, которое делится равномерно на любое известное простое число в списке, а затем возвращает первое вхождение.
Полагаю, это что-то похожее на ленивое пошаговое сито Эратосфена.
Вы можете увидеть некоторые основные попытки оптимизации: я начинаю проверять с p n-1 +2 и перебираю четные числа.
Я попробовал более дословно просеять Эратосфена, просто пропустив Integer.stream через каждое вычисление и, найдя простое число, обернуть Integer.stream в новый Stream.drop_while, который отфильтровал только кратные значения этого простого числа. Но поскольку потоки реализованы как анонимные функции, это искалечило стек вызовов.
Стоит отметить, что я не предполагаю, что вам нужны все предыдущие простые числа для генерации следующего. Просто так получилось, что благодаря моей реализации.
3 ответа
Вам не нужны все предыдущие простые числа, достаточно только тех, которые находятся ниже квадратного корня из текущей производственной точки, при создании композитов из простых чисел с помощью алгоритма сита Эратосфена.
Это значительно снижает требования к памяти. Простые числа - это просто те нечетные числа, которых нет среди композитов.
Каждое простое число p создает цепочку из своих кратных, начиная с квадрата, перечисляемого с шагом 2p (потому что мы работаем только с нечетными числами). Эти мультипликаторы, каждый со своим значением шага, сохраняются в словаре, таким образом формируя приоритетную очередь. В этой очереди приоритетов присутствуют только простые числа до квадратного корня текущего кандидата (такое же требование к памяти, как и для сегментированного сита E.).
Символично, что SoE
P = {3,5,7,9,...} \ U {{p2, p2+ 2p, p2+ 4p, p2+ 6p,...} | п в П}
Каждое (нечетное) простое число порождает поток своих кратных чисел; когда все эти потоки объединены вместе, мы получаем все (нечетные) составные числа; и простые числа - все шансы без композитов (и 2).
В Python (возможно, читается как исполняемый псевдокод),
def postponed_sieve(): # postponed sieve, by Will Ness
yield 2; yield 3; yield 5; yield 7; # original code David Eppstein,
D = {} # ActiveState Recipe 2002
ps = (p for p in postponed_sieve()) # a separate Primes Supply:
p = ps.next() and ps.next() # (3) a Prime to add to dict
q = p*p # (9) when its sQuare is
c = 9 # the next Candidate
while True:
if c not in D: # not a multiple of any prime seen so far:
if c < q: yield c # a prime, or
else: # (c==q): # the next prime's square:
add(D,c + 2*p,2*p) # (9+6,6 : 15,21,27,33,...)
p=ps.next() # (5)
q=p*p # (25)
else: # 'c' is a composite:
s = D.pop(c) # step of increment
add(D,c + s,s) # next multiple, same step
c += 2 # next odd candidate
def add(D,x,s): # make no multiple keys in Dict
while x in D: x += s # increment by the given step
D[x] = s
Как только простое произведено, оно может быть забыто. Отдельное первичное предложение берется из отдельного вызова одного и того же генератора, рекурсивно, для поддержания словаря. И основной источник для этого взят от другого, также рекурсивно. Каждый из них должен поставляться только до квадратного корня от его производственной точки, поэтому в целом требуется очень мало генераторов, а их размеры асимптотически незначимы (sqrt( sqrt( N)), так далее).
Для любого числа k вам нужно только попробовать деление с простыми числами до и включая √k. Это потому, что любое простое число больше √k необходимо умножить на простое число меньше √k.
Доказательство:
√k * √k = k so (a+√k) * √k> k (для всех 0). Отсюда следует, что (a+√k) делит k тогда и только тогда, когда есть другой делитель, меньший, чем √k.
Это обычно используется для ускорения поиска простых чисел.
Я написал программу, которая генерирует простые числа по порядку, без ограничений, и использовал ее для суммирования первого миллиарда простых чисел в моем блоге. Алгоритм использует сегментированное сито Эратосфена; дополнительные простые числа просеивания рассчитываются для каждого сегмента, поэтому процесс может продолжаться бесконечно, пока у вас есть место для хранения простых чисел просеивания. Вот псевдокод:
function init(delta) # Sieve of Eratosthenes
m, ps, qs := 0, [], []
sieve := makeArray(2 * delta, True)
for p from 2 to delta
if sieve[p]
m := m + 1; ps.insert(p)
qs.insert(p + (p-1) / 2)
for i from p+p to n step p
sieve[i] := False
return m, ps, qs, sieve
function advance(m, ps, qs, sieve, bottom, delta)
for i from 0 to delta - 1
sieve[i] := True
for i from 0 to m - 1
qs[i] := (qs[i] - delta) % ps[i]
p := ps[0] + 2
while p * p <= bottom + 2 * delta
if isPrime(p) # trial division
m := m + 1; ps.insert(p)
qs.insert((p*p - bottom - 1) / 2)
p := p + 2
for i from 0 to m - 1
for j from qs[i] to delta step ps[i]
sieve[j] := False
return m, ps, qs, sieve
Вот ps список простых чисел просеивания меньше текущего максимума и qs смещение наименьшего кратного соответствующего ps в текущем сегменте. advance функция очищает битаррай, сбрасывает qs, расширяет ps а также qs с новыми просеивающими простыми числами, затем просеивает следующий сегмент.
function genPrimes()
bottom, i, delta := 0, 1, 50000
m, ps, qs, sieve := init(delta)
yield 2
while True
if i == delta # reset for next segment
i, bottom := -1, bottom + 2 * delta
m, ps, qs, sieve := \textbackslash
advance(m, ps, qs, sieve, bottom, delta)
else if sieve[i] # found prime
yield bottom + 2*i + 1
i := i + 1
Размер сегмента 2 * delta произвольно установлен в 100000. Этот метод требует O(sqrt(n)) место для простых чисел просеивания плюс постоянное пространство для решета.
Это медленнее, но экономит место для генерации кандидатов с колесом и проверки кандидатов на первичность.
function genPrimes()
w, wheel := 0, [1,2,2,4,2,4,2,4,6,2,6,4,2,4, \
6,6,2,6,4,2,6,4,6,8,4,2,4,2,4,8,6,4,6, \
2,4,6,2,6,6,4,2,4,6,2,6,4,2,4,2,10,2,10]
p := 2; yield p
repeat
p := p + wheel[w]
if w == 51 then w := 4 else w := w + 1
if isPrime(p) yield p
Может быть полезно начать с сита и переключиться на колесо, когда сито становится слишком большим. Еще лучше продолжать просеивать с некоторым фиксированным набором простых чисел, когда набор становится слишком большим, а затем сообщать только об этих значениях. bottom + 2*i + 1 которые проходят тест на первичность.