素数抽出アルゴリズムの解読

素数抽出アルゴリズムの解読:

前回記事のコメントにて、改善版のアルゴリズムを頂いたのでそれの解読の記録。

主に上記記事の自作コードとの比較です。

私はこの記事の中で何度自分のコードを無駄呼ばわりすることになるのでしょうか。

注意事項

• 記載いただいたコードは本来引用記法にするべきと思いますが、インデントが崩れるためコード記法にしています。
• 記載いただいたコードのうち、処理時間計測の部分は省略させて頂きました。
• ※未解読の部分は【考え中】としています。自力で気付きたいので、解説と思しきコメントを頂いてもすぐには返信コメントができません。ご了承下さい。

@shiracamus さんのPython

元コメント

import datetime 
 
MAX = 100000 
 
# print(2) 
# print(3) 
count = 2 
 
for i in range(5, MAX + 1, 2): 
    if all(i % j != 0 for j in range(3, int(i**0.5) + 1, 2)): 
        # print(i) 
        count += 1 
 
print("total count : " + str(count)) 

自作コードとの相違点

1. ループのiの初期数が5
2. 
   iのインクリメントが+2
3. 割り算の母数の最大値が int(i**0.5)+1
   

ループのiの初期数が5

自作コードでは4でした。4は当然素数ではないので不要です。アホでした。

iのインクリメントが+2

偶数を除外できます。

偶数の除外自体は自作コードでも行っていたのですが、インクリメントが+1のfor文の中でif (i%2 == 0):continueという分岐を設けるというものでした。

頂いたコードのようにインクリメントを+2にすることで、単純にループの回数を半分にできます。大幅短縮。

割り算の母数の最大値が int(i**0.5)+1

int(i**0.5)はiの平方根の少数切捨て。

自作コードではround(i/2)としていましたが、この場合、商が2になるまで計算していることになります。これでは無駄な計算が生じます。

以下詳細解説

例としてi=97の場合、

• 97/2 = 48 余り1
• 97/48 = 2 余り1

の2つの計算は、素数を求めたい（=余りが0になるかどうかを知りたい）という目的からすると重複しているので、二番目の計算は無駄です。

ではどこまでなら無駄にならないのかというと、境目はここ。

• 97/8 = 11 余り 9
• 97/9 = 9 余り 16
• 97/10 = 9 余り 7
• 
  97/11 = 8 余り 9　★97/8とやっていることは同じ（=無駄）

割り算の母数が97の平方根を超えてしまえば、あとは既に行っている計算の逆バージョンとなり、余りが欲しいという目的からすると無駄。

97の平方根は9.8488...ですので、整数で言うのであれば10まで。

これを式にすると、int(i**0.5)+1となるわけです。

@hogefuga さんのJavaScript①

元コメント

const MAX = 100000; 
 
let count = 2; 
 
for (let i = 5, d = -1; i <= MAX; i += d + 3, d = -d) { 
 
   let isPrime = true; 
 
   const max = i ** 0.5 | 0; 
   for (let j = 5, e = -1; j <= max; j += e + 3, e = -e) { 
 
      if (!(i % j)) { 
         isPrime = false; 
         break; 
      } 
 
   } 
 
   if (isPrime) { 
      count++; 
   } 
 
} 
 
console.log(`total count: ${ count }`); 

自作コードとの相違点

沢山ありすぎてもう何がなんだか

1. 割り算の母数の最大値が i ** 0.5 | 0【考え中】
2. ループfor (let i = 5, d = -1; i <= MAX; i += d + 3, d = -d)
   

割り算の母数の最大値が i ** 0.5 | 0 【考え中】

意味としてはiの平方根の少数切り上げになります。@shiracamus さんのPythonと同じですね。Math.floor(sqrt(i))です。

【考え中】なぜbit演算でそのような値が出るのかまだいまいちわかっていません。

ループfor (let i = 5, d = -1; i <= MAX; i += d + 3, d = -d)

動作としては、

• 
  @shiracamus さんのPythonと同様、5からスタート
• インクリメントは、+2と+4を交互に行う

というものとなります。

それで得られる数列はこちら。

5 7 11 13 17 19 23 25 29 31 35 37 41 43 47 49 53 ... 

偶数だけではなく、3の倍数まで排除した数列になっています。

かなり素数率高い！！すげええええ！！！！

私なりの詳細解説

変数dは延々と-1と+1をスイッチしていきますので、iのインクリメントは

+=(-1+3), +=(1+3), +=(-1+3), +=(1+3),...

すなわち

+=2, +=4, +=2, +=4,...

となります。

なぜそれで3の倍数が排除できるのか？

一旦偶数はアリとして、3の倍数だけを排除したい場合を考えてみましょう。

偶数だけを排除したい場合はインクリメントを+2にすれば済みましたが、同じようにインクリメントを+3にしてしまうと素数の取りこぼしが発生します。

5 8 11 14 17 20 23  

確かに3の倍数は排除できていますが、これだけでも既に素数である7, 13, 19が取りこぼされています。

偶数と違い、「3の倍数以外の数」は2種類あることを考慮しなければなりません。

2種類とは、3で割った時の「余りが1の数」と「余りが2の数」です。

上記の配列では、5に始まった「余りが2の数」しか拾えていないわけです。

では2種類を同時に拾うにはどうしたらよいのかというと、

• 「余りが1の数」には1を足す（「余りが2の数」になる）
• 「余りが2の数」には2を足す（「余りが1の数」になる）

を繰り返せばよいのです。

上記の配列を修正すると、スタートは5=「余りが2の数」ですので、
+=2, +=1, +=2, +=1, +=2, +=1,...

としていくと、

5 7 8 10 11 13 14 16 17 19 20 22 23 25 26 28 29 31 32 34 35 

という数列が得られます。これなら素数を取りこぼすことなく3の倍数だけを排除できています。

さて。

お分かり頂けただろうか。

この数列、等間隔に偶数が出現していることに。

5 7 (8 10) 11 13 (14 16) 17 19 (20 22) 23 25 (26 28) 29 31 (32 34) 35 

+1と+2を繰り返せば、「奇数+1で偶数」「偶数+2で偶数」が並ぶのは考えてみれば当たり前ですね。頂いたコードを見るまで考え付きもしませんでしたけど

等間隔なのであればforループで何とかなるはずです。

上記配列の偶数部分をすっ飛ばすのであれば、インクリメントは、
+=2, (+=1+2+1), +=2, (+=1+2+1), +=2, (+=1+2+1),...

となります。計算すると
+=2, +=4, +=2, +=4, +=2, +=4,...

です。

頂いたコードのインクリメントそのままです。
ファンタスティック！！

私だったらここまでたどり着いても、+2のfor文の中で場合分けして一回おきに更に+2、とかしてしまいそうですが、そんな無駄をせずスマートに書いたのが

for (let i = 5, d = -1; i <= MAX; i += d + 3, d = -d)

なわけです。Huuuuuuuu!!!

@hogefuga さんのC++ 【考え中】

元コメント

#include<iostream> 
 
constexpr int count_prime(const int max) { 
   int count = 2; 
   for (int i = 5, d = -1, r = 3, s = 9; i <= max; i += d + 3, d = -d) { 
      bool is_prime = true; 
      if (i >= s) { 
         r++; 
         s = r * r; 
      } 
      for (int j = 5, e = -1; j <= r; j += e + 3, e = -e) { 
         if (i % j == 0) { 
            is_prime = false; 
            break; 
         } 
      } 
      if (is_prime) { 
         count++; 
      } 
   } 
   return count; 
} 
 
int main() { 
   constexpr int count = count_prime(100000); 
   static_assert(count == 9592); 
   std::cout << "count: " << count << std::endl; 
} 

@hogefuga さんのJavaScript② 【考え中】

元コメント

const MAX = 100000; 
 
let count = 2; 
for (let i = 5, d = 1, r = 3, s = 9, n = 7; i <= MAX; i += (d = -d) + 3) 
   for (let j = 5, e = ((i >= s && (++r, s += n, n += 2)), 1); j <= r || !++count; j += (e = -e) + 3) 
      if (!(i % j)) break; 
 
console.log(`total count: ${count}`);