Lucas Lehmer Test

Mersenne素数

一般に(2の整数乗)-1をMersenne 数といい, 整数をpとすればMpと書く. Mpが素数の時はMersenne素数という. (2の合成数乗)-1は素数にならない. 例えば24-1は2進法で表すと1が4個並ぶ1111となる. これは11が2回並んでいるので, 11*100+11つまり11*101と分解できる. (2の合成数乗)-1は1がその合成数個並んでいる, それで上のように分解出来てしまう. つまり

2km-1=(2m-1)(2m(k-1) + 2m(k-2) + ... + 2m(k-k)).

(2の素数乗)-1だとこういう分解は出来ない. そこでpが素数のMersenne数Mpは素数になる可能性がある.

素数のpに対するMpの素数性をみてみると
 p     Mp  素数性
 2      3   素数
 3      7   素数
 5     31   素数
 7    127   素数
11   2047  合成数 23*89
13   8191   素数
このあとp=17, 19, 31, 61, ...と素数のMpが続く.

かつてM8191が素数かどうかという問題があって. 上の表を見てみよう. p=2の時Mp=3でM3は素数である. p=3の時Mp=7でM7は素数である. p=5の時Mp=31でM31は素数である. p=7の時Mp=127で, 上に記述はないが, M127も素数である. つまりMpが素数ならM(Mp)も素数になるかということである. 次のMersenne素数はp=13でMp=8191. そこで同様にM8191は素数だろうという仮説であった. しかしこれはイリノイ大学のIlliac I計算機で調べたところM8191は合成数であった.

既知のMersenne素数の表は例えばKnuth他著, 有沢誠他訳「コンピュータの数学」p.109参照 (訳書には原著の発刊より後に見つかったMersenne素数も記述してある.) それ以降は後述のホームページ参照.

Mpの素数性はLucas-Lehmerテストで調べる. (Knuth: The Art of Computer Programming, Vol 2 (3rd Ed) p.409など参照)

Lucas-Lehmerテスト

奇数のpに対し, S(1)=4, S(n+1)=S(n)2-2 mod 2p-1とし, S(p-1)=0 mod 2p-1ならMp=2p-1は素数である.
p=3のとき, Sは
4, 0 ->素数

p=5のとき, Sは
4, 14, 8, 0 ->素数

p=7のとき, Sは
4, 14, 67, 42, 111, 0 ->素数

p=11のとき, Sは
4, 14, 194, 788, 701, 119, 1877, 240, 282, 1736 ->合成数

p=13のとき, Sは
4, 14, 194, 4870, 3953, 5970, 1857, 36, 1294, 3470, 128, 0 ->素数
ところでこのLucasテストは2進法の計算機に向いている. つまり2p-1で割った剰余をとる計算であるが, S(n)は2進法でp桁なので, S(n)2-2は2p桁になって

<- 2進法で上のp桁 -><- 2進法で下のp桁 ->

のように出来た時, 上のp桁をp桁右シフトして下のp桁に合わせ
 <- 2進法で下のp桁 ->
+<- 2進法で上のp桁 ->
と足せばよい. これは何故かというと剰余の定理により, の多項式 f(x)をx-aで割った剰余はf(a)である. いまxを2pだとすると上の
<- 2進法で上のp桁 -><- 2進法で下のp桁 ->

<- 2進法で上のp桁 -> * x + <- 2進法で下のp桁 ->
のこと. これを x - 1 で割った剰余は
<- 2進法で上のp桁 -> * 1 + <- 2進法で下のp桁 ->
となる. したがって実際の割算はしないで済む.

Mersenne数の話題は下のホームページなどにも載っている.
http://www.mersenne.org/prime.htm
http://www.utm.edu/research/primes/mersenne.shtml

プログラム

このプログラムにはLucas Lehmer testの本体と, 印刷ルーチンP1および多倍長の乗算, 加算のサブルーチンがある.

多倍長の乗算サブルーチンはC(C(g))(下の桁)〜C(C(1g))(上の桁)に C(C(2g))(下の桁)〜C(C(3g))(上の桁)をかけ, 結果をC(4g)(下の桁)〜C(5g)(上の桁)におく. 各長語に符号のビットを除き35ビット詰めてある.

     a+l  a+l-2     a+i       a+2  a
          #### .... #### .... #### ####


     b+m  b+m-2     b+j       b+2  b
          #### .... #### .... #### ####


c+n  c+n-2  c+i+j+2 c+i+j     c+1  c
     #### .... #### #### .... #### ####
上の図でa+iとb+jをかけた2倍長の部分積を, その下半分がc+i+j, 上半分がc+i+j+2の語に足せばよい.

多倍長の加算サブルーチンはC(C(g))〜 と C(C(1g))〜を足し, C(2g)(下の桁)〜C(3g)(上の桁)におく.

Lucas Lehmer testの本体はまずpを読み, pが0なら停止, そうでなければMpの素数性をしらべ, 素数ならpを, 合成数ならcを印字し, 次のpを読み込む.

加算, 乗算の作業番地はd(390)からで, 60長語分あるからpとしては20 x 35 = 700まではテストできる. しかしこのアルゴリズムはO(p3)なので, p=200くらいまでが限界である.

0dから始まる10語の作業番地の内容は次のとおり.

0c:    -35      2
2c:      0      1
4c:      0      2
6c:      w      p  (w = 2 * (p div 35))
8c:  q bit strobe  (q = p mod 35)
0gからのパラメータの場所の内容は次のとおり,
0g:      d   d+w
2g:      d   d+w
3g:    d+w  d+3w

    0p=116,

0p: 長語整数印字サブルーチンP1
    0t:
    0m=0r,

  0m:   a   40  ┐      多倍長乗算ルーチン
  1     x   14r ┘ plant link
  2     p    4g ┐
  3     x    7r │
  4     s    5g │
  5     z   11r │ 積の場所をクリア
  6     p    2  │
  7     tl (  ) │
  8     p    7r │
  9     a   40  │
 10     jl   3r ┘
 11     p     g
 12     x   29r
 13     s    1g
 14     z  (  )   link
 15     p    2g
 16     x   30r
 17     s    3g 
 18     z   47r   内側のループ終り
 19     p   29r
 20     s     g   i 
 21     a   30r
 22     s    2g   i+j
 23     a    4g   c+i+j
 24     x   28r



 25     x   33r
 26     a   40
 27     x   36r
 28     ql (  )   c+i+j
 29     pl (  )   a+i
 30     wl (  )   b+j
 31     tl    o
 32     q    2
 33     tl (  )   c+i+j
 34     pl    o
 35     zl  44r
 36     al (  )   c+i+j+2
 37     ll2013
 38     tl    o ┐
 39     p   36r │
 40     x   33r │ carry処理
 41     a   40  │
 42     x   36r │
 43     jl  32r ┘
 44     p   30r
 45     a   40
 46     jl  16r
 47     p   29r
 48     a   40
 49     jl  12r

    0t:
    0a=0r,

    0a: a   40  ┐       多倍長加算ルーチン
  1     x   11r ┘ plant link
  2     p    2
  3     tl    o   carry クリア
  4     p     g
  5     x   13r
  6     p    1g
  7     x   14r
  8     p    2g
  9     x   18r
 10     s    3g
 11     z  (  )   link
 12     pl    o
 13     al (  )   a+i



 14     al (  )   b+i
 15     ll2013
 16     tl    o
 17     q    2
 18     tl (  )   c+i
 19     p   13r
 20     a   40
 21     x   13r
 22     p   14r
 23     a   40
 24     x   14r
 25     p   18r
 26     a   40
 27     jl   9r

    0t:
    0q=0r,
    32:212980q,     +を読んだ時, 0q番地へ飛ぶ定数
    0c=380,         定数と作業場所の原点
    0d=390,         多倍長の作業場所の原点

    0c: 262109,2,   左は-35,右は2
        0,1,        長語の1
        0,2,        長語の2
                    6cに語数w, 7cに素数p, 8cストロボ

    0q: pl    o   Lucas-Lehmer test
  1     zl    c   0の入力で停止
  2     tl   6c
  3     it
  4     jl    p   入力した素数を印字
  5     p  139r
  6     t     g
  7     t    2g
  8     p    7c   pをアキュムレータへ
  9     al    c   pから35を引き続け 語数とシフト数を得る
 10     kl   9r
 11     s     c   シフト数q
 12     x   35r
 13     x   82r
 14     x   96r
 15     x   99r
 16     l   18
 17     t    6c   語数wを格納
 18     a     g   d+w
 19     t    1g
 20     t    3g
 21     t    4g
 22     x   58r
 23     a    6c
 24     a    6c   d+3w
 25     t    5g
 26     s    6c   d+2w
 27     s   40    d+2w-2
 28     x   91r
 29     x   92r
 30     s    6c   d+w-2
 31     x   95r
 32     x  100r
 33     x  128r
 34     pl   2c
 35     l  (  )   q
 36     sl   2c
 37     tl   8c   ストロブ 右にqビットの1
 38     pl   4c
 39     al   4c
 40     tl    d   s0=4に設定
 41     p     g
 42     a   40  ┐
 43     x   47r │
 44     s    1g │
 45     z   50r │
 46     p    2  │
 47     tl (  ) │ 上位の作業場所を0にする
 48     p   47r │
 49     jl  42r ┘
 50     p    7c
 51     s    5c
 52     t    7c   カウンタをp-2に設定
 53     zl 117r   カウンタ=0なら素数の判定へ
 54     it
 55     jlm       多倍長乗算で s*s を計算
 56     p    1g   積はd+w〜d+3wにできる
 57     x   63r ┐
 58     pl (  ) │ d+w
 59     sl   4c │ s*s から 2 を引く
 60     rl  35  │
 61     tl    o │
 62     q    2  │
 63     tl (  ) │ d+w,d+w+2,...
 64     p   63r │
 65     a   40  │
 66     x   63r │
 67     x   70r │
 68     pl    o │ この辺carry処理
 69     zl  72r │



 70     al (  ) │
 71     jl  60r ┘
 72     p     g
 73     x   84r ┐
 74     p   91r │ d+2w-2
 75     x   80r │
 76     a   40  │ 積の上半分をd〜d+wへ
 77     x   81r │
 78     s    5g │ d+3w?
 79     z   90r │
 80     ql (  ) │ d+2w-2
 81     pl (  ) │ d+2w
 82     rl (  ) │
 83     q    2  │
 84     tl (  ) │ d
 85     p   84r │
 86     a   40  │
 87     x   84r │
 88     p   81r │
 89     jl  75r ┘
 90     pl   8c   積の中央部の残骸を消す
 91     cl (  )   d+2w-2
 92     tl (  )   d+2w-2
 93     it
 94     jla       多倍長加算へ
 95     pl (  )   d+w-2
 96     rl (  )   q
 97     tl    o   p桁からのあふれ
 98     p    2
 99     ll (  )   q
100     tl (  )   d+w-2 あふれ以外格納
101     p     g
102     x  106r
103     x  110r
104     pl    o ┐
105     zl 114r │
106     al (  ) │ d
107     ll2013  │ あふれがあれば足し込む
108     tl    o │
109     q    2  │
110     tl (  ) │ d
111     p  106r │
112     a   40  │
113     jl 102r ┘
114     p    7c
115     s    3c   カウンタ-1
116     jl  52r
117     p     g   素数判定(pビット全部1か)
118     x  121r
119     s   95r
120     z  127r
121     pl (  )   d
122     al   2c
123     kl 130r   全部が1ではない, 合成数
124     p  121r
125     a   40
126     jl 118r
127     pl   8c
128     bl (  )   d+w-2
129     zl 132r   最上語が全部1, 素数
130     p   91r   cをアキュムレータに置く
131     jl 133r
132     it        pをアキュムレータに置く
133     tllt      pまたはcを印字
134     p33
135     tllt      CR
136     l2 
137     tllt      LF
138     jl  40    次のデータを読みにいく
139     0d


        3+       ;p テストデータ 200までの素数
        5+       ;p ;pはMersenne素数
        7+       ;p
        11+
        13+      ;p
        17+      ;p
        19+      ;p
        23+
        29+
        31+      ;p
        37+
        41+
        43+
        47+
        53+
        59+
        61+      ;p
        67+
        71+
        73+
        79+
        83+
        89+      ;p



        97+
        101+
        103+
        107+     ;p
        109+
        113+
        127+     ;p
        131+
        137+
        139+
        149+
        151+
        157+
        163+
        167+
        173+
        179+
        181+
        191+
        193+
        197+
        199+
        0+