r-statistics-fanの日記

統計好きの現場の臨床医の覚書のようなもの

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恐るべし。