素数の計算（１）

 ソフィー・ジェルマン素数の続き

件の記事の中には，「p=999671で 2p+1=1999343，が約数となる。 2^p−1=2^99671−1 は10進法で300932桁となったが，開発用等ではないスペックのかなり低いパソコンでも2～3分で計算してくれた」とあったので追試しようとして，はまってしまった。

Mathematicaでは，このあたりのソフィー・ジェルマン素数を次のように一瞬で計算した。

In[1]:= Timing[ Do[p = Prime[2^16 + i];  If[Mod[p, 4] == 3 && PrimeQ[2*p + 1],  Print[p, " : ", Mod[2^p - 1, 2*p + 1]]], {i, 12900, 13000}]]
999611 : 0
999623 : 0
999671 : 0
1000151 : 0
1000211 : 0
1000403 : 0
Out[1]= {0.007209, Null}

In[2]:= Timing[Print[p = Prime[10^8], " ", PrimeQ[p]]]
2038074743 True
Out[2]= {0.000185, Null}

そこで，10^8番目の素数をpython3で計算して判定すると，

from sympy import *
import time
start_time=time.time()
p=prime(10**8)
end_time=time.time()
print(p,isprime(p),end_time-start_time)

2038074743 True 2.994921922683716

同様に，10^8番目の素数をjulia 1.5で計算して判定すると， 

using Primes
n=10^8
function sgp(n)
  p=prime(n)
  println(isprime(p)," ",p)
end

@time sgp(n)
true 2038074743
420.535971 seconds (333.45 M allocations: 4.970 GiB, 0.24% gc time)

なお，p=999671は，prime(78476)であり，この計算時間はjuliaでも0.153秒だった。それにしてもこの遅さはいったい何ということでしょう・・・？