gmpのisprimeはどこまで正確か~pseudoprime
前回の記事でgmpのisprimeがあまりにも正確というのを書いた。
デフォルトのrep=40だと、何年回してもprobably primeなのに実際は
そうでない例(pseudoprime)を見つけられないと思う。
しかしreps = を少ない数にすればpseudoprimeを見つけられそう。
rでやってみる。
library(gmp) for (re in 1:3){ i <- 1 while(TRUE){ i <- i + 2 temp <- isprime(i, reps = re) if (temp == 0 | temp == 2) next if (length(factorize(i)) > 1) break } cat(re, i, "\n") }
1 1537381
2 1943521
3 465658903
> isprime(1537381, rep = 1)
[1] 1
> isprime(1943521, rep = 2)
[1] 1
> isprime(465658903, rep = 3)
[1] 1
> factorize(1943521)
Big Integer ('bigz') object of length 3:
[1] 61 151 211
> factorize(1537381)
Big Integer ('bigz') object of length 2:
[1] 877 1753
> factorize(465658903)
Big Integer ('bigz') object of length 2:
[1] 15259 30517
rep=1でも約150万、2で約200万、3では一気に4.6億まで正しい。
恐らくrep=4以上はpseudoprimeの起こる場所を、手法を完全に
理解して限定しないと無理っぽい。4以上は総当りだと時間がかかる
のであきらめて、ググって以下を発見した。
http://www.hoegge.dk/gmp/gmp501.htm
http://www.hoegge.dk/gmp/spsp.txt
使っているプログラムが異なるが、手法は一緒っぽいので
大きく違うことはないだろう。
##引用開始
mpz_spsp1: 1,537,381 = 877*1753
mpz_spsp2: 1,943,521 = 61*151*211
mpz_spsp3: 465,658,903 = 15259*30517
mpz_spsp4: 239,626,837,621 = 346141*692281
mpz_spsp5: 2,028,576,353,203 = 1007119*2014237
mpz_spsp8: 14,910,591,535,003 = 2730439*5460877
mpz_spsp16: 76963568720428328358838771
##引用終わり
temp <- "76963568720428328358838771" factorize(temp) isprime(temp, rep =16) isprime(temp, rep =17) > temp <- "76963568720428328358838771" > factorize(temp) Big Integer ('bigz') object of length 2: [1] 12406737582493 6203368791247 > isprime(temp, rep =16) [1] 1 > isprime(temp, rep =17) [1] 0
rep=16でも、すでに素のRで桁落ちせず数値として認識できる範囲を超えている。
(factorizeなどは数値ではなく”文字列”として与えないと計算してくれなかった。)
デフォルトのrep=40では、想像もつかない数まで正確ということになる。
isprimeが想像以上に正確なのにあらためて驚いた。
Miller-Rabin probabilistic primary tests恐るべし。