Эффективное сито эратосфена в питоне
Этот очень короткий и простой код в #Python пытается смоделировать "Сито Эратосфена" для первых N натуральных чисел с ограничениями (0) краткости сценария; (1) сведение к минимуму циклов if и for/while; (2) эффективность с точки зрения процессорного времени.
import numpy as np
N = 10**5
a = np.array(range(3,N,2))
for j in range(0, int(round(np.sqrt(N),0))):
a[(a!=a[j]) & (a%a[j] == 0)] = 0
a = a[a!=0]
a = [2]+list(a)
На Intel Core I5 он возвращает простые числа среди первых:
- N = 100000 за 0,03 секунды;
- N = 1 000 000 за 0,63 секунды;
- N = 10 000 000 за 22,2 секунды.
Кто-то хотел бы поделиться более эффективными кодами с точки зрения процессорного времени в рамках вышеупомянутых ограничений?
1 ответ
Фактическое сито Эратосфена NumPy выглядит так:
def sieve(n):
flags = numpy.ones(n, dtype=bool)
flags[0] = flags[1] = False
for i in range(2, n):
# We could use a lower upper bound for this loop, but I don't want to bother with
# getting the rounding right on the sqrt handling.
if flags[i]:
flags[i*i::i] = False
return numpy.flatnonzero(flags)
Он поддерживает массив "возможно простых" флагов и напрямую сбрасывает флаги, соответствующие кратным числам простых чисел, без необходимости проверять делимость, особенно для чисел, которые не делятся на простое число, обрабатываемое в настоящее время.
То, что вы делаете, это пробное деление, где вы просто проходите и проверяете, делятся ли числа делителями-кандидатами. Даже при хорошей реализации пробного деления необходимо выполнять больше операций и более дорогостоящих операций, чем сито. Ваша реализация делает еще больше работы, чем эта, потому что она рассматривает делители-кандидаты, не являющиеся простыми числами, и потому, что она продолжает выполнять тесты делимости для чисел, которые она уже должна знать, являются простыми.
Я решил поиграть с этим и создал дополнительную оптимизированную версию NumPy, опубликованную @user2357112, которая поддерживает Monica, которая использует Numba JIT для ускорения работы.
import numba
import numpy
import timeit
import datetime
@numba.jit(nopython = True, parallel = True, fastmath = True, forceobj = False)
def sieve (n: int) -> numpy.ndarray:
primes = numpy.full(n, True)
primes[0], primes[1] = False, False
for i in numba.prange(2, int(numpy.sqrt(n) + 1)):
if primes[i]:
primes[i*i::i] = False
return numpy.flatnonzero(primes)
if __name__ == "__main__":
timestart = timeit.default_timer()
print(sieve(1000000000))
timestop = timeit.default_timer()
timedelta = (timestop - timestart)
print(f"Time Elapsed: {datetime.timedelta(seconds = timedelta)}")
else:
pass
На своем ноутбуке я отсеиваю простые числа в первом миллионе (1e9
) натуральные числа в
0:00:10.378686
секунд. JIT обеспечивает здесь как минимум порядок производительности; следующий самый быстрый ответ на момент написания взял
0:01:27.059963
минут. К сожалению, у меня в этой системе (Mac) нет графического процессора Nvidia и Cuda, иначе я бы использовал это.
1.94 секунды для 10.000.000
def sieve_eratosthene(limit):
primes = [True] * (limit+1)
iter = 0
while iter < limit**0.5 :
if iter < 2:
primes[iter]= False
elif primes[iter]:
for i in range(iter*2, limit+1, iter):
primes[i] = False
iter+=1
return(x for x in range(number+1) if primes[x])
Вот простой способ сделать это с помощью NumPy. Это несколько похоже на идею OP об индексировании вместо повторного цикла внутри основного цикла, но без проверки деления, только нарезка.
Он также похож на то, что user2357112 поддерживает ответ Моники, но он учитывает только нечетные числа, что делает его быстрее.
Это типичное сито для нечетных чисел: /questions/26957514/sokraschenie-ispolzovaniya-pamyati-pri-proektirovanii-sita-iz-eratosfena-v-s/26957523#26957523.
- Создайте массив bool размера n.
- Перебирайте нечетные числа до квадратного корня из n.
Примечание: числа используются как индексы как взаимозаменяемые - Пометить составные числа (кратные простым числам) как ложные.
В конце концов, у нас есть массив bool, который мы можем использовать, чтобы проверить, является ли число ниже n простым (кроме четных чисел, вы можете проверить их без массива, используя & 1 или около того)
Средн. время для n = 20000000: 0,1063 с
import numpy as np
n = 20000000
isprime = np.ones(n, dtype=np.bool)
# odd only sieve
i = 3
while (i * i < n):
if isprime[i]:
isprime[i * i:n:2 * i] = False
i += 2
# test
print(isprime[97]) # should be True