午後のひとときに、素数を列挙するプログラムを考えてみる。
素数を列挙するプログラムはいろいろと作ってきたが、今回はウィルソンの定理を使って素数判定を行うことに重きを置くこととします。
ウィルソンの定理とは、
ある2以上の自然数pにおいて、
(p－1)!≡p－1 (mod p)
を満たすとき、pは素数である。
もう一個手前で止めると、
(p－2)!≡1 (mod p)
ということも言えます。
簡単に言えば、階乗とmodを使って素数判定をするということです。
例えば、多倍長電卓LM
このソフトは、桁数を自由に設定して、C言語風にプログラムが書けるので、結構重宝しています。
例えば、この電卓を使って、100!の10進数での桁数を求めたいとして、
ceil(log10(100!))
と入力すれば、
>ceil(log10(100!))
= 158.
といった値を返してくれ、10進数で158桁必要だということが解る。
まぁ、余裕を持って200桁分表示出来るようにするには、
オプション(O)→設定→有効桁数
にて、10進表示桁数(100)の欄に200と入力して、[OK]
これで、この電卓では100!を余裕で表示出来る環境が整ったわけです。
100!と入力してみると、
>100!
=
+ 933 26215
44394 41526 81699 23885 62667 00490 71596 82643 81621 46859
29638 95217 59999 32299 15608 94146 39761 56518 28625 36979
20827 22375 82511 85210 91686 40000 00000 00000 00000 00000.
といったように、下位から5桁ずつ区切って、1行50桁で表示してくれます。
さて、今回のウィルソンの定理をつかって、100以下の素数を列挙させるプログラムを書いてみると、
for (p=2; p<=100; p++) {
if ( (p-1)!%p == p-1 ) {
printf("%d\n",p);
}
}
こんなに簡単に記述することが出来、実行結果は、
>for (p=2; p<=100; p++) {
if ( (p-1)!%p == p-1 ) {
printf("%d\n",p);
}
}
2
3
5
7
11
13
17
19
23
29
31
37
41
43
47
53
59
61
67
71
73
79
83
89
97
となる。
さて、もっと桁数を増やせば、もっと大きな数も対応出来るが、桁数を増やすと重くなってしまうので、どこまで増やせるかは環境に依存してしまうだろう。
私は高速な計算はC言語で記述することが多いので、C言語バージョンも作りたい。
#include <stdio.h>
#include <stdlib.h>
int main()
{
long long p, q, r;
for (p=2; p<100000; p++) {
q = 1;
for (r=2; r<=p-1; r++) {
q *= r;
q %= p;
}
if ( q == 1 || q == p-1 ) printf("%lld\n", p);
}
return EXIT_SUCCESS;
}
C言語で64ビット型のlong longを使ったとして、階乗計算出来る限界は20!くらいだろうから、それでは意味がないので、掛けた先から剰余を求めるようにした。
このプログラムはまったく枝刈りをしていないたたき台で、次は枝刈りの可能性を考えてみる。
#include <stdio.h>
#include <stdlib.h>
int main()
{
long long p, q, r;
printf("2\n");
for (p=3; p<100000; p+=2) {
q = 1;
for (r=2; r<=(p-2)/2; r++) {
q *= r*(p-r);
q %= p;
if ( q == 0 || q == 1 || q == p-1 ) break;
}
if ( q == 1 || q == p-1 ) printf("%lld\n", p);
}
return EXIT_SUCCESS;
}
変更点1は赤文字
変更点2は青文字
変更点3は緑文字
変更点1は、素数判別される値の枝刈りであり、誰もがやるであろう高速化です。
変更点2は、ウィルソンの定理だから出来る枝刈りです。
(p－1)!まで求めて、p－1と比較するのではなく、
(p－2)!まで求めて、1と比較するのと、たった1個しか違わないのですが、
そちらを採用すると、pが奇数しか検査しないとなると、2からp－2までの偶数個の積と考えられて、前後どちからでも計算させることが出来る。
2からp－1では奇数個になってしまうので、中央を別個考えなければ行けないので、プログラムをひとつややこしくすることになる。
ただし、これが高速化されるのかは不明ではある。
変更点3は、階乗計算の途中であっても、剰余が0になれば合成数、剰余が1かp－1になれば素数が確定すると考えていて、確率的にこれらの3つの値に早い段階でなれば、余計なループ回数が減るという算段である。
今回は外側の2つの積からやっているが、内側の2つの積のほうが速いかもしれないし、内と外も一緒にやるというのもありかもしれない。
そう思うのは、出来るだけ約数の多い値を早めに剰余を求められれば、素数か合成数かの結果も早く出ると考えている。
これに関しては、数学的な考証をする必要があるのだけれども、考証をせずに組み込んでみました。
出力結果を見比べるという検証も出来るかとは思うが、実際に正しいのかは、後々考えるとしよう。
ではでは