１．はじめに

或る会合に２４人の人たちが参加しました。
この２４人の中に、誕生日が同じ日である人が１組いる確率はどれくらいか？

答え
一人の人物について、特定の日の生まれである確率は1/365です。
２人目の人が誕生日が同じでない確率は364/365=1-1/365になります。
同様に、３人目の人が前の２人と誕生日が同じでない確率は363/365=1-2/365になります。
４人目の人が前の３人と誕生日が同じでない確率は362/365=1-3/365になります。
従って、２４人の誕生日が同じでない確率ｐは

p=364/365*363/365*362/365*.....342/365
 =(1-1/365)*(1-2/365)*(1-3/365)*(1-23/365)
 ≒0.4616

従って、誕生日が同じである人が１組はいる確率は1-0.4616=0.5384　なので
２４人の中に５０％くらいの確率で、誕生日が同じである人が１組はいることになります。
どうでしょうか、直感的には納得できませんが、これが真実なのです！
これが有名なBirtday Paradox　です。
これを素因数分解に応用したのがPollardという人でモンテカルロ法と呼ばれています。
Birtday Paradoxを一般化すると次の様になります。

n通りの値があり、m回試行して同じ値が起きない確率は次のようになります。
確率をpとすると
p=(1-1/n)(1-2/n)(1-3/n).....(1-m/n)≒e^(-m^2/(2n))
e^(-m^2/(2n))=p

従って、m=sqrt(2nlog(1/p))

特に確率がp=1/2の時はm=sqrt(2nlog(2))　になるから、m≒1.1774sqrt(n)　となります。
同様にp=1/3の時はm=sqrt(2nlog(3))=1.4823sqrt(n)
p=1/4の時はm=sqrt(2nlog(4))=1.6651sqrt(n) となります。
どちらにしても確率はsqrt(n)に比例する事が判ります。

さて、次の様な数列を考えましょう。
a(1)=2
a(i)=a(i-1)^2+1 mod n
a(2i)=(a(2i-2)^2+1 mod n)^2+1 mod n

a(i)とa(2i)が同じ値をとる時、GCD(a(2i)-a(i),n)=1,nでなければnの素因数が見つかった事になる。

例１

n=103*1009=103927　として a(i),a(2i)を求めてみると次の様に１４回目にGCD(a(2i)-a(i),n)=103となって
n=103927の素因数103が見つかった。

  i    a(i)   a(2i)

  1      2      5 
  2      5    677
  3     26  94852
  4    677  49802
  5  42622 102606
  6  94852   7364
  7  45442  33175
  8  49802  40303
  9  21350  60829
 10 102606   4212
 11  82210  31674
 12   7364  80667
 13  82530  70297
 14  33175  64487

a(i),a(2i)を mod 103 で考えると次の様な結果になり、１４回目でa(i)=a(2i)=9となっている事がわかる。
この事実がa(i),a(2i)を mod 103927で考えた時、１４回目で素因数103が見つかった理由である。

 　i    a(i)   a(2i)

  1      2      5
  2      5     59
  3     26     92
  4     59     53
  5     83     18
  6     92     51
  7     19      9
  8     53     30
  9     29     59
 10     18     92
 11     16     53
 12     51     18
 13     27     51
 14      9      9


この結果をまとめると、或るｎをn=p*qと素因数分解する時、分解に要する試行回数はBirtday Paradoxに
よってsqrt(n)に比例するが、一方nの素因数pの試行回数はsqrt(p)に比例するので、結局nの分解に要する
試行回数はn^(1/4)に比例することが判ったことになります。

また、試行回数はnの大きさよりも素因数であるpの大きさ（試行回数）に依存することがわかります。
すなわち、n1=p1*p2,n2=q1*q2 として n1>n2,p1p1 とするとn1はn2より少ない試行回数で
素因数p1が見つかると予想されます。
計算例で、2^256+1が2^512+1,2^1024+1 に比べて試行回数が多いことが示されています。



       
UBASICで実際にプログラムを作りました、次の様になります。

　　プログラム名：moncarno.ub

   10   input "N=";N:input "seed=";seed
   20   a1=seed:a2=(a1^2+1)@N
   30   G=1:cnt=0
   40   while G<=1
   45     cnt=cnt+1:if cnt@10000 = 0 then print "cnt=";cnt
   50     W=1
   60     for I=1 to 10
   70       a1=(a1^2+1)@N
   80       a2=((a2^2@N+1)^2+1)@N
   90       W=W*(a1-a2)@N
  100     next
  110     G=gcd(W,N)
  120   wend
  130   if G > 1 and G < N then print "factor is ";G
  140   :else print "try new seed" 
  150   if modpow(2,G-1,G) <> 1 then print G," is not prime"
　160   :else print G," is maybe prime"
  200   end



seedは通常２で良いでしょう。２で分解できない場合は３、４、５、・・と変えて試すと良いでしょう。





3.計算例

N=2^64+1= 18446744073709551617 とすると

factor is  274177 	 cnt= 810 
試行回数８１０回で分解。
実際は、N=274177 * 67280421310721 
この数は、１８８０年Landryが８２歳の時に数ヶ月かけて、素因数274177を発見しました。
パソコンも電卓も無い時代に、さぞや大変だったと思います。
今や１秒もかからないで分解されてしまいます。


N=2^512+1
=13407807929942597099574024998205846127479365820592393377723561443721764030073546976801874298166903427690031858186486050853753882811946569946433649006084097 

factor is  2424833 	 cnt= 1570 
最終的には、N=2424833 * 7455602825647884208337395736200454918783366342657 * P99  


 
N=2^1024+1
=179769313486231590772930519078902473361797697894230657273430081157732675805500963132708477322407536021120113879871393357658789768814416622492847430639474124377767893424865485276302219601246094119453082952085005768838150682342462881473913110540827237163350510684586298239947245938479716304835356329624224137217 

factor is  45592577 	 cnt= 9010 
最終的には、N=45592577 * 6487031809 * 4659775785220018543264560743076778192897 * P252 


N=2^256+1
 = 115792089237316195423570985008687907853269984665640564039457584007913129639937 
factor is  1238926361552897 	 cnt= 14817000 

最終的には、N=1238926361552897 * P62 
2^512+1,2^1024+1 に比べて試行回数が多いのは、素因数が大きいのが原因であることが判るでしょう。



N=2^37-1 = 137438953471 
factor is  223 	 cnt= 10 


N=2^41-1= 2199023255551 
factor is  13367 	 cnt= 50



N=2^83-1= 9671406556917033397649407 
factor is  167 	 cnt= 10 


N= 11111111111111111 
factor is  2071723 	 cnt= 2870 

パソコンの性能にもよりますが、よほど運が悪くない限り数秒で素因数が見つかります。
例えば、2^1024+1　は３００桁位の数ですが素因数が45592577と非常に小さいため数秒で
見つかりました。
反対に、2^256+1　の場合は素因数が1238926361552897と大きいため数分かかりました。

この様に素因数が小さい場合は、試行回数は少なくて済みますが一般的には分解されるまで
素因数の大きさは判らない訳ですから、処理時間を考えるとＮは２０桁から２５桁位が限度
ではないでしょうか。