Некоторое время назад я использовал (молниеносно) простейшее в python, которое нашел здесь: Самый быстрый способ перечислить все простые числа ниже N
Чтобы быть точным, эта реализация:
def primes2(n):
""" Input n>=6, Returns a list of primes, 2 <= p < n """
n, correction = n-n%6+6, 2-(n%6>1)
sieve = [True] * (n/3)
for i in xrange(1,int(n**0.5)/3+1):
if sieve[i]:
k=3*i+1|1
sieve[ k*k/3 ::2*k] = [False] * ((n/6-k*k/6-1)/k+1)
sieve[k*(k-2*(i&1)+4)/3::2*k] = [False] * ((n/6-k*(k-2*(i&1)+4)/6-1)/k+1)
return [2,3] + [3*i+1|1 for i in xrange(1,n/3-correction) if sieve[i]]
Теперь я могу немного понять идею оптимизации, автоматически пропуская кратные 2, 3 и т. Д., Но когда дело доходит до переноса этого алгоритма на C ++, я застреваю (у меня хорошее понимание Python и разумный / плохое понимание C ++, но достаточно хорошее для рок-н-ролла).
То, что я сейчас прокрутил сам, это (isqrt() - это просто целочисленная функция квадратного корня):
template <class T>
void primesbelow(T N, std::vector<T> &primes) {
T sievemax = (N-3 + (1-(N % 2))) / 2;
T i;
T sievemaxroot = isqrt(sievemax) + 1;
boost::dynamic_bitset<> sieve(sievemax);
sieve.set();
primes.push_back(2);
for (i = 0; i <= sievemaxroot; i++) {
if (sieve[i]) {
primes.push_back(2*i+3);
for (T j = 3*i+3; j <= sievemax; j += 2*i+3) sieve[j] = 0; // filter multiples
}
}
for (; i <= sievemax; i++) {
if (sieve[i]) primes.push_back(2*i+3);
}
}
Эта реализация достойна и автоматически пропускает кратные 2, но если бы я мог портировать реализацию Python, я думаю, что она могла бы быть намного быстрее (50% -30% или около того).
Для сравнения результатов (в надежде, что этот вопрос будет успешно дан), текущее время выполнения с N=100000000, g++ -O3 на Q6600 Ubuntu 10.10 составляет 1230 мс.
Теперь я хотел бы получить некоторую помощь в понимании того, что делает приведенная выше реализация Python, или в том, что вы перенесли бы его для меня (хотя и не так полезно).
EDIT
Некоторая дополнительная информация о том, что мне трудно.
У меня проблемы с такими методами, как корректирующая переменная и вообще как она получается. Ссылка на сайт, объясняющий различные оптимизации Eratosthenes (кроме простых сайтов, которые говорят: «ну, вы просто пропустите кратные 2, 3 и 5», а затем получите файл с 1000 строчками C), было бы здорово.
Не думаю, что у меня возникли бы проблемы со 100% прямым и буквальным портом, но, в конце концов, это для обучения, которое было бы совершенно бесполезным.
EDIT
Посмотрев на код в оригинальной numpy версии, он на самом деле довольно прост для реализации и с некоторыми соображениями, не слишком сложными для понимания. Это версия C ++, которую я придумал. Я публикую его здесь в полной версии, чтобы помочь другим читателям в случае, если им понадобится довольно эффективное простое изображение, не содержащее двух миллионов строк кода. Это простое число делает все простые числа меньше 100000000 примерно за 415 мс на той же машине, что и выше. Это в 3 раза быстрее, чем я ожидал!
#include <vector>
#include <boost/dynamic_bitset.hpp>
// http://vault.embedded.com/98/9802fe2.htm - integer square root
unsigned short isqrt(unsigned long a) {
unsigned long rem = 0;
unsigned long root = 0;
for (short i = 0; i < 16; i++) {
root <<= 1;
rem = ((rem << 2) + (a >> 30));
a <<= 2;
root++;
if (root <= rem) {
rem -= root;
root++;
} else root--;
}
return static_cast<unsigned short> (root >> 1);
}
// /2029352/samyi-bystryi-sposob-perechislit-vse-prostye-chisla-nizhe-n#2029364
// https://stackoverflow.com/questions/5293238/porting-optimized-sieve-of-eratosthenes-from-python-to-c/5293492
template <class T>
void primesbelow(T N, std::vector<T> &primes) {
T i, j, k, l, sievemax, sievemaxroot;
sievemax = N/3;
if ((N % 6) == 2) sievemax++;
sievemaxroot = isqrt(N)/3;
boost::dynamic_bitset<> sieve(sievemax);
sieve.set();
primes.push_back(2);
primes.push_back(3);
for (i = 1; i <= sievemaxroot; i++) {
if (sieve[i]) {
k = (3*i + 1) | 1;
l = (4*k-2*k*(i&1)) / 3;
for (j = k*k/3; j < sievemax; j += 2*k) {
sieve[j] = 0;
sieve[j+l] = 0;
}
primes.push_back(k);
}
}
for (i = sievemaxroot + 1; i < sievemax; i++) {
if (sieve[i]) primes.push_back((3*i+1)|1);
}
}