「プログラマにも理解できるように問題が提示されるように望む」 (ろうそくの科学より)。 というわけで、ガウスジョルダンが理解できず、数一を2度落した男が贈ります。

メルセンヌ素数 

  メルセンヌ素数とは、 2^P-1 の形の素数です。(Pは自然数)
2進数で示すと全て 1 の立った形になります。
メルセンヌ素数プログラムを自作したい方へ 

新メルセンヌ素数を発見できる速度のプログラムを作りたい場合は、
Richard Crandallさん作成のlucdwt.c
が参考になります。

http://www.perfsci.com/free-software.asp#giantint
にはそれ以外にも楕円曲線法やストラッセン乗算法など巨大数を扱う
Cソースがいろいろあります。
Linux上でUBASICを使う 

linux に DOSBOX をインストールし DOSBOX を起動
$dosbox
ubasicのある場所が/tmpなら/tmpをC:ドライブに指定
>mount c /tmp
ubasic を起動
>C:ubv32.exe
プログラムを /tmp/TEST.UB という名前で用意しておいて load
>load "test"
実行
>run

フーリエ変換を使った多倍長乗算 

  4 要素のフーリエ変換を使って多倍長の乗算をしてみます。

 4 要素なので 4 乗した時に 1 になる 4 乗根 (0,1) を使います。
 (8 要素の場合は当然 8 乗根を使います)

 12 の 2 乗を計算してみます。

 まず 12 を多倍長の表現に変換します。

        ( 0, 0)         ( 0, 0)         ( 1, 0)         ( 2, 0)

次に 4 要素の順フーリエ変換を行います。

( 0, 1)*( 0, 0)+(-1, 0)*( 0, 0)+( 0,-1)*( 1, 0)+( 1, 0)*( 2, 0) =( 2,-1)
(-1, 0)*( 0, 0)+( 1, 0)*( 0, 0)+(-1, 0)*( 1, 0)+( 1, 0)*( 2, 0) =( 1, 0)   
( 0,-1)*( 0, 0)+(-1, 0)*( 0, 0)+( 0, 1)*( 1, 0)+( 1, 0)*( 2, 0) =( 2, 1)   
( 1, 0)*( 0, 0)+( 1, 0)*( 0, 0)+( 1, 0)*( 1, 0)+( 1, 0)*( 2, 0) =( 3, 0)

次に各要素を 2 乗して、 4 要素なので正規化する為に 4 で割ります。

( 2,-1)*( 2,-1)/4=( 0.75,-1)
( 1, 0)*( 1, 0)/4=( 0.25, 0)
( 2, 1)*( 2, 1)/4=( 0.75, 1)
( 3, 0)*( 3, 0)/4=( 2.25, 0)

次に 4 要素の逆フーリエ変換を行います。
( 0,-1)*( 0.75,-1)+(-1, 0)*( 0.25, 0)+( 0, 1)*( 0.75, 1)+( 1, 0)*( 2.25, 0)=( 0, 0)
(-1, 0)*( 0.75,-1)+( 1, 0)*( 0.25, 0)+(-1, 0)*( 0.75, 1)+( 1, 0)*( 2.25, 0)=( 1, 0)
( 0, 1)*( 0.75,-1)+(-1, 0)*( 0.25, 0)+( 0,-1)*( 0.75, 1)+( 1, 0)*( 2.25, 0)=( 4, 0)
( 1, 0)*( 0.75,-1)+( 1, 0)*( 0.25, 0)+( 1, 0)*( 0.75, 1)+( 1, 0)*( 2.25, 0)=( 4, 0)

これで回答、

        ( 0, 0)         ( 1, 0)         ( 4, 0)         ( 4, 0)

が求まりました。実際の計算で要素数が増えてくると誤差が出て来てちょうど整数には
ならないので四捨五入してやります。(あまりにも誤差が大きい場合は要素数を増やす)
フーリエ変換には FFT が使えるので桁が増えた場合には素直に計算する場合に比較し
て大変高速になります。この方法は東大計算機センターの金田さんによりπの世界記録
をだすのにも使用されています。

ちょっと補足

( 1, 0) は、複素数の実部が 1 虚部が 0 の数を表しています。

( 0, 1)*(  ,  )+(-1, 0)*(  ,  )+( 0,-1)*(  ,  )+( 1, 0)*(  ,  ) 
(-1, 0)*(  ,  )+( 1, 0)*(  ,  )+(-1, 0)*(  ,  )+( 1, 0)*(  ,  ) 
( 0,-1)*(  ,  )+(-1, 0)*(  ,  )+( 0, 1)*(  ,  )+( 1, 0)*(  ,  ) 
( 1, 0)*(  ,  )+( 1, 0)*(  ,  )+( 1, 0)*(  ,  )+( 1, 0)*(  ,  ) 

ここに出て来る定数は回転子と
いいます。右 中央の 4 乗根の (0,1)  から始めて 上側にその 2 乗 3 乗 と続き
それぞれを 左側に 2 乗 3 乗 4 乗した数になっています。逆フーリエ変換の時は
 (0,1) の逆数の (0,-1) から始めます。
  
  実際は FFT のプログラムを流用してコーディングされると思うのですが世の中に
流通している FFT は正規化をしないものした後の値を返却するものなどいろいろで
す、入手したプログラムとここに書いた値の差を見て多倍長乗算のプログラムを作
成してください。

８要素の FT の回転子は

    0              1              2            3             4            5               6           7
0(-0.57, 0.57) ( 0   , 1   ) (-0.57, 0.57) (-1   , 0   ) (-0.57,-0.57) ( 0   ,-1   ) ( 0.57,-0.57) ( 1   , 0   )
1( 0   , 1   ) (-1   , 0   ) ( 0   ,-1   ) ( 1   , 0   ) ( 0   , 1   ) (-1   , 0   ) ( 0   ,-1   ) ( 1   , 0   )
2(-0.57, 0.57) ( 0   ,-1   ) ( 0.57, 0.57) (-1   , 0   ) ( 0.57,-0.57) ( 0   , 1   ) (-0.57,-0.57) ( 1   , 0   )
3(-1   , 0   ) ( 1   , 0   ) (-1   , 0   ) ( 1   , 0   ) (-1   , 0   ) ( 1   , 0   ) (-1   , 0   ) ( 1   , 0   )
4(-0.57,-0.57) ( 0   , 1   ) ( 0.57,-0.57) (-1   , 0   ) ( 0.57, 0.57) ( 0   ,-1   ) (-0.57, 0.57) ( 1   , 0   )
5( 0   ,-1   ) (-1   , 0   ) ( 0   , 1   ) ( 1   , 0   ) ( 0   ,-1   ) (-1   , 0   ) ( 0   , 1   ) ( 1   , 0   )
6( 0.57,-0.57) ( 0   ,-1   ) (-0.57,-0.57) (-1   , 0   ) (-0.57, 0.57) ( 0   , 1   ) ( 0.57, 0.57) ( 1   , 0   )
7( 1   , 0   ) ( 1   , 0   ) ( 1   , 0   ) ( 1   , 0   ) ( 1   , 0   ) ( 1   , 0   ) ( 1   , 0   ) ( 1   , 0   )

 FFT 風に書くと

a(0)    t(0)=a(0)+a(4)    a(0)=t(0)+      t(4)    t(0)=a(0)+             a(4)    a(0)=t(0)                                                                         
a(1)    t(1)=a(0)-a(4)    a(1)=t(0)-      t(4)    t(1)=a(0)-             a(4)    a(1)=t(4)                                                                         
a(2)    t(2)=a(1)+a(5)    a(2)=t(1)+(0,1)*t(5)    t(2)=a(1)+( 0   ,1   )*a(5)    a(2)=t(2)                                                                         
a(3) → t(3)=a(1)-a(5) → a(3)=t(1)-(0,1)*t(5) → t(3)=a(1)-( 0   ,1   )*a(5) → a(3)=t(6)                                                                         
a(4)    t(4)=a(2)+a(6)    a(4)=t(2)+      t(6)    t(4)=a(2)+( 0.57,0.57)*a(6)    a(4)=t(1)                                                                         
a(5)    t(5)=a(2)-a(6)    a(5)=t(2)-      t(6)    t(5)=a(2)-( 0.57,0.57)*a(6)    a(5)=t(5)                                                                         
a(6)    t(6)=a(3)+a(7)    a(6)=t(3)+(0,1)*t(7)    t(6)=a(3)+(-0.57,0.57)*a(7)    a(6)=t(3)                                                                         
a(7)    t(7)=a(3)-a(7)    a(7)=t(3)-(0,1)*t(7)    t(7)=a(3)-(-0.57,0.57)*a(7)    a(7)=t(7)

代数体を使った多倍長乗算 

  フーリエ変換と同じように代数体を使っても多倍長の乗算が
できます。

   17 による剰余演算を使って 4 要素の多倍長乗算をしてみます。
 4 要素なので 17 を法とした演算で 4 乗した時に 1 になる 
4 乗根 4 を使います。

 これも 12 の 2 乗を計算してみます。

 まず 12 を多倍長の表現に変換します。

        0         0          1         2

次に 代数体上で 4 要素の順変換を行います。
(全て 17 を法とする剰余演算であることに注意)

  4  *  0 + 16 * 0 +13 * 1 + 1 * 2  = 15 
 16  *  0 +  1 * 0 +16 * 1 + 1 * 2  =  1    
 13  *  0 + 16 * 0 + 4 * 1 + 1 * 2  =  6 
  1  *  0 +  1 * 0 + 1 * 1 + 1 * 2  =  3 

次に各要素を 2 乗して、 4 要素なので正規化する
為に 代数体上の 4 の逆数である 13 をかけます。
(除算とはかけて 1 になる数をかける事と同値)

 15 *15 * 13 =  1
  1 * 1 * 13 = 13
  6 * 6 * 13 =  9
  3 * 3 * 13 = 15

次に 代数体上で 4 要素の逆の変換を行います。

 13  *  1 + 16 *13 + 4 * 9 + 1 * 15  =  0 
 16  *  1 +  1 *13 +16 * 9 + 1 * 15  =  1 
  4  *  1 + 16 *13 +13 * 9 + 1 * 15  =  4 
  1  *  1 +  1 *13 + 1 * 9 + 1 * 15  =  4 

これで回答

     0  1  4  4

が求まりました。基数が10進だろうと5進だろうと 12*12 は
 144 だということがわかります。代数体の演算にも FFT と
同じような演算量節約法が使えますし、法を 2^n+1 の素数
にとることによりシフトと加減算で簡単に剰余演算を行うこ
とができます。(法は素数な必要があります)

代数体を使った多倍長乗算の演算の節約 

  代数体を使った多倍長演算にどのように FFT のような演算の節約が適用
できるか考えます。

   今度は 17 による剰余演算を使って 8 要素の多倍長乗算によって考えてみます。
 8 要素なので 17 を法とした演算で 8 乗した時に 1 になる 
8 乗根 2 を使います。

 これも 12 の 2 乗を計算してみます。

 まず 12 を多倍長の表現に変換します。

      0        0        0        0        0        0        1        2

 次に 代数体上で 8 要素の正の変換を行います。
 (全て 17 を法とする剰余演算であることに注意)

 0  2 * 0 +  4 * 0 +  8 * 0 + 16 * 0 + 15 * 0 + 13 * 0 +  9 * 1 +  1 * 2  = 11 
 1  4 * 0 + 16 * 0 + 13 * 0 +  1 * 0 +  4 * 0 + 16 * 0 + 13 * 1 +  1 * 2  = 15 
 2  8 * 0 + 13 * 0 +  2 * 0 + 16 * 0 +  9 * 0 +  4 * 0 + 15 * 1 +  1 * 2  =  0 
 3 16 * 0 +  1 * 0 + 16 * 0 +  1 * 0 + 16 * 0 +  1 * 0 + 16 * 1 +  1 * 2  =  1 
 4 15 * 0 +  4 * 0 +  9 * 0 + 16 * 0 +  2 * 0 + 13 * 0 +  8 * 1 +  1 * 2  = 10 
 5 13 * 0 + 16 * 0 +  4 * 0 +  1 * 0 + 13 * 0 + 16 * 0 +  4 * 1 +  1 * 2  =  6 
 6  9 * 0 + 13 * 0 + 15 * 0 + 16 * 0 +  8 * 0 +  4 * 0 +  2 * 1 +  1 * 2  =  4 
 7  1 * 0 +  1 * 0 +  1 * 0 +  1 * 0 +  1 * 0 +  1 * 0 +  1 * 1 +  1 * 2  =  3 

 これを FFT 風にしてみると

a(0)    t(0)=8*a(0)+ 9*a(4)    a(0)=4*t(0)+13*t(4)    t(0)=1*a(0)+16*a(4)    a(0)=t(0)                                                                         
a(1)    t(1)=  a(4)+   a(0)    a(1)=  t(4)+   t(0)    t(1)=  a(4)+   a(0)    a(1)=t(4)                                                                         
a(2)    t(2)=4*a(1)+13*a(5)    a(2)=4*t(1)+13*t(5)    t(2)=1*a(1)+16*a(5)    a(2)=t(2)                                                                         
a(3) → t(3)=  a(5)+   a(1) → a(3)=  t(5)+   t(1) → t(3)=  a(5)+   a(1) → a(3)=t(6)                                                                         
a(4)    t(4)=2*a(2)+15*a(6)    a(4)=1*t(2)+16*t(6)    t(4)=1*a(2)+16*a(6)    a(4)=t(1)                                                                         
a(5)    t(5)=  a(6)+   a(2)    a(5)=  t(6)+   t(2)    t(5)=  a(6)+   a(2)    a(5)=t(5)                                                                         
a(6)    t(6)=1*a(3)+16*a(7)    a(6)=1*t(3)+16*t(7)    t(6)=1*a(3)+16*a(7)    a(6)=t(3)                                                                         
a(7)    t(7)=  a(7)+   a(3)    a(7)=  t(7)+   t(3)    t(7)=  a(7)+   a(3)    a(7)=t(7)                                                                         

  これで上の計算と同じ結果が得られます。かなりの演算量の節約になっています。
  これでプログラム(FORTRAN77!)を書くと

        integer ia(8)
        integer iw1(12),iw2(12)
        data ia / 1, 2, 0, 0, 0, 0, 0, 0/
        n=17
c       n=257
        write(*,3) ia 
        call initfft(iw1,iw2,n,1)
        write(*,1) iw1
        write(*,2) iw2
        call fft(ia,iw1,iw2,n)
        do i=1,8
         ia(i)=mod(ia(i)*ia(i),n)
c 8 で割る(8 を掛けて 1 になる数を掛ける)
         ia(i)=mod(ia(i)*15,n)
c        ia(i)=mod(ia(i)*225,n)
        enddo
        call initfft(iw1,iw2,n,-1)
        write(*,1) iw1
        write(*,2) iw2
        call fft(ia,iw1,iw2,n)
        write(*,4) ia 
 1      format('iw1   =',12i4)
 2      format('iw2   =',12i4)
 3      format('input =', 8i4)
 4      format('output=', 8i4)
        end
        subroutine fft(ia,iw1,iw2,n)
        integer ia(8),it(8)
        integer iw1(12),iw2(12)
        i1=1
        do j=1,3
         do i=1,4
          it(i*2-1)=mod(ia(i)+ia(i+4),n)
          itt=mod(iw1(i1)*ia(i),n)+mod(iw2(i1)*ia(i+4),n)
          it(i*2)=mod(itt,n)
          i1=i1+1
         enddo
         do i=1,8
          ia(i)=it(i)
         enddo
        enddo
        ia(1)=it(1)
        ia(2)=it(5)
        ia(3)=it(3)
        ia(4)=it(7)
        ia(5)=it(2)
        ia(6)=it(6)
        ia(7)=it(4)
        ia(8)=it(8)
        end
        subroutine initfft(iw1,iw2,n,if)
        integer iw1(12),iw2(12)
        i1=1
        i2=n-1
c 原始 8 乗根 (逆変換の時は逆数)
        i3=2
        if(if.eq.-1) i3=9
c       i3=4   
c       if(if.eq.-1) i3=193
        i5=1
        i6=4
        i9=1
        do j=1,3
         i10=i1
         i11=i2
         do i7=1,i6
          do i8=1,i5
            iw1(i9)=i10
            iw2(i9)=i11
            i9=i9+1
          enddo
          i10=mod(i10*i3,n)
          i11=mod(i11*i3,n)
         enddo
         i5=i5*2
         i6=i6/2
         i3=mod(i3*i3,n)
        enddo
        end 

  となります。(実際は mod の計算など全てが高速になるようにプログラム
しなければいけません)

  このプログラムで答えがあうのは、桁が 8 桁を越えないのはもちろん
ですが、筆算で書いたときの中間和が 17 以上にならない必要もあります。
  例えば

       3  3
    *  3  3
---------------
       9  9
    9  9
---------------
    9 18  9 
      18 で 17 以上なのでこの場合は答えがあいません。( 9 1 9 となる) 

おまけ 17 を 法とする 九九表

     |  1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16
  ---+---------------------------------------------------------------
   1 |  1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16
   2 |  2   4   6   8  10  12  14  16   1   3   5   7   9  11  13  15
   3 |  3   6   9  12  15   1   4   7  10  13  16   2   5   8  11  14
   4 |  4   8  12  16   3   7  11  15   2   6  10  14   1   5   9  13
   5 |  5  10  15   3   8  13   1   6  11  16   4   9  14   2   7  12
   6 |  6  12   1   7  13   2   8  14   3   9  15   4  10  16   5  11
   7 |  7  14   4  11   1   8  15   5  12   2   9  16   6  13   3  10
   8 |  8  16   7  15   6  14   5  13   4  12   3  11   2  10   1   9
   9 |  9   1  10   2  11   3  12   4  13   5  14   6  15   7  16   8
  10 | 10   3  13   6  16   9   2  12   5  15   8   1  11   4  14   7
  11 | 11   5  16  10   4  15   9   3  14   8   2  13   7   1  12   6
  12 | 12   7   2  14   9   4  16  11   6   1  13   8   3  15  10   5
  13 | 13   9   5   1  14  10   6   2  15  11   7   3  16  12   8   4
  14 | 14  11   8   5   2  16  13  10   7   4   1  15  12   9   6   3
  15 | 15  13  11   9   7   5   3   1  16  14  12  10   8   6   4   2
  16 | 16  15  14  13  12  11  10   9   8   7   6   5   4   3   2   1

上の変換はさらに

a(0)    t(0)=8*(a(0)-a(4))   a(0)=4*(t(0)-t(4))   t(0)=1*(a(0)-a(4))   a(0)=t(0)                                                                         
a(1)    t(1)=   a(4)+a(0)    a(1)=   t(4)+t(0)    t(1)=   a(4)+a(0)    a(1)=t(4)                                                                         
a(2)    t(2)=4*(a(1)-a(5))   a(2)=4*(t(1)-t(5))   t(2)=1*(a(1)-a(5))   a(2)=t(2)                                                                         
a(3) → t(3)=   a(5)+a(1) → a(3)=   t(5)+t(1) → t(3)=   a(5)+a(1) → a(3)=t(6)                                                                         
a(4)    t(4)=2*(a(2)-a(6))   a(4)=1*(t(2)-t(6))   t(4)=1*(a(2)-a(6))   a(4)=t(1)                                                                         
a(5)    t(5)=   a(6)+a(2)    a(5)=   t(6)+t(2)    t(5)=   a(6)+a(2)    a(5)=t(5)                                                                         
a(6)    t(6)=1*(a(3)-a(7))   a(6)=1*(t(3)-t(7))   t(6)=1*(a(3)-a(7))   a(6)=t(3)                                                                         
a(7)    t(7)=   a(7)+a(3)    a(7)=   t(7)+t(3)    t(7)=   a(7)+a(3)    a(7)=t(7)                                                                         

と書けそうです。

2^P+1 による乗余の計算 

  2 進法の計算機では、2^P による乗余の計算は  2^P-1 との and だけで
実現できます。
  2^P+1 による乗余は、2^P-1 との and から 2^P で 割った結果(シフト
で実現できます)を引くことで実現できます。結果が負になった場合は、2^P+1
を加算します。

2^(2^n)+1 の形の数 

  2^(2^n)+1 の形の数はフェルマー数と呼ばれて Ｆn で表されます。
Ｆ1 〜 Ｆ4 は素数ですがそれ以上の数に素数は見つかっていません。

2^31+1 の n 乗根 

  2^32+1 ぐらいが法に適当なのですが上に書いたように 2^32+1 は素数
ではないので大きな n に対する n 乗根を見つけることができません。
F5 = 232+1 = 4294967297 = 641 x 6700417
  なので、641 と 6700417 を法とする演算で調べてみましたが、128 乗
根が最大で 256 乗根は存在しないようです。

Ｆ5以上の数を法に使うには 

  上に書いたように 2^32+1 以上では適当な n 乗根を見つけることはで
きません。そこで再帰的に演算を適用することを考えてみましょう。
  例えば 2^1048576 位の演算をしようと思えば 2^1024+1 を法にして
1024 要素の順変換(2^1024+1 の 1024 乗根はすぐに見つかります)をす
ると 2^1024 の大きさの要素が 1024 個できます。それぞれに今度は 
2^32+1 を法として 32 要素の順変換をします。すると 2^32 の大きさ
の要素が 1024*32 個できます。この要素それぞれを乗算して 逆変換をす
ると乗算結果が求まります。この時に各演算の回転子は 2^(2^n)+1 を
法にする限り 2 の巾になるので乗算はシフトだけで済み命令で直接扱え
る長さを越えるような乗算も高速に可能です。

http://www.craig-wood.com/nick/armprime/math.html
には、2^64-2^32+1 を法に使った整数型 DWT を作成された方の方法があります。

また別の方法 

  今度は再帰的に演算を適用するのではなくて、他の 2^(2^n)+1 とお
互いに素(最大公約数が 1) な法で 同じ要素数で順変換します。2 つの法
はお互いに素なので 2 つの結果から 2 の法の積までの結果が求められる
はずです。

後先生 

http://www.cs.t-kougei.ac.jp/nsim/index.html

　円周率の世界記録で有名な後先生の研究室です。
FFT を多重に適用することによって浮動小数点の精度を節約する方法など
役にたちそうです。

n 乗根の存在条件 

  ある素数 p を法にした場合の原始 n 乗根が存在する必要十分条件は、
       p = 1 mod n 
です。

2^p-1 を割る整数について 

  2^p-1(p は素数) の形をもつ整数でメルセンヌ素数にならない合成数は
q = +/-1 (mod 8) and q = 2kp + 1 の形の約数しかもちません。

  2kp+1 で、2^p-1 が割り切れるかを確かめる場合は直接除算をせずに
2kp+1 を法とする剰余演算で 2^p-1 を計算して 0 になるかを確かめる。
  2kp+1 を 3,5 などの小さな素数で割って合成数を排除するのも有効です。
この方法をルーカスチェックの前にかけることで高速化できます。

また、p が 4 で割ってあまりが 3 になり、 2*p+1 が素数の場合は、
2*p+1 で 2^p+1 が割り切れる。

固定小数点を使ったFFT 

  FFT の回転子の実部、虚部は、-1〜1 の範囲の数なので固定小数点を
使った演算が可能です。MMX命令が使えるかもしれません。(mersenne.org でも
K6用のプログラムは固定小数点を使っているそうです。)

剰余の演算 

  素数を法に使った計算では剰余の演算速度が全体の性能に大きく効いて
きます。こんなことが「計算の効率化とその限界」に書いてありました。
        p = 2^l - d    d が 2^l  より十分に小さい
場合に x が 2^l 以上の場合に 2^l を越える分を xu 未満を xl と
すると 2^l = d mod p から
        x = d*xu + xl mod p
が成立ちます。x がまだ p より大きい場合はこの計算を繰り返します。
d^2 が p 未満 の場合には 3 回で収束することが知られています。
  テストプログラムです。
        integer p,l,d,x,xu,xl,xa
        p=11213  
        l=14
        d=2**l-p
        x=2**17-1
        write(*,*) ' p=',p,' l=',l,' d=',d,' d**2=',d*d
        write(*,*) ' mod(x,p)= ',mod(x,p)
        xa=x
        do i=1,7
         xu=rshift(xa,l)
         xl=and(xa,2**l-1)
         xa=d*xu + xl
         write(*,*) ' xu= ',xu,' xl= ',xl,' xa= ',xa
        enddo
        end 
 d を小さくするために p と l の差が小さい必要があるので p には
2 の冪乗になるべく近い素数を選ぶ必要があります。

  2^16 以下で 128 乗根を持つなるべく大きい素数を探すと
        i=2**16+1
 1      continue
        i=i-128   
        do j=3,2**8,2
          if(mod(i,j).eq.0) then
            write(*,*) i,' devide by ',j 
            goto 1
          endif
        enddo
        write(*,*) ' prime=  ',i
        end
21164:~#a.out
 65409 can devide  3
 65281 can devide  97
 65153 can devide  11
 65025 can devide  3
 64897 can devide  7
 64769 can devide  239
 64641 can devide  3
  prime=   64513
がみつかりました。

小数点以下を扱う為には 

  ここまでは、整数のみを扱ってきましたが小数点以下の数を扱うにはどうし
たらよいでしょうか。今の方法で例えば 4 桁の数 a b c d を 
 a*x^3 + b*x^2 + c*x^1 + d  
と表して x に 10 が入った場合が 10 進法でした。では、 x に 0.1 を入れる
とどうなるでしょうか。 d.cba という小数点以下を表す数になるはずです。
大きい整数の計算も小数点以下の計算も繰り上がり処理の方向が違うだけでそ
れまでは全く同じ計算でできそうです。これで、円周率の多倍長計算もできそ
うですね。

FFTW を使った多倍長乗算プログラム 

        real*8 a(4)
        real*8 t,err,tt
	integer*8 n
	integer*8 plan1
	integer*8 plan2
        include "/usr/local/include/fftw3.f"
c 4 要素の FFT を使う
	n=4
c FFT の定義
        call dfftw_plan_r2r_1d(plan1,n,a,a,FFTW_R2HC, FFTW_ESTIMATE)
        call dfftw_plan_r2r_1d(plan2,n,a,a,FFTW_HC2R, FFTW_ESTIMATE)
c  0  0  23 12 の二乗の計算をする
	a(1) = 12
	a(2) = 23
	a(3) = 0
	a(4) = 0
c 実数型 FFT 正変換
        call dfftw_execute(plan1)
c 各要素を二乗する
	a(1)=a(1)*a(1)
	a(n/2+1)=a(n/2+1)*a(n/2+1)
	do i=1,n/2-1
	 t=a(i+1)
	 tt=a(n-i+1)
	 a(n-i+1)=2.0*t*tt
	 a(i+1)=(t+tt)*(t-tt)
	enddo
c 実数型 FFT 逆変換
        call dfftw_execute(plan2)
c 四捨五入で結果を丸める
        do i=1,n
         a(i)=a(i)/n
         t=a(i)
         a(i)=dnint(a(i))
 	 tt=dabs(a(i)-t)
	 if(tt.gt.err)then
	  err=tt
	 endif
        enddo
	write(*,*) a(1),a(2),a(3),a(4)
c FFT 定義の開放
        call dfftw_destroy_plan(plan1)
        call dfftw_destroy_plan(plan2)
        end

fftw の Guru Interface 

	{
	fftw_iodim dims[1];
	dims[0].n = n;
	dims[0].is = 1;
	dims[0].os = 1;
        forw =  fftw_plan_guru_split_dft(1,dims,0,(fftwf_iodim *) 0,x,xw,x,xw,FFTW_ESTIMATE);
        back =  fftw_plan_guru_split_dft(1,dims,0,(fftwf_iodim *) 0,xw,x,xw,x,FFTW_ESTIMATE);
	fftw_execute(forw);
	fftw_execute(back);
	}
で 実部、虚部が別配列の複素数変換。

netlib の dfftpack を使った多倍長乗算プログラム 

        real*8 a(4),work(16)
        real*8 t,err,tt
	integer*8 n
c 4 要素の FFT を使う
	n=4
c FFT の初期化
	call drffti(n,work)
c  0  0  23 12 の二乗の計算をする
	a(1) = 12
	a(2) = 23
	a(3) = 0
	a(4) = 0
c 実数型 FFT 正変換
	call drfftf(n,a,work)
c 各要素を二乗する
	a(1)=a(1)*a(1)
	a(n)=a(n)*a(n)
	 do i=2,n-1,2
	 t=a(i)
	 tt=a(i+1)
	 a(i+1)=2.0*t*tt
	 a(i)=(t+tt)*(t-tt)
	 enddo
c 実数型 FFT 逆変換
	call drfftb(n,a,work)
c 四捨五入で結果を丸める
        do i=1,n
         a(i)=a(i)/n
         t=a(i)
         a(i)=dnint(a(i))
 	 tt=dabs(a(i)-t)
	 if(tt.gt.err)then
	  err=tt
	 endif
        enddo
	write(*,*) a(1),a(2),a(3),a(4)
        end

IBDWT(irrational base discrete wavelet transform)無理数を底とする離散ウェーブレット変換 

http://www.wolfram.com/news/mersenne.ja.html
>IBDWTの考え方は巨大数を無理数を底としたものにまで拡げることです．
>The idea of the IBDWT is to expand giant numbers in an irrational base.

これを発見しメルセンヌ素数界隈に投げ込んだことにより、Richard Crandall 氏
の名前は、GIMPS の新素数発見のプレスリリースに永遠に記載されることになった。
(Crandall氏が Mathematica でいろいろ試しながらアルゴリズムを確立して
それを Cソースに起こし(luctwd.c)、それを GIMPS の GW がアセンブラでチューニング)

IBDWTの具体的なアルゴリズムは
「Prime numbers: a computational perspective By Richard E. Crandall, Carl Pomerance」
に記載されています。
(Σ...とかなので自分には理解不能ですがこの本はクヌースの本を面白いと思う人にはお薦めです)

参考文献

　このページの内容を書くのに使った参考文献です。 

          数学セミナー増刊「計算の効率化とその限界」
          別冊サイエンス「コンピュータ数学」
          和田秀男「コンピュータと素因子分解」「数のはなし」
          クヌース「準数値算法 4 算術演算」
	  一松信「暗号の数理」「数のエッセイ」
	  共立出版「素数大百科」

大型コンピュータによって発見されたメルセンヌ素数

Pyeardiscoverermachine(時間は一つの素数候補をテストする時間)
5211952RAPHAEL M. ROBINSONSWAC(1m)(Standards Western Automatic Computer) NBS(National Bureau of Standards)で作られた。1951年稼働開始。
6071952RAPHAEL M. ROBINSONSWAC(1.5m)
12791952RAPHAEL M. ROBINSONSWAC(13.5m)
22031952RAPHAEL M. ROBINSONSWAC(59m)
22811952RAPHAEL M. ROBINSONSWAC(1h6m)
32171957HANS RIESEL(Riesel prime で有名)BESK(Binary Electronic Sequence Calculator)(5h30m)
42531961ALEXANDER HURWITZIBM-7090 360 アーキティクチャ発表以前で計算機は事務用と技術計算用に別れていた。
44231961ALEXANDER HURWITZIBM-7090(50m)
96891963DONALD B. GILLIESILLIAC-II(1h23m)
99411963DONALD B. GILLIESILLIAC-II(1h30m)
112131963DONALD B. GILLIESILLIAC-II(2h15m)
199371971BRYANT TUCKERMANIBM 360/91 汎用機が世界最速の計算機も兼ねていた時代。(35.01m)
217011978Landon Curt Noll,Laura Nickel(now Ariel Glenn) ２人は高校生であった。CDC-CYBER-174(7h40m20s)
232091979Landon Curt NollCDC-CYBER-174 当時、世界で３番目に高速な計算機であった。(8h39m37s)
444971979HARRY L. NELSON,DAVID SLOWINSKICRAY-1(マシンサイクルが 13.5 nsec 1マシンサイクル毎に乗算と加算の結果を得られた)
862431982DAVID SLOWINSKICRAY 1S(1h36m22s)
1105031988Walter N. Colquitt & Luther Welsh,Jr.SX-2
1320491983DAVID SLOWINSKICRAY X-MP(30m)
2160911985DAVID SLOWINSKICRAY X-MP/24
7568391992DAVID SLOWINSKI,Paul GageCRAY 2
8594331994DAVID SLOWINSKI,Paul GageCRAY C916
12577871996DAVID SLOWINSKI,Paul GageCRAY T94

メルセンヌ素数探索プロジェクトについて

　メルセンヌ素数とは 2^N-1(N は自然数:素数でないとメルセンヌ素数に ならないことが知られている)の形をした素数です。メルセンヌ牧師という 人がこの素数に強い興味を持ったことからこの名がつきました。現代にお いてこの数が注目されているのはフェルマの小定理に似た、ルーカスチェ ックという素数判定法が2進法を使った計算機の上ではメルセンヌ素数 に対して有効に働くためです。 理論上はルーカスレーマチェックと同じ位に効率的なある形の素数に 対して有効な素数判定法があり、メルセンヌ素数の桁が毎回大きくなりすぎる ことなどから、計算機時代の初期には最大の素数がメルセンヌ数でないこともありましたが、 下記のプロジェクトなどによってしばらくはメルセンヌ数の天下が続くと思われます。 (実は Generalized Fermat prime の発見プログラムの方がはるかに単純)

あなたもオイラーと並んで数学史上に永遠に名をとどめてみませんか


メルセンヌ素数というのは計算機の速度が進歩すれば忘れ去られてしまう、 「RC64」とか「コンピュータがチェスで人間に勝った」とかと違って、 整数という人類が不変の興味を持って来た一分野でありあなたがそれを 発見すれば貴方の名前はしばらくは忘れ去られることはないでしょう。

パソコンメルセンヌ素数探索のメッカ http://www.mersenne.org/

こちらの Web page では、ペンティアムの互換機を 御持ちで数週間の夜間JOB(昼間は止めておけるプログ ラムになっている) を実行して頂ける協力者の方を求 めています。 (ただし根気良くマシンを回し続けられるかなり暇な方)

あなたのコンピュータディスプレイに稲妻のように 世界最大の素数が輝く瞬間を体験してみませんか。

このなかにある ML にはその筋では有名な高 校生で当時世界最大の素数を発見したクルトノルや Cray のストロビィンスキーなども参加しています。

GIMPSに参加している計算機の規模 

2004年1月現在 GIMPS での mersenne 素数探索に参加しているコンピュータは

http://mersenne.org/primenet/

によると

              ------- Aggregate CPU Statistics, P90 Units* -------
 
                  Last 7 Days Average             Cumulative Today
                 from 11 Jan 2004 06h           from 17 Jan 2004 06h
                ----------------------  ----------------------------------
  Test Type     CPU yr/day    GFLOP/s    CPU years  CPU yr/day    GFLOP/s
  ------------  ----------  ----------  ----------  ----------  ----------
  Lucas-Lehmer   1001.883   12060.372      36.275     886.607   10672.712
  Factoring        30.873     371.646       1.337      32.668     393.245
                ----------  ----------  ----------  ----------  ----------
  TOTALS         1032.757   12432.018      37.612     919.275   11065.957

で、P90 換算で 30万台程度のようだ。(11Tflops!)
この能力は、次のメルセンヌ素数を探すにはそれまでのメルセンヌ素数を全て
検査する時間に等しい時間がかるといわれる、
2^20996011-1 までのルーカスチェックを１年で終えることができるようだ。

# P4 3Ghz 換算だと 3000台になります。

(*Measured in calibrated P5 90Mhz, 32.98 MFLOP units: 25658999 FPO / 0.778s using 256k FFT.)
                 
という記述があり、P4 3Ghz マシンは P90 より 100倍程度速いのでなんと 3Gflops
出ていることになる。(FFT のサイズが大きくなるとその差は50倍程度につまる)


それに対して distributed.net の RC5-72 は

Current Information

2,951,343 Blocks were completed yesterday (0.000268% of the keyspace)
at a sustained rate of 34 Blocks/sec.

12,675,921,664,278,527 Keys were completed yesterday (0.000268% of the keyspace)
at a sustained rate of 146,712,056,300 Keys/sec.

ということで、P90 換算で 150 万台程度。

意外に GIMPS が健闘していることがわかります。


RC5 系の整数型計算と mersenne 素数の浮動小数点
計算での、P90 と P4 マシンの性能差はどうでしょうか。

http://www.orange.co.jp/~masaki/rc5/ratej.html

の RC64 のベンチマークによると

132.3  Pentium-90 の 132.3 に対して Pentium 4-2260 の  4083.3  
と性能差は 30倍位です。

mersenne 素数の計算もだいたい同じ位です。
浮動小数点の性能向上の方が著しいと思ったのですが
そうでもないですね。
prime95 の Mflops 値 

http://www.fftw.org/speed/ 

によると fft の 要素数(N) と計算時間と Mflops 値 の関係は

mflops = 5 N log2(N) / (time for one FFT in microseconds)

/ 2 for real-data FFTs 

Mflops: Million floating-point operation per second (毎秒100万回の浮動小数点演算)

ここに書かれている演算量は理論値で実際に計算されている浮動小数点演算の量は
もう少し少ないことが多い。

256k 要素だと 1 回 の fft に

  5 * 2^18 * log2(2^18) / 2 = 90 * 2^18 /2  = 12M 位の浮動小数点演算

* mersenne 素数の計算には実数型の fft が使われる。

http://mersenne.org/primenet/

によると

*P90 CPU time according to Woltman/Kurowski formulation. Calibrated by
benchmark P5 90Mhz, 32.98 MFLOP units: 25658999 FLOP/0.778s (256k FFT).

256k の FFT で  25M 位の計算を 0.778s で行うと 32Mflops としている。

http://www.mersenne.org/bench.htm

によると

┏━━━━━━━━━━━━━━━━━━━━━━━┯━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓
┃                     CPU                      │ Iteration Times For Various Exponent Ranges and FFT sizes        ┃
┠───────┬───┬───┬───┬───┼───┬───┬───┬───┬───┬───┬────┬────┨
┃              │      │      │      │      │6.52M │7.76M │9.04M │10.33M│12.83M│15.3M │ 17.85M │ 34.56M ┃
┃              │Speed │Memory│  L2  │  L2  │  to  │  to  │  to  │  to  │  to  │  to  │   to   │   to   ┃
┃     Type     │(MHz) │Speed │Cache │Cache │7.76M │9.04M │10.33M│12.83M│15.3M │17.85M│ 20.4M  │ 39.50M ┃
┃              │      │      │ Size │Speed │(384K)│(448K)│(512K)│(640K)│(768K)│(896K)│(1024K) │(2048K) ┃
┠───────┼───┼───┼───┼───┼───┼───┼───┼───┼───┼───┼────┼────┨
┃   P4         │ 3060 │      │ 512  │ Full │ 0.013│ 0.016│ 0.018│ 0.023│ 0.028│ 0.034│   0.037│   0.068┃
┠───────┼───┼───┼───┼───┼───┼───┼───┼───┼───┼───┼────┼────┨
┃   Pentium    │ 100  │  66  │ 256  │ Bus  │ 0.988│ 1.188│ 1.328│ 1.773│ 2.139│ 2.558│   2.866│   -----┃
┗━━━━━━━┷━━━┷━━━┷━━━┷━━━┷━━━┷━━━┷━━━┷━━━┷━━━┷━━━┷━━━━┷━━━━┛

となっていて、ここで言う Iterration Times というのは FFT 2回分の時間らしい。(乗算の前後で 2 回 FFT を行うため)

 512k での演算量は、

  5 * 2^19 * log2(2^19) /2 *2 = 95 * 2^19 = 50M 回の浮動小数点演算

 P100 の性能は、

 50M / 1.328s = 37.7Mflops

 P4 3Ghz の性能は、

 50M / 0.018s = 2778Mflops

 2048k での演算量は、

  5 * 2^21 * log2(2^21) /2 *2 = 105 * 2^21 = 210M 回の浮動小数点演算

 P4 3Ghz の性能は、

 210M / 0.068s = 3088Mflops

となる。

The time it takes to test one exponent is the exponent being tested multiplied by the iteration time. For example, 
you own a 2.4 GHz Pentium 4 machine and are considering running first-time tests, if the PrimeNet server shows that 
the smallest available first-time exponents are around 24,000,000, find the column marked 20.18M to 25.09M and the 
row containing your CPU to find a per-iteration speed of 0.064 seconds. The time it takes to run one first-time 
test is 24,000,000 * 0.064 seconds, or 1,536,000 seconds, or (divide by 86,400) 17.8 days.

2.4Ghz P4 で P = 24,000,000 だと 20.18M から 25.09M のレンジで 一回に 0.064 sec かかるので

24,000,000 * 0.064sec で約 17.8 days で計算が終わる と書かれている。 

Mlucas を実行した結果

P=20996011 の時 1280k の fft を使い 16bit/要素 1280k*16bit が 20996011 bit の倍(乗算すると桁が倍になるので)なので話があう。


ちなみに、P4 は SSE2 命令を使用すれば 1clock毎に倍精度浮動小数点演算を 2 個 発行できるため、4Ghz の clock の
マシンで ピーク性能は 8Gflops である。
プログラムを動かす 

EB164(21164) で メルセンヌ素数探索プログラムを動かしてみる。

http://hogranch.com/mayer/README.html
http://hogranch.com/mayer/bin/ALPHA_UNIX/Mlucas-2.7b-ev4-4x.tar.gz

から Tru64Unix 用のバイナリをダウンロードする。
ソースも置いてあるが f90 なのでコンパイルする環境がない。

#Linux/alpha 用のバイナリは、NetBSD を簡単に crash させる。

alpha ~/Mlucas$./Mlucas-2.7b-ev4-4x 
  looking for worktodo.ini file...
  no worktodo.ini file found...switching to interactive mode.
 Enter exponent, FFT length in K (set K = 0 for default FFT length) >216091  ← 調べたいメルセンヌ素数
16                                                                           ← FFT の要素数(エラーが出ない範囲で小さい方が速い。大き過ぎてもエラーになるようだ)
 Enter 'y' to run a self-test,  for a full LL test >
 Enter index of radix set to be used for the FFT: 
(See file fft_radix.txt for a list of available choices; enter -1 to get the default)  >0  ← 2 の倍数でない 2 の倍数でない FFT にどんな組合せを使うかの指定(0でいいようだ)
 Enter 'y' to enable per-iteration error checking,  for no error checking >y       ← y にしないとエラーチェックしないので結果が間違ってるかもしれない
  p is prime...proceeding with Lucas-Lehmer test...
M(   216091 ): using FFT length       16K =    16384 8-byte floats.
  this gives an average    13.1891479492188      bits per digit
 INFO: Using real*16 for FFT sincos and DWT weights tables inits.
 Using complex FFT radices           8           8           8          16
M(   216091 ) iteration =     2000 clocks = 00:00:17.678. Res64: 4B3F667335A12495  
M(   216091 ) iteration =     4000 clocks = 00:00:06.629. Res64: B265971DA23E272F  
....
M(   216091 ) iteration =   214000 clocks = 00:00:07.512. Res64: B5E1A392091C916B  
M(   216091 ) iteration =   216000 clocks = 00:00:07.435. Res64: F5AB54BCEDDF95AD  
M(  216091 ) is a known MERSENNE PRIME.

非常に高速である。

現在の記録の 2^20996011-1 に対して
He used a 2 GHz Pentium 4 Dell Dimension PC running for 19 days to prove the number prime.
だそうである。19日で１個調べられるなら流す気がおきますね。

XP1000(500Mhz 21264)で
bash-2.05b$ ./Mlucas-2.7b-gen-4x 
  looking for worktodo.ini file...
  no worktodo.ini file found...switching to interactive mode.
 Enter exponent, FFT length in K (set K = 0 for default FFT length) >20996011
1280
 Enter 'y' to run a self-test,  for a full LL test >
 Enter index of radix set to be used for the FFT: 
(See file fft_radix.txt for a list of available choices; enter -1 to get the default)  >0
 Enter 'y' to enable per-iteration error checking,  for no error checking >y
  p is prime...proceeding with Lucas-Lehmer test...
M( 20996011 ): using FFT length     1280K =  1310720 8-byte floats.
  this gives an average    16.0186851501465      bits per digit
 INFO: Using real*16 for FFT sincos and DWT weights tables inits.
 Using complex FFT radices           5           8           8           8
          16          16
M( 20996011 ) iteration =     2000 clocks = 00:17:06.080. Res64: D36F95D9D040977B

で 123日位かかる計算なので Pen 4 との比較ではだいぶ遅そう。


http://www.mersenne.org/freeware.htm
にある mers.tar.gz
(今は http://www.mersenne.org/various/freeware.htm)

EB164(300Mhz 21164)で

alpha ~/Mlucas$time ./fftlucas 216091
Dec 30 03:40:44 alpha /netbsd: osf1_getsysinfo called with unknown op=67
Dec 30 05:33:07 alpha ntpdate[983]: step time server 133.100.9.2 offset 5.872221 sec
M( 216091 )P, n = 16384, fftlucas v6.50  Crandall

real    302m25.074s
user    301m38.702s
sys     0m0.453s

EB164(300Mhz 21164)で

alpha ~/Mlucas$time ./mersenne1 216091
M( 216091 )P, n = 16384, mersenne1 v6.50  Kline

real    121m38.889s
user    121m22.696s
sys     0m0.933s

EB164(300Mhz 21164)で

alpha ~/Mlucas$time ./mersenne2 216091
M( 216091 )P, n = 16384, mersenne2 v6.50  Kline

real    135m18.800s
user    135m1.832s
sys     0m0.260s

EB164(300Mhz 21164)で

alpha ~/Mlucas$time ./Mlucas-2.7b-ev4-4x 
  looking for worktodo.ini file...
  no worktodo.ini file found...switching to interactive mode.
 Enter exponent, FFT length in K (set K = 0 for default FFT length) >216091
12
 Enter 'y' to run a self-test,  for a full LL test >
 Enter index of radix set to be used for the FFT: 
(See file fft_radix.txt for a list of available choices; enter -1 to get the default)  >0
 Enter 'y' to enable per-iteration error checking,  for no error checking >y
  p is prime...proceeding with Lucas-Lehmer test...
M(   216091 ): using FFT length       12K =    12288 8-byte floats.
  this gives an average    17.5855305989583      bits per digit
 INFO: Using real*16 for FFT sincos and DWT weights tables inits.
 Using complex FFT radices           6           8           8          16
M(   216091 ) iteration =     2000 clocks = 00:00:16.583. Res64: 4B3F667335A12495
...
M(   216091 ) iteration =   214000 clocks = 00:00:10.929. Res64: B5E1A392091C916B
M(   216091 ) iteration =   216000 clocks = 00:00:10.941. Res64: F5AB54BCEDDF95AD
M(  216091 ) is a known MERSENNE PRIME.

real    19m47.505s
user    19m35.922s
sys     0m0.946s

XP1000(500Mhz 21264)で
bash-2.05b$ time ./fftlucas 216091
M( 216091 )P, n = 16384, fftlucas v6.50  Crandall

real    180m36.110s
user    180m23.176s
sys     0m0.191s

XP1000(500Mhz 21264)で
bash-2.05b$ time ./mersenne1 216091
M( 216091 )P, n = 16384, mersenne1 v6.50  Kline

real    82m34.257s
user    82m28.432s
sys     0m0.074s

XP1000(500Mhz 21264)で
bash-2.05b$ time ./Mlucas-2.7b-gen-4x
  looking for worktodo.ini file...
  no worktodo.ini file found...switching to interactive mode.
 Enter exponent, FFT length in K (set K = 0 for default FFT length) >216091
12
 Enter 'y' to run a self-test,  for a full LL test >
 Enter index of radix set to be used for the FFT: 
(See file fft_radix.txt for a list of available choices; enter -1 to get the default)  >0
 Enter 'y' to enable per-iteration error checking,  for no error checking >y
  p is prime...proceeding with Lucas-Lehmer test...
M(   216091 ): using FFT length       12K =    12288 8-byte floats.
  this gives an average    17.5855305989583      bits per digit
 INFO: Using real*16 for FFT sincos and DWT weights tables inits.
 Using complex FFT radices           6           8           8          16
M(   216091 ) iteration =     2000 clocks = 00:00:14.575. Res64: 4B3F667335A12495
...
M(   216091 ) iteration =   216000 clocks = 00:00:04.617. Res64: F5AB54BCEDDF95AD
M(  216091 ) is a known MERSENNE PRIME.

real    8m29.231s
user    8m16.145s
sys     0m0.346s

手持ち最速の Pen マシンでは、
CPU: Celeron (730.90-MHz 686-class CPU)
OSは、FreeBSD
(このマシンの計測データは他のタスクも走っている。メインマシンなので)

まずは、最速の探索プログラムも動くのでその benchmark を流す。

Your timings will be written to the results.txt file.
Compare your results to other computers at http://www.mersenne.org/bench.htm
That web page also contains instructions on how your results can be included.

Timing 36 iterations at 256K FFT length.  Best time: 72.777 ms.
Timing 29 iterations at 320K FFT length.  Best time: 95.155 ms.
Timing 24 iterations at 384K FFT length.  Best time: 114.680 ms.
Timing 20 iterations at 448K FFT length.  Best time: 136.719 ms.
Timing 18 iterations at 512K FFT length.  Best time: 153.295 ms.
Timing 14 iterations at 640K FFT length.  Best time: 198.322 ms.
Timing 12 iterations at 768K FFT length.  Best time: 241.780 ms.
Timing 10 iterations at 896K FFT length.  Best time: 282.542 ms.
Timing 10 iterations at 1024K FFT length.  Best time: 325.200 ms.
Timing 10 iterations at 1280K FFT length.  Best time: 416.351 ms.
Timing 10 iterations at 1536K FFT length.  Best time: 505.853 ms.
Timing 10 iterations at 1792K FFT length.  Best time: 603.524 ms.

http://www.mersenne.org/bench.htm

┏━━━━━━━━━━━━━━━━━━━━━━━┯━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓
┃                     CPU                      │                Iteration Times For Various Exponent Ranges and FFT sizes                 ┃
┠───────┬───┬───┬───┬───┼───┬───┬───┬───┬───┬───┬────┬────┬────┬────┬┬┨
┃              │      │      │      │      │6.52M │7.76M │9.04M │10.33M│12.83M│15.3M │ 17.85M │ 20.4M  │ 25.35M │ 30.15M ││┃
┃              │Speed │Memory│  L2  │  L2  │  to  │  to  │  to  │  to  │  to  │  to  │   to   │   to   │   to   │   to   ├┼┨
┃     Type     │(MHz) │Speed │Cache │Cache │7.76M │9.04M │10.33M│12.83M│15.3M │17.85M│ 20.4M  │ 25.35M │ 30.15M │ 35.1M  ││┃
┃              │      │      │ Size │Speed │(384K)│(448K)│(512K)│(640K)│(768K)│(896K)│(1024K) │(1280K) │(1536K) │(1792K) ├┼┨
┃              │      │      │      │      │      │      │      │      │      │      │        │        │        │        ││┃
┠───────┼───┼───┼───┼───┼───┼───┼───┼───┼───┼───┼────┼────┼────┼────┼┼┨
┃   P4         │ 3060 │      │ 512  │ Full │ 0.013│ 0.016│ 0.018│ 0.023│ 0.028│ 0.034│   0.037│   0.053│   0.065│   0.081││┃
┠───────┼───┼───┼───┼───┼───┼───┼───┼───┼───┼───┼────┼────┼────┼────┼┼┨
┃   Celeron    │ 733  │  66  │ 128  │ Full │ 0.138│ 0.165│ 0.189│ 0.246│ 0.301│ 0.358│   0.407│   0.524│   0.647│   0.775││┃
┠───────┼───┼───┼───┼───┼───┼───┼───┼───┼───┼───┼────┼────┼────┼────┼┼┨
┃   Pentium    │ 100  │  66  │ 256  │ Bus  │ 0.988│ 1.188│ 1.328│ 1.773│ 2.139│ 2.558│   2.866│   3.713│   4.459│   5.354││┃
┗━━━━━━━┷━━━┷━━━┷━━━┷━━━┷━━━┷━━━┷━━━┷━━━┷━━━┷━━━┷━━━━┷━━━━┷━━━━┷━━━━┷┷┛

とあり、大体あっている。現在最速の P4 とは10倍の開きがあるようだ。
(1万円で買ったマシンなので価格程の開きはないようだ)

Celeron (730.90-MHz 686-class CPU) で

テスト中に任意の P についてテストできるモードがあるので試す。

        20.  Help/About PrimeNet Server
Your choice: 7

Exponent to test: 216091

Accept the answers above? (Y): y
Starting primality test of M216091
[Dec 30 17:45] Iteration: 10000 / 216091 [4.62%].  Per iteration time: 0.002 sec.
[Dec 30 17:45] Iteration: 20000 / 216091 [9.25%].  Per iteration time: 0.002 sec.
[Dec 30 17:45] Iteration: 30000 / 216091 [13.88%].  Per iteration time: 0.002 sec.
...
[Dec 30 17:51] Iteration: 200000 / 216091 [92.55%].  Per iteration time: 0.002 sec.
[Dec 30 17:51] Iteration: 210000 / 216091 [97.18%].  Per iteration time: 0.002 sec.
M216091 is prime! WY6: 8368E9FF

7分程のようだ。

Celeron (730.90-MHz 686-class CPU) で

compaq ~/p$time ./fftlucas 216091
M( 216091 )P, n = 16384, fftlucas v6.50  Crandall

real    135m34.955s
user    130m46.688s
sys     0m12.536s

Celeron (730.90-MHz 686-class CPU) で

compaq ~/p$time ./mersenne1 216091
M( 216091 )P, n = 16384, mersenne1 v6.50  Kline

real    121m15.626s
user    93m15.641s
sys     0m15.756s

まとめると

2^216091-1 に対するベンチマーク
(感触として 2^216091-1 に対する性能と現在探査が行われている範囲との性能は比例しそうである。)
		   XP1000(500Mhz 21264)	EB164(300Mhz 21164) 730Mhz Celeron       2Ghz P4
本番プログラム	              -			-   	          7m          1m以下と予想される
Mlucas-2.7b-ev4-4x 	      8m 	       20m		  -
fftlucas	            180m 	      302m		135m
mersenne1	             82m 	      121m    	        121m
現在の中古価格            5万円               1万円             1万円              10万円            

Mlucas-2.7b-ev4-4x は、Pen 用の本番プログラムと比較してチューニングの度合
としてはがんばっていそうである。alpha は、Pen と比較してこの問題に
向いていないという結論が出てしまいそう。

http://www.ocrat.com/mersenne/
にある mayerE25C3.tar.gz  を動かす。

xp1000 $./lucas_mayer 
 no restart file found...looking for range file...
 no range file found...switching to interactive mode.
 Enter p,n (set n=0 for default FFT length) >11213,0
 Enter "y" to run a self-test,  for a full LL test >
 p is prime...proceeding with Lucas-Lehmer test...
 using an FFT length of 512
 this gives an average 21.900391 bits per digit
calculating sin/cos values...done
 M(    11213 ) iteration =     2000 clocks = 00:00:00.000
 M(    11213 ) iteration =     4000 clocks = 00:00:01.000
 M(    11213 ) iteration =     6000 clocks = 00:00:00.000
 M(    11213 ) iteration =     8000 clocks = 00:00:01.000
 M(    11213 ) iteration =    10000 clocks = 00:00:01.000
 M(   11213 ) is a known MERSENNE PRIME.

これもなかなか高速である。

GIMPS Source Code Timings Page
http://hogranch.com/mayer/gimps_timings.html

を見ると Glucas というプログラムがあるようだ。

GLUCAS - Yet Another FFT
http://glucas.sourceforge.net/

からソースをダウンロードして、./configure make する。

glucas.que というファイルを作って中に 調べたい P を入れる。

216091 で試してみると

xp1000 glucas-2.9.0$time ./Glucas 

real    9m56.624s
user    9m52.195s
sys     0m0.079s

glucas.out に結果が入る。

Iter. 200000 ( 92.55%), Err= 0.000, 26.87 user  99% CPU (0.003 sec/iter).
M216091. Saved file at iteration 200704. Res64: F4E92A3E32F77E9F.
M216091. Saved file at iteration 204800. Res64: 78897F6C7A3E982E.
M216091. Saved file at iteration 208896. Res64: 1007948E422A60B9.
Iter. 210000 ( 97.18%), Err= 0.000, 26.88 user  99% CPU (0.003 sec/iter).
M216091. Saved file at iteration 212992. Res64: 240FD3A839CB978D.
[Wed Feb  4 05:08:51 2004]
M216091 is a known Mersenne prime!
Terminated all the queued job in file: glucas.que.

P=20996011 で動かしてみる。

[Wed Feb  4 05:26:54 2004]
Going to work with exponent 20996011
Starting from iteration 1. Exponent 20996011.
M20996011. Saved file at iteration 4096. Res64: BFA31B6520E49912.
M20996011. Saved file at iteration 8192. Res64: 141E6DE86BB2DD30.
Iter.    10000 (  0.05%), Err= 0.148, 5128.38 user 100% CPU (0.513 sec/iter).
M20996011. Saved file at iteration 12288. Res64: 59F8E563958FF7CC.
M20996011. Saved file at iteration 16384. Res64: 239A8FBAC9CE49E6.
Iter.    20000 (  0.10%), Err= 0.141, 5128.05 user 100% CPU (0.513 sec/iter).

alpha で動作する ソース付きとしてはこれが一番高速である。


http://mersenne.org/primenet/status.txt

にステータスのページがある
 
 prime      fact  current          days
exponent    bits iteration  run / to go / exp   date updated     date assigned   a
-------- -- ---- ---------  -----------------  ---------------  ---------------  -
21039677     67   8816895    95.0  65.4  77.4  24-Jan-04 15:54  01-Nov-03 05:02  S

が私のパソコンで実行されている。
67bit まで試行除算されている。
試行除算プログラムを動かす 

これも mers.tar.gz に含まれる試行除算のプログラム
(今は http://www.mersenne.org/various/freeware.htm)

mersfacgmp

/pkgsrc/devel/gmp がインストールされている必要がある。

xp1000 $./mersfacgmp 
12000257
M( 12000257 )C: 96002057
M( 12000257 )J: 24000515 96002057
M( 12000257 )H: 96002057
216091
M( 216091 )J: 432183 42485219329
M( 216091 )U: 42485219329

のように実行する。割り切れた場合は、C: が表示される。
J: が除算を試行した値の最小値と最大値

xp1000 $./mersfacgmp  -e42
216091
M( 216091 )J: 432183 4418462810113
M( 216091 )U: 4418462810113

試行する桁を増やす場合は、-e で指定する。(log みたい)

#何故か FreeBSD/i386 では上手く動かなかった。

mersfaclip

コンパイルには、
http://www.und.edu/org/crypto/crypto/numbers/programs/freelip/freelip_1.1.tar.gz
が必要。

xp1000 $time ./mersfaclip -e45 < 216091 
M( 216091 )J: 432183 35220246822913
M( 216091 )U: 35220246822913

real    1m23.136s
user    1m23.032s
sys     0m0.005s

compaq $time ./mersfaclip -e45 < 216091 
M( 216091 )J: 432183 35220246822913
M( 216091 ): 35220246822913

real    1m37.294s
user    1m36.804s
sys     0m0.032s

このプログラムは FreeBSD/i386 でも動いた。


mersfacgmp を使って

2^20996011-1 〜 2^21012671-1 の 1000個のメルセンヌ数について素因数を調べてみる。

P1000 に p を格納し

xp1000 $./mersfacgmp -e 55 < P1000 > D1000

この範囲で素因数が見付かったのが 539個

素因数の(10進での)桁数毎の個数を見てみると

 8桁    33
 9桁    98
10桁    91
11桁    74
12桁    67
13桁    82
14桁    38
15桁    38
16桁    37
17桁     9

となった。

http://hogranch.com/mayer/README.html

にある  Mfactor.c.gz も試行除算プログラム

cc -O4 Mfacytor.c でコンパイルして P=21008633 に対して実行してみると

xp1000 Mlucas$time ./a.out 21008633 0 1271201920312151
# Mersenne factorer.  Author:  Peter L. Montgomery, pmontgom@cwi.nl
# $Id: Mfactor.c,v 1.13 1997/05/19 05:34:22 wedgingt Exp $
# Compiled by GCC on unknown architecture under unknown OS.
# ARITHMETIC_MODEL = ARITHMETIC_FP64
# NUM_SIEVING_PRIME = 34567,   SIEVE_LENGTH_BITS = 589824
# SIEVE_STEP_BIT = 840
# TRYQ_BLOCKING = 3,     USING_ASSEMBLY = 0,    USING_POPULATION_COUNT = 0
# VECTOR = 0,    SHORT_VECTOR_LENGTH = 20
# Type sieving_prime_t has 32 bits.
# Type wsieve_t has 64 bits, of which all 64 are used.
# Type natural_t is signed and has 64 bits.
# Floating point mantissa has 53 significant bits.
# QBITSMAX = 53
#
# Invoked as ./a.out 21008633 0 1271201920312151
# Running on xp1000 at Wed Feb  4 12:52:18 2004
# No history file specified (use -h parameter).
# Results will be written to standard output only.
#
# p = 21008633.  Searching q in [1, 1271201920312151]
# The isieve loop in tryp sieves an interval of length
#    10408772598497280 = 1.041e+16  in NSIEVE = 96 pieces.
# It will need    0.12 passes to search from 1 to 1271201920312151.
# Your SIEVE_STEP_BIT (or SIEVE_LENGTH_BITS) seems high for your search limits.
# 34567-th sieving prime is 409639.  25593 are below SIEVE_LENGTH_BITS/2,
# 8974 more below SIEVE_LENGTH_BITS, 0 above.
M( 21008633 )C: 1271201920312151  # xp1000  @  Wed Feb  4 12:52:22 2004
# 1323070 values of q tried out of 6915264 sieved.
#  80.87% of potential q's (mod 840*p) were eliminated by sieving.
M( 21008633 )U: 1271201920312151  # xp1000  @  Wed Feb  4 12:52:22 2004

real    0m4.161s
user    0m4.150s
sys     0m0.008s

となる。

mersfacgmp で同じ P で実行すると

xp1000 mers$time ./mersfacgmp -e50 < 21008633 
M( 21008633 )J: 42017267 1127617031503873
M( 21008633 )U: 1127617031503873

real    0m19.736s
user    0m19.725s
sys     0m0.007s

となる Mfactor の方が数倍高速である。
lucdwt.c 

http://www.mersenne.org/gimps/lucdwt.c
(今は http://mersenneforum.org/gimps/lucdwt.c)

GIMPS にある一番単純なルーカスチェックのプログラム
Richard E. Crandall氏がこれを GIMPS に提供し、
Prime95 はこれをチューングしたようです。

$ ./a.out 11213 512 0 1
.........
11209 maxerr: 0.078125
11210 maxerr: 0.085938
11211 maxerr: 0.000000
11213  0

0 になればメルセンヌ素数。
(引数は、テストするp、FFT要素数、途中で停止する場合のくり返し数、くり返し毎に print するか)

lucdwt.c:
    880         for(j=1;j< last;j++)
    881         {
    882                 err = lucas_square(x,n,errflag);
    883                 if (errflag)
    884                 {
    885                         printf("%ld maxerr: %f\n",j,err);
    886                         fflush(stdout);
    887                 }
    888                 x[0] -= 2.0;
    889         }
メインループ二乗するたび -2 してます。

    415 void
    416 squareg(
    417         double  *x,
    418         int     size
    419 )
    420 {
    421         fft_real_to_hermitian(x, size);
    422         square_hermitian(x, size);
    423         fftinv_hermitian_to_real(x, size);
    424 }

FFT して各要素を二乗して逆FFT
Generalized Fermat Prime 

mersenne 素数以外に同じような効率で素数検査できる数に Generalized Fermat Prime
がある。

b^N+1 ( N は 2^n ) という形の数。

http://primes.utm.edu/primes/lists/all.txt

にある

 順位                                   桁数        発見年   性質
  498  1287484^16384+1                  100103 g197 2001 Generalized Fermat
 1198  2^216091-1                        65050 S    1985 Mersenne 31
 1199  87077894^8192+1                   65044 g303 2003 Generalized Fermat
 4309  1430330^8192+1                    50426 g162 2001 Generalized Fermat
 4744  764898^8192+1                     48199 g232 2001 Generalized Fermat
 4995  504216^8192+1                     46716 g6   2000 Generalized Fermat

をテストに使う。

http://perso.wanadoo.fr/yves.gallot/primes/download.html

には、ソース付きのプログラムがある。

genefer80.gz: binary for Linux of the 80-bit version of the program.

は、バイナリのみだが x86 の 80bit版の命令を使って 探査する桁数によっては
高速だという。

Linux 上 Celern 533Mhz で
$ time ./genefer80 
GeneFer 1.3 (80-bit x86-ASM)  Copyright (C) 2001-2003, Yves GALLOT
A program for finding large probable generalized Fermat primes.

Usage: GeneFer -b            run bench
       GeneFer -t            run test
       GeneFer <filename>    test <filename>
       GeneFer               use interactive mode

  1. bench
  2. test
  3. normal

N: 8192
b: 504216
504216^8192+1 is a probable prime. (err = 1.53e-05)

real    16m35.914s
user    16m17.690s
sys     0m0.010s


genefer64.gz : binary for Linux.

は通常の 倍精度を使った版。

Linux 上 Celern 533Mhz で
[msft@funaki tmp]$ time ./genefer64 
GeneFer 1.3 (x86-ASM)  Copyright (C) 2001-2003, Yves GALLOT
A program for finding large probable generalized Fermat primes.

Usage: GeneFer -b            run bench
       GeneFer -t            run test
       GeneFer <filename>    test <filename>
       GeneFer               use interactive mode

  1. bench
  2. test
  3. normal

N: 8192
b: 504216
504216^8192+1 is a probable prime. (err = 9.59e-03)

real    8m57.813s
user    8m49.060s
sys     0m0.000s

こちらの方が速い。まだサイズが小さいのかもしれない。

ソース版の genefer.c を使ってみる。

Linux 上 Celern 533Mhz で
$ cc -O3 ./genefer.c -lm
$ time ./a.out 
GeneFer 1.3 (standard)  Copyright (C) 2001-2002, Yves GALLOT
A program for finding large probable generalized Fermat primes.

Usage: GeneFer -b            run bench
       GeneFer -t            run test
       GeneFer <filename>    test <filename>
       GeneFer               use interactive mode

  1. bench
  2. test
  3. normal

N: 8192
b: 504216

504216^8192+1 is a probable prime. (err = 2.16e-02)

real    18m43.390s
user    18m27.520s
sys     0m0.040s

alpha 21264 500Mhz/XP1000/NetBSD/gcc での結果

N: 8192
b: 504216
504216^8192+1 is a probable prime. (err = 2.73e-02)

real    4m28.451s
user    4m19.482s
sys     0m0.006s

なかなか良い結果。

linux/Celeron 533Mhz での測定結果(user time)

	N		|	genefer64	|	genefer80	|	mprime(GMIPS)
------------------------+-----------------------+-----------------------+---------------------
504216^8192+1	    	|	8m49.060s	|	16m17.690s	|
764898^8192+1	    	|	10m9.740s	|			|
2^216091-1		|	---		|	---		|	12m28.510s
87977894^8192+1		|	12m15.350s	|			|
1287484^16384+1     	|	32m52.700s 	|			|

87977894^8192+1 と 2^216091-1 はほぼ同じ桁数だが
計算時間も ほぼ同じでこのプログラムが良く最適化されていることがわかる。
(本当は P4 で比較しないと意味ないけど)

Athlon Optimized GFN Deep Sieve

には、Windows上で動作する GFN 用のふるいプログラムがある。
通常の Windows 版は、US Yahoo のフォーラムのみで配布されているようだ。

http://www.utm.edu/research/primes/programs/gallot/

には Proth.exe  がある。

 Proth's Theorem (1878): Let N = k^(2^n)+1 with 2^n > k.   

If there is an integer a such that

a^(N-1)/2 = -1 (mod N), 

then N is prime. 

N = k^(2^n)+1 で 2**n > k の時 ヤコビ記号[a,N] = -1 
a^(N-1)/2 = -1 (mod N) なる a が存在すれば N は素数

という Proth's が証明した定理に基づいて素数を探索するプログラム。
(実際には、このタイプ以外の素数も探索できるように拡張されている)

http://netnews.gotdns.org/WallStreet/6351/txt/gfn.f

自分でも作ってみた。

-- P4 でも測定した

Freebsd/P4 1.4GMhz での測定結果(user time)

	N		|	genefer64	|	genefer80	|	mprime(GMIPS)
------------------------+-----------------------+-----------------------+---------------------
504216^8192+1	    	|	1m31.258s	|	4m54.811s 	|
764898^8192+1	    	|	1m34.290s	|			|
2^216091-1		|	---		|	---		|	 1m58.912s
87977894^8192+1		|	2m7.615s  	|			|
2^756839-1		|	---		|	---		|	31m28.277s
10021136^32768+1	|      37m56.906s(ERR) 	|			|
linux 上で TOP5000 に入る素数を見つける方法 

mersenne 素数を除くと k.2^n-1 の形の素数が現在は、一番高速に探索できるようだ。

http://www.utm.edu/research/primes/programs/NewPGen/

より linux 用の newpgen をダウンロードする。

$ ./snewpgen 

NewPGen 2.82 - Use Control-C to manually stop the current sieving.

What do you want to do now ?
        0 : To create a new file.
        1 : To continue sieving an old file.
        2 : To change one or several option(s).
        3 : To get a service.
        4 : To get some help.
        5 : To exit from this program.
Please, enter your choice (0) : 0
Verify results ? (n) : 
Creating a new file ; Output file name = Newpgen  ← 出力ファイル
Using primorial mode ? (n) : 
What type of sieve do you want to do ?
The available types of sieve are :
        0 : k.b^n+1
        1 : k.b^n-1
        2 : Twin: k.b^n
        3 : SG: k.b^n-1 & 2k.b^n-1
        4 : CC: k.b^n+1 & 2k.b^n+1
        5 : BiTwin: k.b^n & 2k.b^n
        6 : Twin/SG: k.b^n & 2k.b^n-1
        7 : Twin/CC: k.b^n & 2k.b^n+1
        8 : LP: k.b^n & k.b^(n-1)+1 & k.b^(n+1)+1
        9 : LM: k.b^n & k.b^(n-1)-1 & k.b^(n+1)-1
        10 : CC 1st kind
        11 : CC 2nd kind
        12 : BiTwin
        13 : AP: b^n+2k-1 (-2 if b is odd)
        14 : 3-tuple: k.b^n & k.b^n+5
        15 : 4-tuple: k.b^n, +5, +7
        16 : k.b^n+1 with k fixed
        17 : k.b^n-1 with k fixed
        18 : SG: k.b^n+1 & 2k.b^n+3
        19 : b^n+k with k fixed
        20 : b^n-k with k fixed
Please, enter your choice (0) : 17 ← k.b^n-1 の k 固定を選択
Base = 2
k = 1192 ← 1002 〜 2000 位の値を選択(小さい方が高速だが小さい値は探索が進んでいる)
nmin = 170000 ← TOP5000 に最低ランクインする n の値
nmax = 999999 ← 最大値
Type : k.b^n-1 with k fixed
Are you ready to start sieving ? (y) : 
CPU capabilities: CMOV: Supported. SSE2: Supported.
Using bitmap : allocating 101.3Kb of RAM...
...succeeded
830000 candidates
Sieving numbers with 51179 to 301033 decimal digits

ふるいが開始される。60秒毎に新しい値がふるえる位の範囲でよさそう。

http://pagesperso-orange.fr/jean.penne/index2.html

から llr をダウンロードする。

$ ./SLLRP4 -m
             Main Menu

         1.  Test/Input Data
         2.  Test/Continue
         3.  Test/Exit

         4.  Options/CPU
         5.  Options/Preferences
         6.  Advanced/Priority

         7.  Help/About
Your choice: 1

Input file (from NewPgen): : Newpgen
Output file (Results): : Results
Line number (1): 1

Accept the answers above? (Y): y
             Main Menu

         1.  Test/Input Data
         2.  Test/Continue
         3.  Test/Exit

         4.  Options/CPU
         5.  Options/Preferences
         6.  Advanced/Priority

         7.  Help/About
Your choice: 3

$ ./SLLRP4 -b

で探索が開始される。
ここで残った値もまだ素数に確定したわけではないので

http://www.utm.edu/research/primes/programs/gallot/

の proth.exe で最後の確認をする。

http://www.mycgiserver.com/~primesearch/mainmenu.jsp

に k= 251-1001 の範囲の探査をしているプロジェクトがある。
mersenne bench mark

Pcpu timesec/iermachineCPUclockprogramOS,compiler,othersnamedate
33,333,333N/A0.012GTX260Gforce1.4GhzMACLucasFFTWlinuxmsft2009.12.1
32,582,657N/A0.084PLAYSTATION3IBM Cell3.2GhzMACLucasFFTWlinuxmsft2009.8.11
20,996,011N/A0.091AOPEN AK86-LAthlon642087MhzmprimeVine linuxmsft2005.1.11
20,996,011N/A0.097IBM NetVista M41 Slim(6844-32J)P41.6GhzmprimeVine linux 2.4.22msft2004.7.3
20,996,011N/A0.507Compaq XP1000alpha 21264500MhzGlucasNetbsd 1.6.1 gcc 2.95.3msft2004.2.8
20,996,011N/A0.644Samsung UP1100alpha 21264600MhzMlucas-2.7b-genRedHat7msft2004.4.18
20,708,011N/A0.049GIGABYTEP4 Prescott2.8GhzPrime95Win2000msft2004.3.25
20,708,011N/A0.049GIGABYTEP4 Prescott2.8GhzmprimeFreeBSD 5.2-RELEASEmsft2004.3.26
20,708,011N/A0.097IBM NetVistaP41.6GhzmprimeVine linux 2.4.22msft2004.7.3
20,708,011N/A0.167GIGABYTEP4 Prescott2.8GhzGlucasWin2000 gcc 3.3.1msft2004.3.25
20,708,011N/A0.283IBM NetVistaP41.6GhzGlucasVine linux 2.4.22 gcc 2.95.3msft2004.7.3
2,976,221N/A0.010IBM NetVistaP41.6GhzmprimeFreeBSD 6.3-RELEASEmsft2008.3.19
2,976,221N/A0.011IBM NetVistaP41.6GhzmprimeVine linux 2.4.22msft2004.7.3
2,976,221N/A0.027GIGABYTEP4 Prescott2.8GhzGlucasWin2000 gcc 3.3.1msft2004.3.25
2,976,221N/A0.032Samsung UP1100alpha 21264600MhzGlucasRedHat7 gcc 3.3.3msft2004.4.18
2,976,221N/A0.034Compaq XP1000alpha 21264500MhzGlucasRedHat7 cccmsft2004.2.22
2,976,2211d13h18m0.044Compaq XP1000alpha 21264500MhzGlucasNetbsd 1.6.1 gcc 2.95.3msft2004.2.8
2,976,2211d22h7m0.047CompaqCeleron730MhzmprimeFreeBSD 4.7-RELEASEmsft2004.2.8
2,976,221N/A0.053Compaq XP1000alpha 21264500MhzMlucas-2.7b-genRedHat7msft2004.4.18
2,976,221N/A0.064Samsung UP1100alpha 21264600MhzMlucas-2.7b-genRedHat7msft2004.4.18
2,976,221N/A0.068CompaqCeleron730MhzGlucasFreeBSD 4.7-RELEASE gcc 2.95.4msft2004.3.25
2,976,221N/A0.069ApplePowerPC G4800MhzGlucasMac OS X 10.2.8 gcc 3.1msft2004.3.11
2,976,221N/A0.079IntelCeleron533MhzmprimeVine linux 2.4.22msft2004.4.06
2,976,221N/A0.103IntelCeleron533MhzGlucasVine linux 2.4.22 gcc 2.95.3msft2004.4.06
2,976,2213d21h45m0.113DEC EB164alpha 21164300MhzGlucasNetbsd 1.6.1 gcc 2.95.3msft2004.2.10
2,976,22126d15h51m0.840CompaqPentium/P54C120MhzGlucasFreeBSD 4.8-RELEASE gcc 2.95.4msft2004.3.9
216,0912m9s0.001AOPEN AK86-LAthlon642.0GhzGlucasGentoo linux gcc 3.4.3msft2004.12.13
216,0912m40sN/APLAYSTATION3IBM Cell3.2GhzMACLucasFFTWlinuxmsft2009.8.11
216,0915m28sN/AIBM NetVistaP41.6GhzGlucasVine linux 2.4.26 gcc 3.3.2msft2004.12.13
216,0916m9s0.002GIGABYTEP4 Prescott2.8GhzGlucasWin2000 gcc 3.3.1msft2004.3.11
216,09111m21s0.003ApplePowerPC G4800MhzGlucasMac OS X 10.2.8 gcc 3.1msft2004.3.11
216,0911h4m0.017CyrixVIA C3533MhzGlucasWin2000 gcc 3.3.1msft2004.3.11
216,0913h0.046CompaqPentium/P54C120MhzGlucasFreeBSD 4.8-RELEASE gcc 2.95.4msft2004.3.11

Riesel prime の探し方 

Riesel prime とは、 k*2n-1, odd k > 1 という形の素数です。

この形の素数を探すには、他の方のプロジェクトに参加するか個人で探す
ことになります。ここでは個人で探す場合について記述します。

過去に探索された範囲、および範囲の予約は、

http://karbon.npage.de/k_%3C_300_15525520.html
http://www.myjavaserver.com/~primesearch/
http://www.15k.org/

ですが、既知素数のリスト

http://primes.utm.edu/primes/download.php             

を見ると上のwwwページではまだ未探索、未発見とされている Riesel prime も沢山あることがわかります。

そこで確実に未探索の領域となると k の大きな値を探すことになりますが、
探索ソフトの llr は k が小さい方が高速(特に k>511) という特徴があります。

http://www.mersenneforum.org/index.php             

でも、いろいろな素数探索の宣言がされているので目を通すと良いようです。
楕円曲線法によるフェルマー数の素因数分解 

 Prime95 に楕円曲線法による素因数分解機能が追加されフェルマー数とメルセンヌ数
の分解プロジェクトもはじまりました。(Prim95 のメニューから選べます)
 乗算にアセンブラによる FFT と IBDWT を用いており従来よりかなり高速になって
います。新素因数の発見が期待されます。

srsieve の紹介 

srsieve は、newpgen より速いらしい篩いソフトです。

../../../www.geocities.com/g_w_reynolds/srsieve

k*2^n-1 で k 固定の場合を篩います。

$./srsieve -a -n 3e5 -N 1e6 -P 1e9 "196608*2^n-1"

196608*2^n-1 で 300000 < n < 1000000 の範囲を
1e9 まで篩う。

篩いは、1e9 まではやります。
(sr2sieve が、1e9 以降からしか開始できないため)

結果は、

$head sr_2.abcd 
ABCD 196608*2^$a-1 [300014] // Sieved to 1000000000 with srsieve
5
25
4
39
20
3
20
22
7

という abc 形式というファイルで得られます。

$./srfile -g sr_2.abcd

で Newpgen 出力フォーマットへ変換できます。

次により高速な sr2sieve で篩いを継続します。

$./sr2sieve --sse2 -i sr_2.abcd -p 1e9 -P 20e9
19978878169 | 196608*2^412147-1                                    
19996720703 | 196608*2^314314-1                                        
p=19996789363, 620774 p/sec, 8144 factors, 100.0% done, 3 sec/factor
sr2sieve 1.7.10 stopped: at p=20000000000 because range is complete.
Found factors for 8144 terms in 31803.538 sec. (expected about 8157.53)

real    530m3.550s
user    262m2.443s
sys     0m4.677s

1e9 から 20e9 までの篩いを指示しています。

$ ls -l
total 1476
-rw-r--r--  1 msft  others    6831 Apr 19 23:57 README
-rw-r--r--  1 msft  others    1726 Apr 19 23:57 README-threads
-rw-r--r--  1 msft  others  317570 Apr 20 09:34 factors.txt
-rwxr-xr-x  1 msft  others   84840 Apr 19 23:57 sr2sieve
-rw-r--r--  1 msft  others   44790 Apr 19 23:57 sr2sieve-1.7.10-linux-x86.tar.gz
-rw-r--r--  1 msft  others    3933 Apr 20 09:34 sr2sieve.log
-rw-r--r--  1 msft  others       0 Apr 19 23:57 sr2work.txt
-rw-r--r--  1 msft  others  156889 Apr 19 23:57 sr_2.abcd
-rwxr-xr-x  1 msft  others   28508 Apr 20 00:03 srfile
-rw-r--r--  1 msft  others  804976 Apr 19 23:57 srsieve.out

篩いの結果は factors.txt というファイルに篩われた数が入る形で得られます。
	

$./srfile --known-factors factors.txt sr_2.abcd

とすると sr_2.abcd から factors.txt にある数を除いて
srsieve.out に出力します。

$./srfile -g srsieve.out
$cat -n t17_b2_k196608.npg |tail
 56436  196608 999924
 56437  196608 999947
 56438  196608 999950
 56439  196608 999951
 56440  196608 999962
 56441  196608 999967
 56442  196608 999972
 56443  196608 999980
 56444  196608 999987
 56445  196608 999994

-g オプションで newpgen の形式に変換します。

$ ./srfile -a t17_b2_k196608.npg 
srfile 0.6.10 -- A file utility for srsieve.
Read 56444 terms for 1 sequence from NewPGen format file `t17_b2_k196608.npg'.
Wrote 56444 terms for 1 sequence to ABCD format file `sr_2.abcd'.

-a オプションで篩いの結果を除いた abc ファイルが作成されます。
これを入力にして sr2sieve を継続できます。

$sr2sieve 1.7.10 -- A sieve for multiple sequences k*b^n+/-1 and b^n+/-k.
Reading `sr_2.abcd' ...
WARNING: --pmin=20000000000 from command line overrides pmin=1000000000 from `sr_2.abcd'
Read 64588 terms for 1 sequence from ABCD format file `sr_2.abcd'.
Split 1 base 2 sequence into 34 base 2^60 subsequences.
Expecting to find factors for about 1085.54 terms in this range.
sr2sieve 1.7.10 started: 300014 <= n <= 999994, 20000000000 <= p <= 30000000000
概素数 

  フェルマーテストなどには全て合格するが、素数検査は桁数が大きすぎる
などで実行できない数を概素数という。

 PRP Records Probable Primes Top 10000

に１万桁以上の概素数が集められて更新されている。

最大の概素数は、

	(2^1148729+2^574365+1)/5  で 345802桁である。(2008/5 現在)

llr を使って Gaussian Mersenne norm(これは素数) を探索する
過程で発見されたようだ。

素数のレコードと比較すると小さく、最速のパソコンを数台お持ちなら
個人でも更新可能と思われる。

ただ、素数性が証明できる数ではいけないので、流通している高速なプログラムの
種類から、llr で 2^n+3 を探すと良いかもしれない。
Fermat Euler Teorem 

http://www2.cc.niigata-u.ac.jp/~takeuchi/tbasic/BackGround/Euler.html

  フェルマーの小定理を一般化したオイラーの定理は強力である。
  この定理は、約数についてフェルマーテストのような素数性検査を可能にする。
  例えば、メルセンヌ数が小さな素数で割れる場合に残った数についての素数性検査
  をメルセンヌ数に対するフェルマーテストのプログラムのループの回数だけを
  変えることで可能にする。

  何故、約数に直接フェルマーテストを行わないかというと(例えば pfgw では可能
  である)、約数では剰余の計算が大変になる、メルセンヌ数にはもともと高速な
  フェルマーテストプログラムが存在するからである。

  ただし、素数性を証明はできないので概素数となる。
まだ続く素数発見者に賞金 

http://www.eff.org/awards/coop

Through the EFF Cooperative Computing Awards, EFF will confer prizes of:
    * $50,000 to the first individual or group who discovers a prime number with at least 1,000,000 decimal digits 
    * $100,000 to the first individual or group who discovers a prime number with at least 10,000,000 decimal digits
    * $150,000 to the first individual or group who discovers a prime number with at least 100,000,000 decimal digits
    * $250,000 to the first individual or group who discovers a prime number with at least 1,000,000,000 decimal digits

まだ賞金は続きます。がんばりましょう。
量子コンピュータ 

  不確定性原理を応用した量子コンピュータというものが発明され、実際に
15 を 3 と 5 に素因数分解することに成功したようです。

  そういえば、ファイマンさんが仁科財団の招きで来日して「計算のエネル
ギーは何処まで減らせるか」という題で講演した時に、日本の物理学者達は
講演内容そっちのけでファイマンさんの過去の業績の質問してました。
メルセンヌ素数リンク

  • Mersenne Prime Search
  • 素数 TOP 5000
  • TOP 5000 登録を待っている素数
  • 2000年:多倍長整数計算研究記録 (Modified:2001-01-18)
  • prime.html
  • Web Teidan: June
  • 素因数分解 @ IDM
  • 哀朱姫小道具製作所
  • 離散フーリエ変換を用いた多倍長乗算の話
  • fftb_export.c.txt
  • FFT (高速フーリエ・コサイン・サイン変換) の概略と設計法
  • 12章 素因数分解アルゴリズム
  • The Largest Known Primes
  • Yves Gallot's Proth.exe: an implementation of Proth's Theorem for Windows
  • 最速の素数判定アルゴリズム
  • 円周率の計算に関する話題
  • PiHex- A distributed effort to calculate Pi
  • FFT と AGM による円周率計算プログラム
  • Mersenne Primes: History, Theorems and Lists
  • mersenneforum.org
  • 素因数分解工房：リンク
  • FMT(Fast Molulo Transformation) Method
  • The Top Twenty: Cullen
  • Pi
  • 各種分散コンピューティングプロジェクトの紹介
  • 各種分散コンピューティングプロジェクトの紹介(英語)
  • 素数探索プロジェクト(森井氏HP)
  • 47*2^n+1 Search
  • Pentium 最適化
  • Largest Known Primes
  • D.H.Lehmer
  • Lehmer's machine
  • Lehmer's machine
  • Lehmer's machine
  • Lehmer's machine
  • Prime
  • Discovery of a Lost Factoring Machine
  • 1926 Bicycle Chain Computing!
  • Miller's Primality Test
  • Generalized Fermat Primes Search
  • Mersenne Prime Digits and Names
  • Faster Integer Mulitplication
  • http://faginfamily.net/barry/Papers/Discrete%20Weighted%20Transforms.pdf
  • http://isthe.com/chongo/src/calc/lucas-calc