メルセンヌさんの素数
このワークシートはMath by Codeの一部です。
＜メルセンヌ素数ってなんだっけ？＞
1644年にフランスのメルセンヌさんが、「M(n)=2n-1型のn<=257の素数は11個だ」
と予想しました。（n=2,3,5,7,13,17,19,31,67,127,257）
そこから、この型の素数をメルセンヌ素数と呼ばれるようになったようです。
1876年にリュカ(Lucas)さんが、127以下のメルセンヌ素数Mnは
n=2,3,5,7,13,17,19,31までは同じですが、
n=61,89,107,127になることを証明したそうです。
1952年以降、その続きがn=512,607,1279,2203,..........と発見が続いています。
メルセンヌ型素数は、mが素数であることは必要条件。
なぜなら、mが合成数だと2m-1も合成数になるから。
このことは整式の分解で説明がつく。
ｍは素数であることが必要になる。しかし、十分とは限らないので要注意。
実際に調べてみよう。
juliaで素数判定関数を作る。
次に、juliaで2から31までの素数リストMを作る。
Mの要素mに対して2m-1が素数になるものをfilterで取り出そう。
＜参考＞
（・素数リストの作り方、素数判定についてはこちら、
・リストにfilterをかける方法についてはこちら）
[IN]julia
#======================
function isPrime(n)
lim = Int(round(n^0.5))
for num in 2:lim
if BigInt(n) % BigInt(num) ==0
return false
end
end
return true
end
M=filter(n->length(filter( m -> n % m==0,1:n))==2, 2:31)
Ms=filter(m->isPrime(2^m-1),M)
println(M)
println(Ms)
#======================[OUT]
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31]
[2, 3, 5, 7, 13, 17, 19, 31]
このように少し調べただけでも
2m-1が素数になる素数ｍを限定できたね。
geogebraでも試してみよう。
Juliaと同じロジックで、ほぼ同じような記述で調べることができる。
メルセンヌ素数判定法
＜リュカ数＞
Mp＝2^p-1
Spは漸化式で、S0=4, Sk=Sk-12-2.
p番めのメルセンヌ数をMp,リュカ数をSpとする。
Sp-2がMpで割り切れるときに、Mpは素数になる
ことが必要十分条件だ。
Spは2乗して２をひくことを繰り返すので、急激に桁が大きくなると予想されるね。
それでも多桁なのでJuliaで多桁用にBigInt()でかこってリュカ数(Lucas)を作ってみよう
#Julia
#[IN]
#====================
#リュカ数を作る
Lucas=BigInt(4^2-2)
S=[Lucas]
for i in 1:8
Lucas=BigInt(Lucas^2-2)
push!(S,Lucas)
end
println(S)
#====================
[OUT]
BigInt[14, 194, 37634, 1416317954, 2005956546822746114, 4023861667741036022825635656102100994, 16191462721115671781777559070120513664958590125499158514329308740975788034, 262163465049278514526059369557563039213647877559524545911906005349555773831236935015956281848933426999307982418664943276943901608919396607297585154, 68729682406644277238837486231747530924247154108646671752192618583088487405790957964732883069102561043436779663935595172042357306594916344606074564712868078287608055203024658359439017580883910978666185875717415541084494926500475167381168505927378181899753839260609452265365274850901879881203714]
＜リュカ・レーマーテスト＞
Sp-2がMpの倍数であることがMpはメルセンヌ素数であることの必要十分条件だった。
S=[14, 194, 37634, 1416317954, 2005956546822746114]
M=[1,3,7,15,31,63,.....]
S1 % M3=14 / 7=2だから、M3=7は素数
S2 % M4=194 % 15 !==0だから、M4は素数ではない。
S3 % M5= 37634 / 31 =1214だから、M5は素数。
こうして、２つ番号ちがいのリュカ数列をメルセンヌ数列で割り切れると
素数と判定できる。割り切れるかどうかは、剰余が０になるかどうかだから、Spそのものではなく、SpをMpで割った剰余で代用して、
けたの爆発を防止しよう。
#Julia
#[IN]
#====================
target=127
M=[]
for x in 3:target
push!(M,BigInt(2)^x - 1)
end
#println(M)
function LucaModM(x)
#x番目のメルセンヌ数を法とした剰余でリュカ数列を作る。
res=4^2-2
for i in 2:x
res = BigInt(res^2-2) % BigInt(M[x])
end
return res
end
for n in 3:target-2
if LucaModM(n)==0
println("M",n+2,"=2^",n+2,"-1=",M[n]," is Prime")
end
end#====================
[OUT]
M5=2^5-1=31 is Prime
M7=2^7-1=127 is Prime
M13=2^13-1=8191 is Prime
M17=2^17-1=131071 is Prime
M19=2^19-1=524287 is Prime
M31=2^31-1=2147483647 is Prime
M61=2^61-1=2305843009213693951 is Prime
M89=2^89-1=618970019642690137449562111 is Prime
M107=2^107-1=162259276829213363391578010288127 is Prime
M127=2^127-1=170141183460469231731687303715884105727 is Prime
質問：リュカ・レーマーテストをrustでかくとどうなりますか。
juliaでは、
・メルセンヌ数列Mpを作る。（pを素数に限らずすべて）
・リュカ数sをメルセンヌ数で割りながら最後の剰余を返すチェック関数。
・剰余の結果を走査して、チェック関数が０になったらメルセンヌ数が素数だと表示する。
この3段階でやってました。
target=127
M=[]
for x in 3:target
push!(M,BigInt(2)^x - 1)
end
#println(M)
function LucaModM(x)
#x番目のメルセンヌ数を法とした剰余でリュカ数列を作る。
res=4^2-2
for i in 2:x
res = BigInt(res^2-2) % BigInt(M[x])
end
return res
end
for n in 3:target-2
if LucaModM(n)==0
println("M",n+2,"=2^",n+2,"-1=",M[n]," is Prime")
end
rustでは、
・Mpを作るときにpが素数であるかどうかを調べてから回すようします。
計算量を減らすために√ｐまで調べます。is_prime関数です。
・メルセンヌ数の作成とリュカ数によるチェックを両方やりましょう。
そのチェック関数がlucas_lehmerです。
・最後に走査するループでチェック関数を呼べば、ループを減らせます。
スピードを実感するために、ターゲットを127ではなく、2203でやってみよう。
わりとすぐに結果が表示されるでしょう。
[IN]rust in Cargo.toml
[dependencies]
num-traits = "0.2"
[IN]rust in main.rs
use num_bigint::BigInt;
use num_traits::{One, Zero};
fn lucas_lehmer(p: u32) -> bool {
if p == 2 { return true; }
//p番目のメルセンヌ数Mp=2^p-1
let m_p = (BigInt::from(2).pow(p)) - BigInt::one();
//リュカ数の初項は4
let mut s = BigInt::from(4);
// p-2 回繰り返す
for _ in 0..p - 2 {
//□番目のリュカ数をp番目のメルセンヌ数で割った余りを求める。
//それで、次のリュカ数とする。 p-2番目まで繰り返す。s = (s^2 - 2) mod M_p
s = (&s * &s - BigInt::from(2)) % &m_p;
}
//割り切れたらtrueを返す。
s == BigInt::zero()
}
fn main() {
let target = 2203;
for p in 2..=target {
// メルセンヌ素数であるためには p 自体が素数である必要がある。
if is_prime(p) {
//メルセンヌ素数であるために リュカ数Sp-2がMpの倍数であることが必要十分。
if lucas_lehmer(p) {
println!("M{} is Prime", p);
}
}
}
}
// 指数 p が素数かどうかを調べるための関数
fn is_prime(n: u32) -> bool {
if n < 2 { return false; }
for i in 2..=((n as f64).sqrt() as u32) {
if n % i == 0 { return false; }
}
true
}
[OUT] コマンドラインでビルドと実行。cargo run --quiet
M2 is Prime
M3 is Prime
M5 is Prime
M7 is Prime
M13 is Prime
M17 is Prime
M19 is Prime
M31 is Prime
M61 is Prime
M89 is Prime
M107 is Prime
M127 is Prime
M521 is Prime
M607 is Prime
M1279 is Prime
M2203 is Prime