音楽室と化学室と美術室とPC室の融合部屋 所謂自由室

趣味と気分で適当に色々やります.なんかあるとたまに更新します.

プログラミング初心者でもフェルマーテストをしてみたい

皆さんフェルマーの小定理って知ってますよね???
知ってると思います!はい!!
ざっくりいうと、とある数naがあり、gcd(n,a) = 1となる時
n素数の時、必ず

a^{n-1} \equiv 1 \pmod{n}
a^{n} \equiv a \pmod{n}

が成り立つというものです。
逆に言うと、これが成り立たなかったら素数ではないということですね。
ここで、a2以上n-1以下の整数を適当に選んで計算してやればいいわけですね。
これで素数判定が簡単に!!って言おうとしてもそうは問屋がおろしません・・・
どんなaを選んでも通り抜けてしまう数値があるみたいで、『カーマイケル数』というみたいです。

561や1729がそれで、結構やねこい存在みたいです。(正直言うと個人的には結構おもしろそうだなーって感じの対象ですが)

 

で、今回このフェルマーテストをC++で組んでみようかなーと思いました。
どうやれば高速でできるかなーって考えた結果、二進法を使おうと思いました。

a^n乗を計算してaに等しかったらいいのですね・・・
そこでnを2進法で表します。
10だと(1010)_2ですね。 つまり2^3 + 2^1ですね
つまり、a^{2^3} * a^{2^1}を計算すればいいので
a = a^1
(a)^2 = a^2
(a^2)^2 = a^4
(a^4)^2 = a^8
...

つまり、r=1から開始して、nを2で割ってやって余りが1ならra^{2^k}をかけてやる。
で、nが0に到達すれば終了 お疲れ様です! って感じですね!!
1の位から逐次的に計算していくので、log_2 n程度の計算量で済む感じですね。凄い。
あと今回は気休め程度にanのgcdを計算し、互いに素でなければ素数判定から弾くようにはしてます。

<hr size="1" />


ちなみに乱数生成のMT.hは

C言語による乱数生成
から取ってきました。Cのrandomライブラリわっかんない・・・;

参考記事

素数判定いろいろ - フェルマーテスト・ミラーラビン素数判定法 - Qiita

 以下ソース

#include < iostream>
#include < time.h>
#include < fstream>
#include < string>
#include < chrono>
#include "MT.h"
#define TESTMAX 100
#define NUM 100000000

unsigned long long Rand(unsigned long long A, unsigned long long B) { return genrand_int32() % (B - A + 1) + A; }

unsigned long long gcd(unsigned long long a, unsigned long long b) { return (a%b==0)?b:gcd(b, a%b); }

unsigned long long FelmatL(unsigned long long a, unsigned long long num) {
unsigned long long result = 1, tmp = a, tnum = num;
while (tnum>0) {
if (tnum % 2 == 1)result = (result*tmp) % num;
tmp = (tmp*tmp) % num;
tnum >>= 1;
}
return result;
}

bool FermatTest(unsigned long long num) {
for (unsigned long long cnt = 0; cnt < TESTMAX; cnt++) {
unsigned long long a = Rand(2, num - 1);
if (gcd(a, num) != 1)return false;
if (FelmatL(a, num) != a)return false;
}
return true;
}

int main(void) {
init_genrand*1;
std::ofstream outf("Prime?.txt");
for (unsigned long long i = 3; i < NUM; i+=2) {
if *2outf << std::to_string(i) << std::endl;
}
outf.close();
std::cout << "Finish" << std::endl;
return 0;
}




出力(全部張ったらフリーズしたので2000までにしておきます)


3
5
7
11
13
17
19
23
29
31
37
41
43
47
53
59
61
67
71
73
79
83
89
97
101
103
107
109
113
127
131
137
139
149
151
157
163
167
173
179
181
191
193
197
199
211
223
227
229
233
239
241
251
257
263
269
271
277
281
283
293
307
311
313
317
331
337
347
349
353
359
367
373
379
383
389
397
401
409
419
421
431
433
439
443
449
457
461
463
467
479
487
491
499
503
509
521
523
541
547
557
563
569
571
577
587
593
599
601
607
613
617
619
631
641
643
647
653
659
661
673
677
683
691
701
709
719
727
733
739
743
751
757
761
769
773
787
797
809
811
821
823
827
829
839
853
857
859
863
877
881
883
887
907
911
919
929
937
941
947
953
967
971
977
983
991
997
1009
1013
1019
1021
1031
1033
1039
1049
1051
1061
1063
1069
1087
1091
1093
1097
1103
1109
1117
1123
1129
1151
1153
1163
1171
1181
1187
1193
1201
1213
1217
1223
1229
1231
1237
1249
1259
1277
1279
1283
1289
1291
1297
1301
1303
1307
1319
1321
1327
1361
1367
1373
1381
1399
1409
1423
1427
1429
1433
1439
1447
1451
1453
1459
1471
1481
1483
1487
1489
1493
1499
1511
1523
1531
1543
1549
1553
1559
1567
1571
1579
1583
1597
1601
1607
1609
1613
1619
1621
1627
1637
1657
1663
1667
1669
1693
1697
1699
1709
1721
1723
1733
1741
1747
1753
1759
1777
1783
1787
1789
1801
1811
1823
1831
1847
1861
1867
1871
1873
1877
1879
1889
1901
1907
1913
1931
1933
1949
1951
1973
1979
1987
1993
1997
1999

*1:unsigned)time(NULL

*2:i-1) % (NUM/1000) == 0)std::cout << i << std::endl;
if (FermatTest(i