Pythonでメルセンヌ素数を計算する

素数にはいろいろな種類のものがあり、その中でもとりわけメルセンヌ素数は興味深いものがあります。

メルセンヌ数とメルセンヌ素数

メルセンヌ数

正の整数kに対して$2^k-1$で表すことができる数をメルセンヌ数といいます。ここでは、ここでは50億までのメルセンヌ数を計算し、リストにします。

k=1
temp=0
mersenn=[]
while temp<5_000_000_000:
temp=2**k-1
mersenn.append(temp)
k+=1
print(mersenn,end=' ')
print(len(mersenn),'個')

[1, 3, 7, 15, 31, 63, 127, 255, 511, 1023, 2047, 4095, 8191, 16383, 32767, 65535, 131071, 262143, 524287, 1048575, 2097151, 4194303, 8388607, 16777215, 33554431, 67108863, 134217727, 268435455, 536870911, 1073741823, 2147483647, 4294967295, 8589934591]

2のべき乗で計算するので、ある程度の大きさになると加速度的に桁数が多くなり、50億まで計算しても、メルセンヌ数はたかだか33個にとどります。

メルセンヌ数とメルセンヌ素数を分ける

メルセンヌ数のうち、素数になるものをメルセンヌ素数といいます。そこで、前節で求めたメルセンヌ数を素数（mersenn_p）と合成数（mersenn_c）に分けます。メルセンヌ数を見るときに、小さい数字から数えての順番が大切になります。そこで、k番目のメルセンヌ数を$M_k$として、kと$M_k$を要素とする辞書形式にします。

import sympy
mersenn_p={}
mersenn_c={}
for cnt,i in enumerate(mersenn):
if sympy.isprime(i):
mersenn_p[cnt+1]=i
else:
mersenn_c[cnt+1]=i
print(mersenn_p)
print(mersenn_c)
print('素　数:',len(mersenn_p),'個')
print('合成数:',len(mersenn_c),'個')

{2: 3, 3: 7, 5: 31, 7: 127, 13: 8191, 17: 131071, 19: 524287, 31: 2147483647}
{1: 1, 4: 15, 6: 63, 8: 255, 9: 511, 10: 1023, 11: 2047, 12: 4095, 14: 16383, 15: 32767, 16: 65535, 18: 262143, 20: 1048575, 21: 2097151, 22: 4194303, 23: 8388607, 24: 16777215, 25: 33554431, 26: 67108863, 27: 134217727, 28: 268435455, 29: 536870911, 30: 1073741823, 32: 4294967295, 33: 8589934591}
素　数: 8 個
合成数: 25 個

メルセンヌ数が素数の時には、その順序kも素数になります。

これは2に対するkを合成数とした場合、k=pq(p.qは1でない正の整数)と表すことができますが、k=pqが合成数の時には、$2^k - 1 = 2^{pq} - 1$

=$ (2^p-1)(2^{pq -p}+\cdots+2^{3p}+2^{2p}+2^p+1)$となり、メルセンヌ数も合成数になるためです。ただし、kが素数であっても、メルセンヌ素数でないものがあります。具体的にはk=11,23,29のときです。

リュカ・レーマーテストによるメルセンヌ素数の計算

リュカ・レーマー数列の計算

メルセンヌ素数かどうかを判断する際、リュカ＝レーマーテストという面白い方法があります。まず、リュカ・レーマー数列を計算します。次に、リュカ・レーマー数列とメルセンヌ数について判定式を用い、余りが0になったときはメルセンヌ素数となるというものです。

リュカ・レーマー数列

$S_1=4$

$S_{k+1}=S_k^2-2$

判定式　ｋは奇素数

${S_{k-1}}\equiv 0 \pmod{P_k}$

この式をもとに計算すると次のようになります。

k=4
ll=[]
for i in range(1,24):
ll.append(k)
k=k**2-2
ll

[4,
 14,
 194,
 37634,
 1416317954,
 2005956546822746114,
 4023861667741036022825635656102100994,
16191462721115671781777559070120513664958590125499158514329308740975788034,
262163465049278514526059369557563039213647877559524545911906005349555773831236935015956281848933426999307982418664943276943901608919396607297585154,
（９番目以降略）

8番目でとんでもない数値になってしまい、手持ちのPCではk=23くらいが限界です。

リュカ・レーマー判定法は、kを奇素数とすると、$S{k-1}$がk番目のメルセンヌ数で割り切れると、k番目のメルセンヌ数はメルセンヌ素数になるというものです。さっそく計算します。

メルセンヌ素数とリュカ・レーマー判定法

前節で求めたメルセンヌ素数をリュカ・レーマー判定法にあてはめます。

for k,p in mersenn_p.items():
if k >23:
break
print(k,p,ll[k-2]%p)

2 3 1
3 7 0
5 31 0
7 127 0
13 8191 0
17 131071 0
19 524287 0

kは奇素数なので3からはじまり、K=3,5,7のときは、メルセンヌ数は7，31、127と順調に素数になっていますが、このときの判定式は０と上手くいっています。次に$M_{11}=2047$はkが素数であるにもかかわらず素数になりませんが、この時の判定式は1736ととなりメルセンヌ素数でないことと一致します。このあと、k=13,17,19は素数になり、k=23のときは判定式が0でなくなり、素数でないことがわかります。

for k,c in mersenn_c.items():
if k >23:
break
print(k,c,ll[k-2]%c)

1 1 0
4 15 14
6 63 23
8 255 149
9 511 205
10 1023 95
11 2047 1736
12 4095 779
14 16383 4193
15 32767 20400
16 65535 25439
18 262143 221468
20 1048575 1036394
21 2097151 840107
22 4194303 1751891
23 8388607 6107895

次に、kが合成数（素数でない）の場合です。k=1は別として、これ以外は全て、判定式の余りは0になっていないことがわかります。とりわけk=11,23はkが素数であるにもかかわらず$p_k$がメルセンヌ素数でないものもしっかりと余りは0でなく、しっかり判別できていることがわかります。

メルセンヌ素数と完全数

メルセンヌ素数と完全数について面白い性質があります。k番目のメルセンヌ数がメルセンヌ素数pであるあとき、$2^{k-1}\times p$は偶数の完全数になるというものです。早速、プログラムを作成します。

約数の合計を計算する関数

整数numに対して、約数の個数と合計を計算します。

def divisor_cs(num):
cnt = 0
sigma = 0
for i in range(1, int(num**0.5)+1):
if num % i == 0:
cnt += 1
sigma += i
if i**2 == num:
continue
cnt +=1
sigma+=(int(num/i))
return cnt,sigma
_,sigma=divisor_cs(5)
sigma

Pythonによる約数の計算と素因数分解

メルセンヌ素数から完全数を計算

前節で求めたメルセンヌ素数の辞書mersenn_pに対して、$2^{k-1} \times p$を計算します。

for k,p in mersenn_p.items():
num=2**(k-1)*p
_,sigma=divisor_cs(num)
print (k,p,num,sigma-num)

2 3 6 6
3 7 28 28
5 31 496 496
7 127 8128 8128
13 8191 33550336 33550336
17 131071 8589869056 8589869056
19 524287 137438691328 137438691328
31 2147483647 2305843008139952128 2305843008139952128

メルセンヌ素数について約数の合計を計算すると、全て完全数であることがわかります。ちなみに、完全数は次の方法で求めることができます。

max=2305843008139952129
for i in range(1, max+1):
_,sigma=divisor_cs(i)
if sigma-i == i:
print (i)

手持ちのPCでは8128あたりまでが限界ですが、完全数は次のものがあります。

6, 28, 496, 8128, 33550336, 8589869056, 137438691328, 2305843008139952128

これは、メルセンヌ数から作った完全数と1対1対応していることがわかります。

メルセンヌ素数に$2^{k-1}$を掛けた値が完全数になるのは次の計算ができるためです。