史上最大の素数が本当に素数なのか確かめてみた

史上最大の素数が更新されたらしいです。

wired.jp

こういうの見つける人ってすごいですよね、、、
でも
この数って本当に素数なの？
という疑問がわきます。

気になります。
どうしましょう。
確かめてみましょう。

素数を見つけるのは大変でも、ある数が素数かどうかを確かめるのはそんなに大変じゃないはず。

どうやって素数を判定するか

素数判定のアルゴリズムはいろいあるみたいですが*1、今回は以下のページを参考に、ミラー・ラビン素数判定法という確率的素数判定法を使います。
d.hatena.ne.jp

ミラー・ラビン素数判定法

ミラー–ラビン素数判定法 - Wikipedia
でも説明されていますが、以下のように判定を行います。

入力：を素数判定対象の奇数の整数、を繰り返し回数として
出力：
1. を2の冪乗（べきじょう）で割って、の形にする
2. 以下を回繰り返す:
(a)[1,n-1]の範囲からをランダムに選ぶ
(b)でかつ[0,s-1]の、範囲の全てのについて なら合成数とする。
3. 素数とする。

実際に史上最大の素数が本当に素数なのか判定する

判定はPythonでやります。
M74207281と名付けられた今回の素数は、こちらからダウンロードできるようです（zipで直接ダウンロードされます。）
解凍して早速開いてみましょう。
21.7MBのテキストファイルのようです。

なんと、100桁ごとに改行されています。
ここの処理はあとでPythonでやるとしましょう。
以下が今回のコードになります。

# -*- coding: utf-8 -*-

import random

#データの読み込み
f = open('C:/Users/Downloads/M74207281/M74207281.txt')
M74207281 = f.read()  # ファイル終端まで全て読んだデータを返す
M74207281 = M74207281.replace('\n','')
M74207281 = M74207281.replace('\r','') #この2行で改行を削除
f.close()

#ミラー・ラビン素数判定法を定義
def is_prime(q,k=50):
    q = abs(q)
    #計算するまでもなく判定できるも(2と2以外の偶数)
    if q == 2: return True
    if q < 2 or q&1 == 0: return False

    #n-1=2^s*dとし（但しaは整数、dは奇数)、dを求める
    d = (q-1)>>1
    while d&1 == 0:
        d >>= 1
    
    #判定をk回繰り返す
    for i in range(k):
        a = random.randint(1,q-1)
        t = d
        y = pow(a,t,q)
        #[0,s-1]の範囲すべてをチェック
        while t != q-1 and y != 1 and y != q-1: 
            y = pow(y,2,q)
            t <<= 1
        if y != q-1 and t&1 == 0:
            return False
    return True

それではM74207281が本当に素数なのか確かめます。

>>>is_prime(int(M74207281))
True

普通に三週間くらいかかってしまいましたが、Trueでした。
世界最大の素数は本当に素数でした。

ちなみに、このプログラムでis_prime関数を定義しているので、いろいろ遊べます。

>>> is_prime(4)
False
>>> is_prime(13)
True

*1:素数判定 - Wikipediaを参照