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+