は、PythonからC++

Question

にエラトステネスのふるいを最適化された移植私はここで見つけるPythonで（高速連射）primesieve使用前にいくつかの時間：今

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]] 

：Fastest way to list all primes below N

は、正確には、この実装を私は自動的に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％程度）することができると思う。

Q6600 Ubuntu 10.10でのN=100000000,g++ -O3の現在の実行時間は1230msです（この質問にはうまく答えられるはずです）。

私は、上記のPython実装が何をしているのか、それとも私のために移植するのかを理解するのに役立つでしょう。

EDIT

私が難しいものについてのいくつかの追加情報。

訂正変数のように使われている手法や一般的にどのようになっているのかよくわかりません。さまざまなEratosthenesの最適化を説明しているサイトへのリンク（「2,3,5の倍数をスキップして1000行のCファイルでスラムする」という単純なサイトを除いて）はすばらしいでしょう。

私は100％直接ポートとリテラルポートに問題はないと思っていますが、結局これはまったく役に立たない学習のためです。

EDIT

オリジナルnumpyのバージョンのコードを見ていたら、実際にはかなり実装が容易と理解することはそれほど難しくなく、いくつかの考え方です。これは私が思いついたC++のバージョンです。私は200万行のコードではない非常に効率的なprimesieveが必要な場合に備えて、読者の皆様の助けになるよう完全版でここに掲載しています。この素数は、上記と同じマシン上で約415ミリ秒で100000000未満のすべての素数を実行します。これは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); 
} 

// https://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188 
// 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); 
    } 
} 

Answer 1

私は可能な限り説明しようとします。 sieveアレイには通常とは異なるインデックス方式があります。 1または5 mod 6に一致する各番号のビットを格納します。したがって、数字6*k + 1は位置2*kに格納され、k*6 + 5は位置2*k + 1に格納されます。 3*i+1|1操作はその逆である：それは、フォーム2*nの数字を取り、6*n + 1に変換して2*n + 1を取り、（+1|1事が5に1と3に0を変換）6*n + 5に変換します。メインループは、5（iが1の場合）で始まる、そのプロパティを持つすべての数値を通じてkを繰り返します。 iはkのsieveに対応するインデックスです。 sieveに最初のスライスを更新すると、k*k/3 + 2*m*k（mの自然数）という形式のインデックスを持つふるい内のすべてのビットがクリアされます。それらのインデックスの対応する番号はk^2で始まり、各ステップで6*kだけ増加します。第2のスライス更新は、インデックスk*(k-2*(i&1)+4)/3（kの場合は数字k * (k+4)、それ以外の場合は1 mod 6およびk * (k+2)に一致する）から開始し、同様に各ステップで6*kの数を増加させます。

は、ここでの説明では、別の試みです：candidatesは少なくとも5であり、どちらか1または5 MOD 6と合同であるすべての数字の集合とします。その集合の2つの要素を掛け合わせると、その集合の中の別の要素が得られます。いくつかのkがcandidatesのsucc(k)をkより大きいcandidatesの次の要素（数値順）とします。その場合には、篩の内側ループは、（sieveの通常インデックス作成を使用して）基本的に：

for k in candidates: 
    for l in candidates where l >= k: 
    sieve[k * l] = False 

理由：要素が sieveに保存されている制限のため、それは同じである

for k in candidates: 
    for (l = k; ; l += 6) sieve[k * l] = False 
    for (l = succ(k); ; l += 6) sieve[k * l] = False 

（現在kは、それが今kとして使用される以前lやとして使用されたいずれかの場合）、ある時点でふるいから（k自体以外）candidatesにkのすべての倍数を除去する

。

Answer 2

脇に、素数を「近似」することができます。おおよその素数Pを呼び出すここではいくつかの式である：

P = 2 * K + 1 //ない割り切れ

P = 6 2によっては、* K + {1,5} 2 //割り切れません、 3

P見つかった数字のセットの5つの

特性、3,2による= 30 * K + {1、7、11、13、17、19、23、29} // divisbleありませんこれらの式によって、Pは素数ではないかもしれないが、すべての素数は集合PあなたがプライムのためにセットPの数字だけをテストすれば、あなたは何も見逃しません。

あなたはこれらの式を定式ことができます：それはあなたのためのより便利であれば

P = X * K + {-i、-j、-k、K、J、I}

。

Here 2で割り切れないPについての式を用いて、この手法を使用し、いくつかのコード、3、5、7

このリンクは、この技術が実際に活用することができる範囲を表すことができるされています。

Answer 3

Howardは、ハワードの応答にピギーバックをかけて、素数性のために2,3,5で割り切れないすべての自然数をテストする必要はありません。配列内の各番号（1を除いて、自己消去する）と、配列内のそれ以降のすべての番号を単純に乗算する必要があります。これらの重複プロダクトは、決定論的乗法プロセスを拡張する点まで、配列内のすべての非素数をあなたに与えます。したがって、配列の最初の非プライムは7の平方または49、2番目、7番目の11、または77などです。ここでの完全な説明：http://www.primesdemystified.com