Ускорить bitstring/bit операции в Python?
Я написал генератор простых чисел, используя Сито Эратосфена и Python 3.1. Код работает корректно и грациозно на 0,32 секунды на ideone.com для генерации простых чисел до 1 000 000.
# from bitstring import BitString
def prime_numbers(limit=1000000):
'''Prime number generator. Yields the series
2, 3, 5, 7, 11, 13, 17, 19, 23, 29 ...
using Sieve of Eratosthenes.
'''
yield 2
sub_limit = int(limit**0.5)
flags = [False, False] + [True] * (limit - 2)
# flags = BitString(limit)
# Step through all the odd numbers
for i in range(3, limit, 2):
if flags[i] is False:
# if flags[i] is True:
continue
yield i
# Exclude further multiples of the current prime number
if i <= sub_limit:
for j in range(i*3, limit, i<<1):
flags[j] = False
# flags[j] = True
Проблема заключается в том, что у меня заканчивается память, когда я пытаюсь создать числа до 1 000 000 000.
flags = [False, False] + [True] * (limit - 2)
MemoryError
Как вы можете себе представить, выделение 1 миллиарда булевых значений ( 1 байт 4 или 8 байтов (см. комментарий) каждый в Python) действительно невозможно, поэтому я просмотрел bitstring. Я полагал, что использование 1 бит для каждого флага будет намного более эффективным для памяти. Однако производительность программы резко снизилась - 24 секунды в режиме ожидания, для простого номера до 1 000 000. Вероятно, это связано с внутренней реализацией bitstring.
Вы можете комментировать/раскомментировать три строки, чтобы увидеть, что я изменил, чтобы использовать BitString в качестве фрагмента кода выше.
Мой вопрос: есть способ ускорить мою программу, с или без bitstring?
Изменить: проверьте исходный код перед отправкой. Я не могу принимать ответы, которые работают медленнее, чем мой существующий код, естественно.
Изменить еще раз:
Я собрал список тестов на своей машине.
Ответы
Ответ 1
Есть несколько небольших оптимизаций для вашей версии. Перевернув роли True и False, вы можете изменить "if flags[i] is False:" на "if flags[i]:". Начальное значение для второго оператора range может быть i*i вместо i*3. Ваша оригинальная версия занимает 0.166 секунды в моей системе. С этими изменениями версия ниже занимает 0.156 секунд в моей системе.
def prime_numbers(limit=1000000):
'''Prime number generator. Yields the series
2, 3, 5, 7, 11, 13, 17, 19, 23, 29 ...
using Sieve of Eratosthenes.
'''
yield 2
sub_limit = int(limit**0.5)
flags = [True, True] + [False] * (limit - 2)
# Step through all the odd numbers
for i in range(3, limit, 2):
if flags[i]:
continue
yield i
# Exclude further multiples of the current prime number
if i <= sub_limit:
for j in range(i*i, limit, i<<1):
flags[j] = True
Это не поможет вам решить проблему с памятью.
Перейдя в мир расширений C, я использовал версию разработки gmpy. (Отказ от ответственности: я один из разработчиков). Версия разработки называется gmpy2 и поддерживает изменяемые целые числа, называемые xmpz. Используя gmpy2 и следующий код, у меня есть время работы 0.140 секунд. Продолжительность работы для лимита в 1 000 000 000 составляет 158 секунд.
import gmpy2
def prime_numbers(limit=1000000):
'''Prime number generator. Yields the series
2, 3, 5, 7, 11, 13, 17, 19, 23, 29 ...
using Sieve of Eratosthenes.
'''
yield 2
sub_limit = int(limit**0.5)
# Actual number is 2*bit_position + 1.
oddnums = gmpy2.xmpz(1)
current = 0
while True:
current += 1
current = oddnums.bit_scan0(current)
prime = 2 * current + 1
if prime > limit:
break
yield prime
# Exclude further multiples of the current prime number
if prime <= sub_limit:
for j in range(2*current*(current+1), limit>>1, prime):
oddnums.bit_set(j)
Нажав оптимизацию и жертвуя ясностью, я получаю время работы от 0,107 до 123 секунд со следующим кодом:
import gmpy2
def prime_numbers(limit=1000000):
'''Prime number generator. Yields the series
2, 3, 5, 7, 11, 13, 17, 19, 23, 29 ...
using Sieve of Eratosthenes.
'''
yield 2
sub_limit = int(limit**0.5)
# Actual number is 2*bit_position + 1.
oddnums = gmpy2.xmpz(1)
f_set = oddnums.bit_set
f_scan0 = oddnums.bit_scan0
current = 0
while True:
current += 1
current = f_scan0(current)
prime = 2 * current + 1
if prime > limit:
break
yield prime
# Exclude further multiples of the current prime number
if prime <= sub_limit:
list(map(f_set,range(2*current*(current+1), limit>>1, prime)))
Изменить: на основе этого упражнения я изменил gmpy2, чтобы принять xmpz.bit_set(iterator). Используя следующий код, время выполнения для всех простых чисел менее 1 000 000 000 составляет 56 секунд для Python 2.7 и 74 секунды для Python 3.2. (Как отмечено в комментариях, xrange быстрее, чем range.)
import gmpy2
try:
range = xrange
except NameError:
pass
def prime_numbers(limit=1000000):
'''Prime number generator. Yields the series
2, 3, 5, 7, 11, 13, 17, 19, 23, 29 ...
using Sieve of Eratosthenes.
'''
yield 2
sub_limit = int(limit**0.5)
oddnums = gmpy2.xmpz(1)
f_scan0 = oddnums.bit_scan0
current = 0
while True:
current += 1
current = f_scan0(current)
prime = 2 * current + 1
if prime > limit:
break
yield prime
if prime <= sub_limit:
oddnums.bit_set(iter(range(2*current*(current+1), limit>>1, prime)))
Редактировать # 2: Еще одна попытка! Я изменил gmpy2, чтобы принять xmpz.bit_set(slice). Используя следующий код, время запуска для всех простых чисел менее 1,000,000,000 составляет около 40 секунд для Python 2.7 и Python 3.2.
from __future__ import print_function
import time
import gmpy2
def prime_numbers(limit=1000000):
'''Prime number generator. Yields the series
2, 3, 5, 7, 11, 13, 17, 19, 23, 29 ...
using Sieve of Eratosthenes.
'''
yield 2
sub_limit = int(limit**0.5)
flags = gmpy2.xmpz(1)
# pre-allocate the total length
flags.bit_set((limit>>1)+1)
f_scan0 = flags.bit_scan0
current = 0
while True:
current += 1
current = f_scan0(current)
prime = 2 * current + 1
if prime > limit:
break
yield prime
if prime <= sub_limit:
flags.bit_set(slice(2*current*(current+1), limit>>1, prime))
start = time.time()
result = list(prime_numbers(1000000000))
print(time.time() - start)
Изменить №3: я обновил gmpy2, чтобы правильно поддерживать нарезку на уровне бит xmpz. Никаких изменений в производительности, но очень хороший API. Я немного подкорректировал, и у меня есть время примерно до 37 секунд. (См. Редактирование № 4 на изменения в gmpy2 2.0.0b1.)
from __future__ import print_function
import time
import gmpy2
def prime_numbers(limit=1000000):
'''Prime number generator. Yields the series
2, 3, 5, 7, 11, 13, 17, 19, 23, 29 ...
using Sieve of Eratosthenes.
'''
sub_limit = int(limit**0.5)
flags = gmpy2.xmpz(1)
flags[(limit>>1)+1] = True
f_scan0 = flags.bit_scan0
current = 0
prime = 2
while prime <= sub_limit:
yield prime
current += 1
current = f_scan0(current)
prime = 2 * current + 1
flags[2*current*(current+1):limit>>1:prime] = True
while prime <= limit:
yield prime
current += 1
current = f_scan0(current)
prime = 2 * current + 1
start = time.time()
result = list(prime_numbers(1000000000))
print(time.time() - start)
Изменить # 4: Я внес некоторые изменения в gmpy2 2.0.0b1, которые нарушают предыдущий пример. gmpy2 больше не обрабатывает True как специальное значение, которое обеспечивает бесконечный источник 1-бит. -1 следует использовать вместо этого.
from __future__ import print_function
import time
import gmpy2
def prime_numbers(limit=1000000):
'''Prime number generator. Yields the series
2, 3, 5, 7, 11, 13, 17, 19, 23, 29 ...
using Sieve of Eratosthenes.
'''
sub_limit = int(limit**0.5)
flags = gmpy2.xmpz(1)
flags[(limit>>1)+1] = 1
f_scan0 = flags.bit_scan0
current = 0
prime = 2
while prime <= sub_limit:
yield prime
current += 1
current = f_scan0(current)
prime = 2 * current + 1
flags[2*current*(current+1):limit>>1:prime] = -1
while prime <= limit:
yield prime
current += 1
current = f_scan0(current)
prime = 2 * current + 1
start = time.time()
result = list(prime_numbers(1000000000))
print(time.time() - start)
Изменить # 5: Я сделал некоторые улучшения для gmpy2 2.0.0b2. Теперь вы можете перебирать все биты, которые либо установлены, либо очищены. Время работы улучшилось на ~ 30%.
from __future__ import print_function
import time
import gmpy2
def sieve(limit=1000000):
'''Returns a generator that yields the prime numbers up to limit.'''
# Increment by 1 to account for the fact that slices do not include
# the last index value but we do want to include the last value for
# calculating a list of primes.
sieve_limit = gmpy2.isqrt(limit) + 1
limit += 1
# Mark bit positions 0 and 1 as not prime.
bitmap = gmpy2.xmpz(3)
# Process 2 separately. This allows us to use p+p for the step size
# when sieving the remaining primes.
bitmap[4 : limit : 2] = -1
# Sieve the remaining primes.
for p in bitmap.iter_clear(3, sieve_limit):
bitmap[p*p : limit : p+p] = -1
return bitmap.iter_clear(2, limit)
if __name__ == "__main__":
start = time.time()
result = list(sieve(1000000000))
print(time.time() - start)
print(len(result))
Ответ 2
Хорошо, так что это мой второй ответ, но поскольку скорость имеет смысл, я подумал, что должен был упомянуть bitarray модуль - даже если он bitstring nemesis:). Он идеально подходит для этого случая, поскольку он не только является расширением C (и, следовательно, быстрее, чем чистый Python имеет надежду на то), но также поддерживает назначение срезов. Это еще не доступно для Python 3.
Я даже не пытался оптимизировать это, я просто переписал версию bitstring. На моей машине я получаю 0,16 секунды для простых чисел менее миллиона.
За миллиард он работает отлично и завершается через 2 минуты 31 секунду.
import bitarray
def prime_bitarray(limit=1000000):
yield 2
flags = bitarray.bitarray(limit)
flags.setall(False)
sub_limit = int(limit**0.5)
for i in xrange(3, limit, 2):
if not flags[i]:
yield i
if i <= sub_limit:
flags[3*i:limit:i*2] = True
Ответ 3
Хорошо, здесь (почти полный) всеобъемлющий бенчмаркинг, который я сделал сегодня вечером, чтобы узнать, какой код работает быстрее всего. Надеюсь, кто-то найдет этот список полезным. Я пропустил все, что занимает больше 30 секунд, чтобы завершить работу на моей машине.
Я хотел бы поблагодарить всех, кто внес вклад. Я получил много понимания от ваших усилий, и я надеюсь, что вы тоже.
Моя машина: AMD ZM-86, 2.40 Ghz Dual-Core, с 4 ГБ оперативной памяти. Это ноутбук HP Touchsmart Tx2. Обратите внимание, что хотя я, возможно, связался с пастебином, , я сравнивал все следующие на моей машине.
Я добавлю тест gmpy2, когда смогу его создать.
Все тесты тестируются в Python 2.6 x86
Возврат списка простых чисел n до 1 000 000: (с использованием Python генераторы)
версия генератора numbase Sebastian (обновлена) - 121 мс @
Mark Sieve + Wheel - 154 мс
версия Robert с разрезом - 159 мс
Моя улучшенная версия с нарезкой- 205 мс
Генератор Numpy с перечислением - 249 мс @
Mark Basic Sieve - 317 мс
улучшение casevh на моем оригинале решение - 343 мс
Модифицированное решение генератора numpy - 407 мс
Мой оригинальный метод в вопрос - 409 мс
Bitarray Solution - 414 мс @
Pure Python с bytearray - 1394 мс @
Скотта BitString решение - 6659 ms @
'@' означает, что этот способ способен генерировать до n < 1 000 000 000 моя установка машины.
Кроме того, если вам не нужны генератор и просто хотите весь список сразу:
решение numpy от RosettaCode - 32 мс @
(Решение numpy также способно генерировать до 1 миллиарда, что заняло 61,6259 секунд. Я подозреваю, что память была заменена один раз, следовательно, в два раза.)
Ответ 4
Связанный вопрос.
Самый быстрый способ перечислить все простые числа ниже N в python
Привет, я тоже ищу код в python для генерации простых чисел до 10 9, насколько это возможно, что сложно из-за проблемы с памятью. до сих пор я придумал это для генерации простых чисел до 10 6 и 10 * 7 (0,171s и 1,764s соответственно на моем старом компьютере), который, кажется, немного быстрее (по крайней мере, на моем компьютере), чем "Моя улучшенная версия с разрезом" из предыдущего сообщения, вероятно, потому, что я использую n//ii +1 (см. Ниже) вместо аналогичной len (флаги [i2:: я < 1]) в другом коде. пожалуйста, исправьте меня, если я ошибаюсь. Любые предложения по улучшению очень приветствуются.
def primes(n):
""" Returns a list of primes < n """
sieve = [True] * n
for i in xrange(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 xrange(3,n,2) if sieve[i]]
В одном из своих кодов Xavier использует флаги [i2:: я < 1 и len (flags [i2:: я < 1]), я вычислил размер явным, используя тот факт, что между ii.. n мы имеем (n-ii)//2 * я кратные 2 * i, так как мы хотим сосчитать ii, и мы суммируем 1 так, что len (sieve [ii:: 2 * i]) равно (ni * i)//( 2 * i) +1 Это сделает код быстрее. Базовый генератор, основанный на коде выше, будет:
def primesgen(n):
""" Generates all primes <= n """
sieve = [True] * n
yield 2
for i in xrange(3,int(n**0.5)+1,2):
if sieve[i]:
yield i
sieve[i*i::2*i] = [False]*((n-i*i-1)/(2*i)+1)
for i in xrange(i+2,n,2):
if sieve[i]: yield i
с небольшим изменением можно написать несколько более медленную версию вышеприведенного кода, которая начинается с сита половину размера сита = [True] * (n//2) и работает для того же n. не знаете, как это отражается в памяти. В качестве примера реализации
измененная версия кода numset rosetta (возможно, быстрее), начиная с половины размера сита.
import numpy
def primesfrom3to(n):
""" Returns a array of primes, 3 <= p < n """
sieve = numpy.ones(n/2, dtype=numpy.bool)
for i in xrange(3,int(n**0.5)+1,2):
if sieve[i/2]: sieve[i*i/2::i] = False
return 2*numpy.nonzero(sieve)[0][1::]+1
Быстрее и больше генератора памяти:
import numpy
def primesgen1(n):
""" Input n>=6, Generates all primes < n """
sieve1 = numpy.ones(n/6+1, dtype=numpy.bool)
sieve5 = numpy.ones(n/6 , dtype=numpy.bool)
sieve1[0] = False
yield 2; yield 3;
for i in xrange(int(n**0.5)/6+1):
if sieve1[i]:
k=6*i+1; yield k;
sieve1[ ((k*k)/6) ::k] = False
sieve5[(k*k+4*k)/6::k] = False
if sieve5[i]:
k=6*i+5; yield k;
sieve1[ ((k*k)/6) ::k] = False
sieve5[(k*k+2*k)/6::k] = False
for i in xrange(i+1,n/6):
if sieve1[i]: yield 6*i+1
if sieve5[i]: yield 6*i+5
if n%6 > 1:
if sieve1[i+1]: yield 6*(i+1)+1
или с еще большим количеством кода:
import numpy
def primesgen(n):
""" Input n>=30, Generates all primes < n """
size = n/30 + 1
sieve01 = numpy.ones(size, dtype=numpy.bool)
sieve07 = numpy.ones(size, dtype=numpy.bool)
sieve11 = numpy.ones(size, dtype=numpy.bool)
sieve13 = numpy.ones(size, dtype=numpy.bool)
sieve17 = numpy.ones(size, dtype=numpy.bool)
sieve19 = numpy.ones(size, dtype=numpy.bool)
sieve23 = numpy.ones(size, dtype=numpy.bool)
sieve29 = numpy.ones(size, dtype=numpy.bool)
sieve01[0] = False
yield 2; yield 3; yield 5;
for i in xrange(int(n**0.5)/30+1):
if sieve01[i]:
k=30*i+1; yield k;
sieve01[ (k*k)/30::k] = False
sieve07[(k*k+ 6*k)/30::k] = False
sieve11[(k*k+10*k)/30::k] = False
sieve13[(k*k+12*k)/30::k] = False
sieve17[(k*k+16*k)/30::k] = False
sieve19[(k*k+18*k)/30::k] = False
sieve23[(k*k+22*k)/30::k] = False
sieve29[(k*k+28*k)/30::k] = False
if sieve07[i]:
k=30*i+7; yield k;
sieve01[(k*k+ 6*k)/30::k] = False
sieve07[(k*k+24*k)/30::k] = False
sieve11[(k*k+16*k)/30::k] = False
sieve13[(k*k+12*k)/30::k] = False
sieve17[(k*k+ 4*k)/30::k] = False
sieve19[ (k*k)/30::k] = False
sieve23[(k*k+22*k)/30::k] = False
sieve29[(k*k+10*k)/30::k] = False
if sieve11[i]:
k=30*i+11; yield k;
sieve01[ (k*k)/30::k] = False
sieve07[(k*k+ 6*k)/30::k] = False
sieve11[(k*k+20*k)/30::k] = False
sieve13[(k*k+12*k)/30::k] = False
sieve17[(k*k+26*k)/30::k] = False
sieve19[(k*k+18*k)/30::k] = False
sieve23[(k*k+ 2*k)/30::k] = False
sieve29[(k*k+ 8*k)/30::k] = False
if sieve13[i]:
k=30*i+13; yield k;
sieve01[(k*k+24*k)/30::k] = False
sieve07[(k*k+ 6*k)/30::k] = False
sieve11[(k*k+ 4*k)/30::k] = False
sieve13[(k*k+18*k)/30::k] = False
sieve17[(k*k+16*k)/30::k] = False
sieve19[ (k*k)/30::k] = False
sieve23[(k*k+28*k)/30::k] = False
sieve29[(k*k+10*k)/30::k] = False
if sieve17[i]:
k=30*i+17; yield k;
sieve01[(k*k+ 6*k)/30::k] = False
sieve07[(k*k+24*k)/30::k] = False
sieve11[(k*k+26*k)/30::k] = False
sieve13[(k*k+12*k)/30::k] = False
sieve17[(k*k+14*k)/30::k] = False
sieve19[ (k*k)/30::k] = False
sieve23[(k*k+ 2*k)/30::k] = False
sieve29[(k*k+20*k)/30::k] = False
if sieve19[i]:
k=30*i+19; yield k;
sieve01[ (k*k)/30::k] = False
sieve07[(k*k+24*k)/30::k] = False
sieve11[(k*k+10*k)/30::k] = False
sieve13[(k*k+18*k)/30::k] = False
sieve17[(k*k+ 4*k)/30::k] = False
sieve19[(k*k+12*k)/30::k] = False
sieve23[(k*k+28*k)/30::k] = False
sieve29[(k*k+22*k)/30::k] = False
if sieve23[i]:
k=30*i+23; yield k;
sieve01[(k*k+24*k)/30::k] = False
sieve07[(k*k+ 6*k)/30::k] = False
sieve11[(k*k+14*k)/30::k] = False
sieve13[(k*k+18*k)/30::k] = False
sieve17[(k*k+26*k)/30::k] = False
sieve19[ (k*k)/30::k] = False
sieve23[(k*k+ 8*k)/30::k] = False
sieve29[(k*k+20*k)/30::k] = False
if sieve29[i]:
k=30*i+29; yield k;
sieve01[ (k*k)/30::k] = False
sieve07[(k*k+24*k)/30::k] = False
sieve11[(k*k+20*k)/30::k] = False
sieve13[(k*k+18*k)/30::k] = False
sieve17[(k*k+14*k)/30::k] = False
sieve19[(k*k+12*k)/30::k] = False
sieve23[(k*k+ 8*k)/30::k] = False
sieve29[(k*k+ 2*k)/30::k] = False
for i in xrange(i+1,n/30):
if sieve01[i]: yield 30*i+1
if sieve07[i]: yield 30*i+7
if sieve11[i]: yield 30*i+11
if sieve13[i]: yield 30*i+13
if sieve17[i]: yield 30*i+17
if sieve19[i]: yield 30*i+19
if sieve23[i]: yield 30*i+23
if sieve29[i]: yield 30*i+29
i += 1
mod30 = n%30
if mod30 > 1:
if sieve01[i]: yield 30*i+1
if mod30 > 7:
if sieve07[i]: yield 30*i+7
if mod30 > 11:
if sieve11[i]: yield 30*i+11
if mod30 > 13:
if sieve13[i]: yield 30*i+13
if mod30 > 17:
if sieve17[i]: yield 30*i+17
if mod30 > 19:
if sieve19[i]: yield 30*i+19
if mod30 > 23:
if sieve23[i]: yield 30*i+23
Ps: Если у вас есть какие-то предложения о том, как ускорить этот код, что-либо от порядка операций операций до предварительного вычисления чего-либо, прокомментируйте.
Ответ 5
Здесь версия, которую я написал некоторое время назад; это может быть интересно сравнить с вашим для скорости. Однако он ничего не делает о космических проблемах (на самом деле они, вероятно, хуже, чем с вашей версией).
from math import sqrt
def basicSieve(n):
"""Given a positive integer n, generate the primes < n."""
s = [1]*n
for p in xrange(2, 1+int(sqrt(n-1))):
if s[p]:
a = p*p
s[a::p] = [0] * -((a-n)//p)
for p in xrange(2, n):
if s[p]:
yield p
У меня есть более быстрые версии, используя колесо, но они намного сложнее. Исходный источник здесь.
Хорошо, здесь версия с использованием колеса. wheelSieve - основная точка входа.
from math import sqrt
from bisect import bisect_left
def basicSieve(n):
"""Given a positive integer n, generate the primes < n."""
s = [1]*n
for p in xrange(2, 1+int(sqrt(n-1))):
if s[p]:
a = p*p
s[a::p] = [0] * -((a-n)//p)
for p in xrange(2, n):
if s[p]:
yield p
class Wheel(object):
"""Class representing a wheel.
Attributes:
primelimit -> wheel covers primes < primelimit.
For example, given a primelimit of 6
the wheel primes are 2, 3, and 5.
primes -> list of primes less than primelimit
size -> product of the primes in primes; the modulus of the wheel
units -> list of units modulo size
phi -> number of units
"""
def __init__(self, primelimit):
self.primelimit = primelimit
self.primes = list(basicSieve(primelimit))
# compute the size of the wheel
size = 1
for p in self.primes:
size *= p
self.size = size
# find the units by sieving
units = [1] * self.size
for p in self.primes:
units[::p] = [0]*(self.size//p)
self.units = [i for i in xrange(self.size) if units[i]]
# number of units
self.phi = len(self.units)
def to_index(self, n):
"""Compute alpha(n), where alpha is an order preserving map
from the set of units modulo self.size to the nonnegative integers.
If n is not a unit, the index of the first unit greater than n
is given."""
return bisect_left(self.units, n % self.size) + n // self.size * self.phi
def from_index(self, i):
"""Inverse of to_index."""
return self.units[i % self.phi] + i // self.phi * self.size
def wheelSieveInner(n, wheel):
"""Given a positive integer n and a wheel, find the wheel indices of
all primes that are less than n, and that are units modulo the
wheel modulus.
"""
# renaming to avoid the overhead of attribute lookups
U = wheel.units
wS = wheel.size
# inverse of unit map
UI = dict((u, i) for i, u in enumerate(U))
nU = len(wheel.units)
sqroot = 1+int(sqrt(n-1)) # ceiling of square root of n
# corresponding index (index of next unit up)
sqrti = bisect_left(U, sqroot % wS) + sqroot//wS*nU
upper = bisect_left(U, n % wS) + n//wS*nU
ind2 = bisect_left(U, 2 % wS) + 2//wS*nU
s = [1]*upper
for i in xrange(ind2, sqrti):
if s[i]:
q = i//nU
u = U[i%nU]
p = q*wS+u
u2 = u*u
aq, au = (p+u)*q+u2//wS, u2%wS
wp = p * nU
for v in U:
# eliminate entries corresponding to integers
# congruent to p*v modulo p*wS
uvr = u*v%wS
m = aq + (au > uvr)
bot = (m + (q*v + u*v//wS - m) % p) * nU + UI[uvr]
s[bot::wp] = [0]*-((bot-upper)//wp)
return s
def wheelSieve(n, wheel=Wheel(10)):
"""Given a positive integer n, generate the list of primes <= n."""
n += 1
wS = wheel.size
U = wheel.units
nU = len(wheel.units)
s = wheelSieveInner(n, wheel)
return (wheel.primes[:bisect_left(wheel.primes, n)] +
[p//nU*wS + U[p%nU] for p in xrange(bisect_left(U, 2 % wS)
+ 2//wS*nU, len(s)) if s[p]])
Что касается колеса: хорошо, вы знаете, что (кроме 2) все простые числа странны, поэтому большинство сит пропускают все четные числа. Аналогично, вы можете пойти немного дальше и заметить, что все простые числа (кроме 2 и 3) соответствуют 1 или 5 по модулю 6 (== 2 * 3), поэтому вы можете уйти с сохранением записей для этих чисел в вашем сите, Следующий шаг - отметить, что все простые числа (кроме 2, 3 и 5) соответствуют одному из 1, 7, 11, 13, 17, 19, 23, 29 (по модулю 30) (здесь 30 == 2 * 3 * 5) и т.д.
Ответ 6
Одним из улучшений скорости, которые вы можете сделать с помощью bitstring, является использование метода set, когда вы исключаете кратные текущего номера.
Итак, жизненный раздел становится
for i in range(3, limit, 2):
if flags[i]:
yield i
if i <= sub_limit:
flags.set(1, range(i*3, limit, i*2))
На моей машине это выполняется примерно в 3 раза быстрее, чем оригинал.
Моя другая мысль заключалась в том, чтобы использовать bitstring для представления только нечетных чисел. Затем вы можете найти несохраненные биты, используя flags.findall([0]), который возвращает генератор. Не уверен, что это будет намного быстрее (найти битовые шаблоны не так просто сделать эффективно).
[Полное раскрытие: я написал модуль битовой строки, так что у меня есть гордость здесь!]
В качестве сравнения я также взял кишки из метода bitstring, чтобы он делал это точно так же, но без проверки диапазона и т.д. Я думаю, что это дает разумный нижний предел для чистого Python, который работает на миллиард элементы (без изменения алгоритма, который я считаю обманом!)
def prime_pure(limit=1000000):
yield 2
flags = bytearray((limit + 7) // 8)
sub_limit = int(limit**0.5)
for i in xrange(3, limit, 2):
byte, bit = divmod(i, 8)
if not flags[byte] & (128 >> bit):
yield i
if i <= sub_limit:
for j in xrange(i*3, limit, i*2):
byte, bit = divmod(j, 8)
flags[byte] |= (128 >> bit)
На моей машине это составляет около 0.62 секунды для миллиона элементов, что означает, что это примерно четверть скорости ответа битаррей. Я считаю это вполне разумным для чистого Python.
Ответ 7
Целочисленные типы Python могут быть произвольного размера, поэтому вам не нужна умная библиотека битовой строки для представления набора бит, всего лишь одно число.
Вот код. Он легко справляется с лимитом в 1 000 000 человек и может даже обрабатывать 10 000 000, не жалуясь слишком много:
def multiples_of(n, step, limit):
bits = 1 << n
old_bits = bits
max = 1 << limit
while old_bits < max:
old_bits = bits
bits += bits << step
step *= 2
return old_bits
def prime_numbers(limit=10000000):
'''Prime number generator. Yields the series
2, 3, 5, 7, 11, 13, 17, 19, 23, 29 ...
using Sieve of Eratosthenes.
'''
yield 2
sub_limit = int(limit**0.5)
flags = ((1 << (limit - 2)) - 1) << 2
# Step through all the odd numbers
for i in xrange(3, limit, 2):
if not (flags & (1 << i)):
continue
yield i
# Exclude further multiples of the current prime number
if i <= sub_limit:
flags &= ~multiples_of(i * 3, i << 1, limit)
Функция multiples_of позволяет избежать затрат на установку каждого отдельного набора по отдельности. Он устанавливает начальный бит, сдвигает его на начальный шаг (i << 1) и добавляет результат сам по себе. Затем он удваивает шаг, так что следующий сдвиг копирует оба бита и так далее, пока не достигнет предела. Это устанавливает все кратные числа до предела в O (log (limit)).
Ответ 8
Один из вариантов, который вы, возможно, захотите посмотреть, - это просто компилировать модуль C/С++, чтобы вы имели прямой доступ к функциям бит-twiddling. AFAIK вы могли бы написать что-то такого характера и только возвращать результаты по завершении вычислений, выполненных на C/С++. Теперь, когда я набираю это, вы также можете посмотреть на numpy, который хранит массивы как родные C для скорости. Я не знаю, будет ли это быстрее, чем модуль bitstring!
Ответ 9
Другой действительно глупый вариант, но это может помочь, если вам действительно нужен большой список простых чисел, доступных очень быстро. Скажем, если вы нуждаетесь в них в качестве инструмента для ответа на проблемы проекта Эйлера (если проблема в том, что просто оптимизировать ваш код как игру ума, это не имеет значения).
Используйте любое медленное решение для создания списка и сохраните его в исходном файле python, говорит primes.py, который просто содержит:
primes = [ a list of a million primes numbers here ]
Теперь, чтобы использовать их, вам просто нужно сделать import primes, и вы получите их с минимальным размером памяти на скорости дискового ввода-вывода. Сложность - количество простых чисел: -)
Даже если вы использовали плохо оптимизированное решение для создания этого списка, это будет сделано только один раз, и это не имеет большого значения.
Возможно, вы можете сделать это еще быстрее, используя pickle/unpickle, но вы получите идею...
Ответ 10
Вы можете использовать сегментированное сито Эратосфена. Память, используемая для каждого сегмента, повторно используется для следующего сегмента.
Ответ 11
Вот некоторый код Python3, который использует меньше памяти, чем решения bitarray/bitstring, опубликованные на сегодняшний день, и около 1/8 памяти Роберта Уильяма Хэнкса primesgen(), при этом работает быстрее, чем primesgen() (немного быстрее на 1,000,000, используя 37 Кбайт памяти, в 3 раза быстрее, чем primesgen() на 1,000,000,000, используя до 34 МБ). Увеличение размера блока (переменная часть кода) использует больше памяти, но ускоряет работу программы, до предела - я выбрал значение, чтобы его вклад в память составлял около 10% от решетки n//30 байт. Это не так эффективно для памяти, как бесконечный генератор Will Ness (fooobar.com/info/9543/...), потому что он записывает (и возвращает в конце в сжатой форме) все рассчитанные простые числа.
Это должно работать корректно, если вычисление квадратного корня получается точным (около 2 ** 51, если Python использует 64-битные удвоения). Однако вы не должны использовать эту программу для поиска больших чисел!
(Я не пересчитывал смещения, просто взял их из кода Роберта Уильяма Хэнкса. Спасибо, Роберт!)
import numpy as np
def prime_30_gen(n):
""" Input n, int-like. Yields all primes < n as Python ints"""
wheel = np.array([2,3,5])
yield from wheel[wheel < n].tolist()
powers = 1 << np.arange(8, dtype='u1')
odds = np.array([1, 7, 11, 13, 17, 19, 23, 29], dtype='i8')
offsets=np.array([[0,6,10,12,16,18,22,28],[6,24,16,12,4,0,22,10],
[0,6,20,12,26,18,2,8], [24,6,4,18,16,0,28,10],
[6,24,26,12,14,0,2,20], [0,24,10,18,4,12,28,22],
[24,6,14,18,26,0,8,20], [0,24,20,18,14,12,8,2]],
dtype='i8')
offsets = {d:f for d,f in zip(odds, offsets)}
sieve = 255 * np.ones((n + 28) // 30, dtype='u1')
if n % 30 > 1: sieve[-1] >>= 8 - odds.searchsorted(n % 30)
sieve[0] &= ~1
for i in range((int((n - 1)**0.5) + 29) // 30):
for odd in odds[(sieve[i] & powers).nonzero()[0]]:
k = i * 30 + odd
yield int(k)
for clear_bit, off in zip(~powers, offsets[odd]):
sieve[(k * (k + off)) // 30 :: k] &= clear_bit
chunk = min(1 + (n >> 13), 1<<15)
for i in range(i + 1, len(sieve), chunk):
a = (sieve[i : i + chunk, np.newaxis] & powers)
a = np.flatnonzero(a.astype('bool')) + i * 8
a = (a // 8 * 30 + odds[a & 7]).tolist()
yield from a
return sieve
Примечание: если вам нужна реальная скорость и требуется 2 ГБ оперативной памяти для списка простых чисел до 10 ** 9, вы должны использовать pyprimesieve (на https://pypi.python.org/, используя primesieve http://primesieve.org/).
