LoginSignup
0
0

More than 5 years have passed since last update.

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

Last updated at Posted at 2018-12-05

前回記事のコメントにて、改善版のアルゴリズムを頂いたのでそれの解読の記録。
主に上記記事の自作コードとの比較です。
私はこの記事の中で何度自分のコードを無駄呼ばわりすることになるのでしょうか。

注意事項

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

2018.12.16 解読完了しました。大変勉強になりました。
コメントくださった @shiracamus さん、 @hogefuga さん、ありがとうございました!!

@shiracamus さんのPython

元コメント

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演算でそのような値が出るのかまだいまいちわかっていませんでしたが、こちらの記事JavaScript イディオム集 に、v | 0 について解説がありました!!
スッキリ!

ループ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;
}

自作コードとの相違点

もう悲しくなってきたので省略していいですか。

解読結果(カウントロジック部分)

まずfor文の定義部分は、i = 5, d = -1; i <= max; i += d+3, d=-dだけ見るとJavaScript①と同様。
それで得られるiの数列は以下です。5以上で偶数と3の倍数を除いた整数ですね。

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

今回は変数がr =3, s = 9と2つ追加されているので、上記数列にこの2変数がどう作用しているかが問題となります。
なんかインクリメントされたり二乗にされたりしていますが、実際に使われるところは二重for文の定義でのj <= r;のみです。

rとsの動きを見ますと、iがs = r*rを超えた時のみrをインクリメントしています。
なので愚直に書き出すと、

i < 3*3 :
    j <= r = 3
i < 4*4 :
    j <= r = 4
i < 5*5 :
    j <= r = 5

という動きに。あれ?これってつまり…

jは素数かを確認するための割り算の母数。上述JavaScript①ではj <= i**0.5 | 0となっていました、iの平方根を超える数で割っても計算の無駄、という話ですね。
rとsの役割もやっていることはこの「平方根(の繰り上げ整数)を超える数では割らない」という制限になります。

こういう書き方もできるんですね。

assert

さりげなく追加されていますが、このassert、MAXが大きくなればなるほど効いてきそうです。
今回(MAX10万)程度ならまだいいんですが、これが大きくなってくると結果として出てきたcountも当然大きくなり、合ってるのかどうかを目視確認するのも面倒になってきますので。
このassertがあるおかげで、確認は(素数の数が判明済のMAXであれば)正解をコピペしてassertしてもらえば、イケてるかどうかが一目瞭然になるわけです。
悲しい話ですがコピペは目視よりも確実です。

なおsは上記のiの判定一回にしか使われておらず、必ずrの二乗になりますので、省略してこう書くことも可能と思われます。手元にC++の環境がないので試せてないのですが。

【未確認】

for (int i = 5, d = -1, r = 3; i <= max; i += d + 3, d = -d) {
      bool is_prime = true;
      if (i >= (r*r)) {
         r++;
     }

@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}`);

解読

ご本人もコメントで難解と仰ってました通り、これは手こずりました。

for文一つ目の定義は上記のC++に似ています。相違点は、

  • dの初期値が-1ではなく1
  • 新キャラn = 7

の2点。ここまではまだなんとかなりそうですが、二重目のfor文の定義がすごい。

for ( let j = 5, e = ((i >= s && (++r, s += n, n += 2)), 1); j <= r || !++count; j += (e = -e) + 3)

ここでほぼ全てを解決しているのか、for文の中でやっているのは割り算だけです。
ちょっと何言ってるかよくわからないので分解していきましょう。

    let j = 5, 

ここはいつも通り。(最初に書いた自分のコードなど忘れました)
問題はこの後。

eの初期化と見せかけて色々な処理を詰め込んだ部分

    e = ((i >= s && (++r, s += n, n += 2)), 1)

あまりにも分からないので動かしながら見ます。
結論から言うとrはJavaScript①で言うところのi ** 0.5 | 0に相当します。話は同じでした。

どういうわけなのかコードを簡略化して調べてみましょう。

const MAX = 50;
for (let i = 5, d = 1, r = 3, s = 9, n = 7; i <= MAX; i += (d = -d) + 3) {
  let e = ((i >= s && (++r, s += n, n += 2)), 1);
  console.log(`i:${i} r:${r} s:${s} n:${n} e:${e}`);
}

結果は

"i:5 r:3 s:9 n:7 e:1"
"i:7 r:3 s:9 n:7 e:1"
"i:11 r:4 s:16 n:9 e:1"
"i:13 r:4 s:16 n:9 e:1"
"i:17 r:5 s:25 n:11 e:1"
"i:19 r:5 s:25 n:11 e:1"
"i:23 r:5 s:25 n:11 e:1"
"i:25 r:6 s:36 n:13 e:1"
"i:29 r:6 s:36 n:13 e:1"
"i:31 r:6 s:36 n:13 e:1"
"i:35 r:6 s:36 n:13 e:1"
"i:37 r:7 s:49 n:15 e:1"
"i:41 r:7 s:49 n:15 e:1"
"i:43 r:7 s:49 n:15 e:1"
"i:47 r:7 s:49 n:15 e:1"
"i:49 r:8 s:64 n:17 e:1"

あれ?s = r*rですよね?このコードだと、++r; s+=n; n+=2。これでrの二乗になるの?
整数の二乗数の出現間隔は、必ずある数(ここではn)を+2していった数列で表せるの?
数式で表せば(n+1)(n+1) - n*n = n + (n+1)なの?

結論から言うとそうでした。
えっこれ数学というか算数の常識だったりするんですか?Tips的な?そりゃ計算すりゃすぐ証明できましたが。

(n+1)(n+1) - n*n = n + (n+1)の証明:

左辺:n^2 + 2n +1 - n^2
        = 2n + 1
        = n + (n+1)     //

わざわざ載せるのも恥ずかしい計算...。

jの最大数定義になぜかcountが登場する

countは素数の数です。iがjで割り切れなかった時にカウントアップされる結果です。なんでこんなところにいるんですか。
しかしcountがインクリメントされるにはちゃんと条件がついているのです。if文とは違う書き方というだけです。

最大数の定義というのはbooleanを返す式で表され、ご存知のように式の結果がfalseになればループ変数は定義に従って変化していきます(今回はjが+2または+4されます)。
んでその式ですが、j <=r || !++count;となっています。ちゃんとboolを返す評価式です。ここで問題になるのは式実行順序ですね。

  • JavaScriptの||は短絡評価をするため、j <= rがtrueであれば後半は評価されない(=カウントアップもされない)
  • !++countの評価順は、インクリメント演算子が頭についているので、まずカウントアップが行われる。
  • ++countは常に0ではない正の整数であるため、!++countの結果は常にfalseである。
  • 要約:「j > rであれば++countしてから次のiに行く」

最大のjまで調べきって、かつ、ブロック内で調べているi%jの結果も0ではない。
ならばそのiは素数と断定できるわけで、そこで初めてカウントアップを行い、次のiに移るというわけです。

なお当然これもコードを読んだだけでわかっただけではありません。
ちょっと分解して全部コンソールに出しました。想像力がないなら具体例を見るしかないのです。
const MAX = 50;
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 <= s; j += (e = -e) + 3) {
      if (j <= r) {
        if (!(i % j)) {
          console.log(`i:${i} j:${j} r:${r} ccount:${count}`);
          break;
        } else {
          console.log(`i:${i} j:${j} r:${r} count:${count}`);
        }
      } else {
          console.log(`i:${i} j:${j} r:${r} count:${count}`);
          ++count;
          break;
      }
   }
}

結果

"i:5 j:5 r:3 count:2"
"i:7 j:5 r:3 count:3"
"i:11 j:5 r:4 count:4"
"i:13 j:5 r:4 count:5"
"i:17 j:5 r:5 count:6"
"i:17 j:7 r:5 count:6"
"i:19 j:5 r:5 count:7"
"i:19 j:7 r:5 count:7"
"i:23 j:5 r:5 count:8"
"i:23 j:7 r:5 count:8"
"i:25 j:5 r:6 ccount:9"
"i:29 j:5 r:6 count:9"
"i:29 j:7 r:6 count:9"
"i:31 j:5 r:6 count:10"
"i:31 j:7 r:6 count:10"
"i:35 j:5 r:6 ccount:11"
"i:37 j:5 r:7 count:11"
"i:37 j:7 r:7 count:11"
"i:37 j:11 r:7 count:11"
"i:41 j:5 r:7 count:12"
"i:41 j:7 r:7 count:12"
"i:41 j:11 r:7 count:12"
"i:43 j:5 r:7 count:13"
"i:43 j:7 r:7 count:13"
"i:43 j:11 r:7 count:13"
"i:47 j:5 r:7 count:14"
"i:47 j:7 r:7 count:14"
"i:47 j:11 r:7 count:14"
"i:49 j:5 r:8 count:15"
"i:49 j:7 r:8 ccount:15"

愚直の愚は愚かと書きます。

0
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
0
0