[size=150][b][size=150]このワークシートは[url=https://www.geogebra.org/m/twxxx3yq]Math by Code[/url]の一部です。[br][/size][br]＜メルセンヌ素数ってなんだっけ？＞[br][br][/b]1644年にフランスのメルセンヌさんが、「M(n)=2[sup]n[/sup]-1型のn<=257の素数は11個だ」[br]と[color=#9900ff]予想[/color]しました。（n=2,3,5,7,13,17,19,31[color=#ff0000],67,[/color]127,[color=#ff0000]257[/color]）[br]そこから、この型の素数をメルセンヌ素数と呼ばれるようになったようです。[br]1876年にリュカ(Lucas)さんが、127以下のメルセンヌ素数Mnは[br]n=2,3,5,7,13,17,19,31までは同じですが、[br]n=[color=#0000ff]61,89,107[/color],127になることを証明したそうです。[br]1952年以降、その続きがn=[color=#0000ff]512,607,1279,2203[/color],..........と発見が続いています。[br][br]メルセンヌ型素数は、mが素数であることは必要条件。[br]なぜなら、mが合成数だと2[sup]m[/sup]-1も合成数になるから。[br]このことは整式の分解で説明がつく。[br][b]ｍは素数であることが必要になる。しかし、十分とは限らないので要注意。[br]実際に調べてみよう。[br]juliaで素数判定関数を作る。[br]次に、juliaで2から31までの素数リストMを作る。[br]Mの要素mに対して2[sup]m[/sup]-1が素数になるものをfilterで取り出そう。[br]＜参考＞[br]（・[b][color=#0000ff]素数リスト[/color][/b]の作り方、[color=#0000ff][b]素数判定[/b][/color]については[url=https://www.geogebra.org/m/jxmueqdt]こちら、[br][/url]　・リストに[color=#0000ff][b]filterをかける方法[/b][/color]については[url=https://www.geogebra.org/m/vre8rh6h]こちら[/url]）[br][br][br][/b][color=#0000ff][IN]julia[br][/color][/size][color=#0000ff]#======================[br][/color]function isPrime(n)[br] lim = Int(round(n^0.5))[br] for num in 2:lim[br] if BigInt(n) % BigInt(num) ==0 [br] return false[br] end [br] end[br] return true[br]end[br]M=filter(n->length(filter( m -> n % m==0,1:n))==2, 2:31)[br]Ms=filter(m->isPrime(2^m-1),M)[br]println(M)[br]println(Ms)[br][color=#0000ff]#======================[/color][OUT][br][2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31][br][2, 3, 5, 7, 13, 17, 19, 31][br][br]このように少し調べただけでも[br]2[sup]m[/sup]-1が素数になる素数ｍを限定できたね。[br][br]geogebraでも試してみよう。[br]Juliaと同じロジックで、ほぼ同じような記述で調べることができる。

＜リュカ数＞[br]Mp＝2^p-1[br]Spは漸化式で、S[sub]0[/sub]=4, S[sub]k[/sub]=S[sub]k-1[/sub][sup]2[/sup]-2.[br][br]p番めのメルセンヌ数をMp,リュカ数をSpとする。[br][color=#0000ff][b][size=200]S[sub]p-2[/sub]がM[sub]p[/sub]で割り切れるときに、Mpは素数になる[br][/size][/b][/color]ことが必要十分条件だ。[br]Spは2乗して２をひくことを繰り返すので、急激に桁が大きくなると予想されるね。[br]それでも多桁なのでJuliaで多桁用にBigInt()でかこってリュカ数(Lucas)を作ってみよう[br][color=#0000ff]#Julia[br][/color][color=#0000ff]#[IN][br][/color][color=#0000ff]#====================[br][/color]#リュカ数を作る[br]Lucas=BigInt(4^2-2)[br]S=[Lucas][br]for i in 1:8[br] Lucas=BigInt(Lucas^2-2)[br] push!(S,Lucas)[br]end[br]println(S)[br][color=#0000ff]#====================[br][OUT][br][/color]BigInt[14, 194, 37634, 1416317954, 2005956546822746114, 4023861667741036022825635656102100994, 16191462721115671781777559070120513664958590125499158514329308740975788034, 262163465049278514526059369557563039213647877559524545911906005349555773831236935015956281848933426999307982418664943276943901608919396607297585154, 68729682406644277238837486231747530924247154108646671752192618583088487405790957964732883069102561043436779663935595172042357306594916344606074564712868078287608055203024658359439017580883910978666185875717415541084494926500475167381168505927378181899753839260609452265365274850901879881203714][br][br]＜リュカ・レーマーテスト＞[br]S[sub]p-2[/sub]がMpの倍数であることがMpはメルセンヌ素数であることの必要十分条件だった。[br] S=[14, 194, 37634, 1416317954, 2005956546822746114][br]M=[1,3,7,15,31,63,.....][br]S1 % M3=14 / 7=2だから、M3=7は素数[br]S2 % M4=194 % 15 !==0だから、M4は素数ではない。[br]S3 % M5= 37634 / 31 =1214だから、M5は素数。[br]こうして、２つ番号ちがいのリュカ数列をメルセンヌ数列で割り切れると[br]素数と判定できる。割り切れるかどうかは、剰余が０になるかどうかだから、[color=#0000ff][b][u]Spそのものではなく、SpをMpで割った剰余で代用[/u][/b][/color]して、[br]けたの爆発を防止しよう。[br][color=#0000ff]#Julia[br][/color][color=#0000ff]#[IN][br][/color][color=#0000ff]#====================[/color][br]target=127[br]M=[][br]for x in 3:target[br] push!(M,BigInt(2)^x - 1)[br]end[br]#println(M)[br]function LucaModM(x)[br] #x番目のメルセンヌ数を法とした剰余でリュカ数列を作る。[br] res=4^2-2[br] for i in 2:x[br] res = BigInt(res^2-2) % BigInt(M[x])[br] end[br] return res[br]end[br]for n in 3:target-2[br] if LucaModM(n)==0[br] println("M",n+2,"=2^",n+2,"-1=",M[n]," is Prime")[br] end[br]end[color=#0000ff]#====================[/color][br][color=#0000ff][OUT][br]M5=2^5-1=31 is Prime[br]M7=2^7-1=127 is Prime[br]M13=2^13-1=8191 is Prime[br]M17=2^17-1=131071 is Prime[br]M19=2^19-1=524287 is Prime[br]M31=2^31-1=2147483647 is Prime[br]M61=2^61-1=2305843009213693951 is Prime[br]M89=2^89-1=618970019642690137449562111 is Prime[br]M107=2^107-1=162259276829213363391578010288127 is Prime[br]M127=2^127-1=170141183460469231731687303715884105727 is Prime[/color][br][br][br]

[color=#9900ff][b][u][size=150]質問：リュカ・レーマーテストをrustでかくとどうなりますか。[br][/size][/u][/b][/color][br]juliaでは、[br]・メルセンヌ数列Mpを作る。（pを素数に限らずすべて）[br]・リュカ数sをメルセンヌ数で割りながら最後の剰余を返すチェック関数。[br]・剰余の結果を走査して、チェック関数が０になったらメルセンヌ数が素数だと表示する。[br]この3段階でやってました。[br][br][color=#0000ff][b]target=127[br][/b]M=[][br]for x in 3:target[br] push!(M,BigInt(2)^x - 1)[br]end[br]#println(M)[br]function [b]LucaModM[/b](x)[br] #x番目のメルセンヌ数を法とした剰余でリュカ数列を作る。[br] res=4^2-2[br] for i in 2:x[br] res = BigInt(res^2-2) % BigInt(M[x])[br] end[br] return res[br]end[br]for n in 3:target-2[br] if LucaModM(n)==0[br] println("M",n+2,"=2^",n+2,"-1=",M[n]," is Prime")[br] end[br][br][/color]rustでは、[br]・Mpを作るときにpが素数であるかどうかを調べてから回すようします。[br]計算量を減らすために√ｐまで調べます。[b]is_prime関数[/b]です。[br]・メルセンヌ数の作成とリュカ数によるチェックを両方やりましょう。[br]そのチェック関数が[b]lucas_lehmer[/b]です。[br]・最後に走査するループでチェック関数を呼べば、ループを減らせます。[br]スピードを実感するために、ターゲットを[b]127ではなく、2203[/b]でやってみよう。[br]わりとすぐに結果が表示されるでしょう。[br][br][color=#0000ff][IN]rust in Cargo.toml[br][dependencies][br]num-traits = "0.2"[br][IN]rust in main.rs[br][br]use num_bigint::BigInt;[br]use num_traits::{One, Zero};[br][br]fn [b]lucas_lehmer[/b](p: u32) -> bool {[br] if p == 2 { return true; }[br] //p番目のメルセンヌ数Mp=2^p-1[br] let m_p = (BigInt::from(2).pow(p)) - BigInt::one();[br] //リュカ数の初項は4[br] let mut s = BigInt::from(4);[br][br] // p-2 回繰り返す[br] for _ in 0..p - 2 {[br] //□番目のリュカ数をp番目のメルセンヌ数で割った余りを求める。[br] //それで、次のリュカ数とする。 p-2番目まで繰り返す。s = (s^2 - 2) mod M_p[br] s = (&s * &s - BigInt::from(2)) % &m_p;[br] }[br] //割り切れたらtrueを返す。[br] s == BigInt::zero()[br]}[br][br]fn [b]main[/b]() {[br] let target = 2203;[br] for p in 2..=target {[br] // メルセンヌ素数であるためには p 自体が素数である必要がある。[br] if is_prime(p) {[br] //メルセンヌ素数であるために リュカ数Sp-2がMpの倍数であることが必要十分。[br] if lucas_lehmer(p) {[br] println!("M{} is Prime", p);[br] }[br] }[br] }[br]}[br][br]// 指数 p が素数かどうかを調べるための関数[br]fn [b]is_prime[/b](n: u32) -> bool {[br] if n < 2 { return false; }[br] for i in 2..=((n as f64).sqrt() as u32) {[br] if n % i == 0 { return false; }[br] }[br] true[br]}[br][OUT] コマンドラインでビルドと実行。cargo run --quiet[br][b]M2 is Prime[br]M3 is Prime[br]M5 is Prime[br]M7 is Prime[br]M13 is Prime[br]M17 is Prime[br]M19 is Prime[br]M31 is Prime[br]M61 is Prime[br]M89 is Prime[br]M107 is Prime[br]M127 is Prime[br]M521 is Prime[br]M607 is Prime[br]M1279 is Prime[br]M2203 is Prime[br][br][br][/b][/color]