１だけならぶ素数
このワークシートはMath by Codeの一部です。
＜レプ・ユニット素数ってなんだっけ？＞
レプ・ユニット数は、反復する１、つまり、repeatedUnitsの略した言い方。
レプ・ユニット素数はどのくらいあるでしょう。
111＝999/9=(1000-1)/9=(10^3-1)/9と表せるように、
一般に(10^n-1)/9で111...1のように1だけが並ぶ整数が作れます。
nが9973以下では、n=2,19,23,317,1031の５個だけだといわれています。
実験してみましょう。
juliaで素数判定関数をまた使います。次に、juliaでn=9999までの整数に対して
レプ・ユニット数(10^n-1)/9に素数filterをかけよう。
[IN]julia
#======================
N=[x for x in 2:30]
Rp=map(n->(n,((10^n-1)÷9)),N)
println(N)
println(Rp)
#======================
[OUT]
[2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30]
[(2, 11), (3, 111), (4, 1111), (5, 11111), (6, 111111), (7, 1111111), (8, 11111111), (9, 111111111), (10, 1111111111), (11, 11111111111), (12, 111111111111), (13, 1111111111111), (14, 11111111111111), (15, 111111111111111), (16, 1111111111111111), (17, 11111111111111111), (18, 111111111111111111), (19, -938527119301061290), (20, 862919959050249102), (21, 430646668853801415), (22, 207190227713669347), (23, 22264046724521073), (24, 222640467245210737), (25, 176766442039934975), (26, -281973810012822641), (27, -770099869716054016), (28, 497554224488149447), (29, 876265784057149667), (30, 564104918922807068)]
＜参考＞
（・素数リストの作り方、素数判定についてはこちら、
・リストにfilterをかける方法についてはこちら）
残念なら桁溢れが19けたで起きている。
juliaの整数intは19ケタが限界なので、
BigInt()を使って、ケタの限界をこえるように改造しよう。
1が３の倍数個あると３の倍数になるからとばそう。
1が2より大きい偶数個（2n)あると、左右にわけた1がn個ならんだ数の倍数になるからとばそう。
いやいや、
それだけではない。
また、m,nに対するレプ・ユユニット数Rm,Rnがあるとき、ｍがｎの倍数ならば、RmはRnの倍数になる。
すると、mが合成数ならば、Rmはｍの因子ｎのRnを因数にもつから、合成数になる。
レプユニット素数Rpを探すには、ｐが素数であることは必要だが、十分でなない。
これって、メルセンヌ素数と同様に、ｍが素数であることは必要だが十分ではないときと同じだね。
質問：m,nに対するレプ・ユユニット数Rm,Rnがあるとき、ｍがｎの倍数ならば、RmはRnの倍数になる。これはどうやって証明しますか。
レプユニット素数の候補数をならべよう。
Pythonの場合はJuliaより無名関数やmap,filterでリストを作ると、記述が長めになる。
だから、約数個数を求める関数divCountを作って、１行に詰め込みすぎないようにしよう。
2から100までの整数のリストMで、dicCountが2個になる整数Pが素数のリストとなる。
素数個の１を並べることで、レプ・ユニット素数の候補ができる。
#[IN]Python
#=================================
M=list(range(2,100))
divCount=lambda n:len(list(filter(lambda m:n % m==0, list(range(1,n+1)))))
P=[x for x in M if divCount(x)==2]
Rp=list(map(lambda n:(n,((10**n-1)//9)),P))
print(P)
print(Rp)
#==================================
Juliaでは、divcountを関数にしなくても、filterの中の無名関数で簡単にかける。
#[IN]Julia
#=================================
P=filter(n->length(filter( m -> n % m==0,1:n))==2, 2:100)
#べき乗のBigInt()は底につける。
Rp=map(n->(n,(BigInt(10)^n-1)÷BigInt(9)), P)
println(P)
println(Rp)
#=================================
[OUT]
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]
[(2, 11),
(3, 111),
(5, 11111),
(7, 1111111),
(11, 11111111111),
(13, 1111111111111),
(17, 11111111111111111),
(19, 1111111111111111111),
(23, 11111111111111111111111),
(29, 11111111111111111111111111111),
(31, 1111111111111111111111111111111),
(37, 1111111111111111111111111111111111111),
(41, 11111111111111111111111111111111111111111),
(43, 1111111111111111111111111111111111111111111),
(47, 11111111111111111111111111111111111111111111111),
(53, 11111111111111111111111111111111111111111111111111111),
(59, 11111111111111111111111111111111111111111111111111111111111),
(61, 1111111111111111111111111111111111111111111111111111111111111),
(67, 1111111111111111111111111111111111111111111111111111111111111111111),
(71, 11111111111111111111111111111111111111111111111111111111111111111111111),
(73, 1111111111111111111111111111111111111111111111111111111111111111111111111),
(79, 1111111111111111111111111111111111111111111111111111111111111111111111111111111),
(83, 11111111111111111111111111111111111111111111111111111111111111111111111111111111111),
(89, 11111111111111111111111111111111111111111111111111111111111111111111111111111111111111111),
(97, 1111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111)]
100以下のレプ・ユニット素数の候補が即座に表示できた。
juliaのいいところは,mapやfilterのリストへの操作が簡潔にかけるところだ。
デメリットはPythonのように、何もしなくても多桁（多倍長整数）対応になるわけではなく、
BigIntでくるむ必要があるため、可読性が落ちることだね。
2007年までに9個しか見つかっていないらしい。
pythonで検証してみる。
2,19の2個だけで数分かかる。
次は23の1個だけでもかなり時間がかかる。
317を見つけるには何かしら手段をさらに工夫が必要となりそうだ。
[IN]Python
#===========================================================
def isPrime(n):
lim = int(round(n**0.5))
for num in range(2,lim):
if n % num ==0:
return False
return True
M=list(range(2,20))
divCount=lambda n:len(list(filter(lambda m:n % m==0, list(range(1,n+1)))))
P=[x for x in M if divCount(x)==2]
R=list(map(lambda n:(n,((10**n-1)//9)),P))
Rp=list(filter(lambda n:isPrime((10**n-1)//9),P))
print(P)
print(R)
print(Rp)
#===========================================================
[2, 3, 5, 7, 11, 13, 17, 19]
[(2, 11), (3, 111), (5, 11111), (7, 1111111), (11, 11111111111), (13, 1111111111111), (17, 11111111111111111), (19, 1111111111111111111)]
[2, 19]
geogebraではn=19けた以上は探さずに終了している。
n=19というので、juliaが19けたが桁溢れしたことを思い出すね。
juliaでも同様なコードをかいて実行しても、時短にはならなかった。残念。
質問：rustでn=317までのレプ・ユニット素数を検索するには、どうしましょう。
素数リスト、レピュニット数リスト、レピュニット素数リストの順に作るのは同じですが、
最後の素数判定でスピードアップの工夫をします。
rustが高速なのと、素数判定が高速なので、あっという間に結果がでます。
約数2個の整数、つまり、素数のリストをつくりましょう。
2からlast＝330までの範囲を対象に
約数2個というフィルターをかけて集めます。
let prime: Vec = (2..=last)
.filter(|&n| (1..=n).filter(|&m| n % m == 0).count() == 2)
.collect();
println!("素数リスト: {:?}", prime);
次に１が素数個ならんだレピュニット数のリストを作りましょう。
素数を対象に、10の素数乗ー１を9で割った商に変換して、集めます。
// Rn = (10^n - 1) / 9
let repunit: Vec = prime
.iter()
.map(|&x| (ten.pow(x as u32) - &one) / &nine)
.collect();
println!("レピュニット数リスト（Rn）: {:?}", repunit);
最後に、レピュニット数Rnが素数かどうかを判定しますが。
１が330個ならんだ数が素数かどうかまで調べる計算は、
通常の素数判定（約数２個とか、２から平方根までに約数が1個もないとか）では
時間が膨大になり、PCが永遠の眠りにつきます。
そこで、乱択法を使いましょう。
円の面積をモンテカルロ法で調べるように、乱数を使って近似的に調べます。
Miller-Rabin法です。
nが素数ならば、
n - 1を２で割り切れるまで何度も割りましょう。
その割り算回数をk、最後の奇数の商をqとすると、n-１ = 2^k * q です。
2 <= a < n-1 の範囲でランダムな整数 a を選択して、
a^q ≡ 1 (mod n) または、
a^[(2^(i-1))*q] ≡-1(mod n)となるiが１からkの中に１つある。
このことの対偶から、
ランダムな整数 aで、mod nで
a^q !≡ 1 (mod n) and 1からkのどのiでもa^[(2^(i-1))*q] ≡-1(mod n)
ならば、nが合成数。このときのaを合成数の証人といいます。
Miller-Rabin法では証人aが4分3の合成数以上確保できるので、10回連続失敗する可能性は100万分の１程度と少ない（しかし、ゼロではない。）
合成数の判定がでないと、相当の確率で素数と言える。しかし、素数でない可能性もの残る。
アルゴリズムとしては、≡1, ≡-1のすべての判定をすべて不合格なら
最後に合成数と言えるので、「たぶん素数か」の問いに対してfalseを返す。
反対に、≡1, ≡-1が言えた時点で次の試行に進みどれも合格したときにtrueにたどり着く。
[IN]rust in Cargo.toml
[dependencies]
rand = "0.8"
num-bigint = { version = "0.4", features = ["rand"] } # rand機能を有効にする
num-traits = "0.2"
num-iter = "0.1"
[IN]rust in main.rs
use num_bigint::{BigInt, RandBigInt, ToBigInt};
use num_traits::{One, Zero};
use rand::thread_rng; // rng ではなく thread_rng
fn main() {
let one: BigInt = One::one();
let ten: BigInt = 10.to_bigint().unwrap();
let nine: BigInt = 9.to_bigint().unwrap();
let last = 330; //n=330までやってみる。
let prime: Vec< i64 > = (2..=last)
.filter(|&n| (1..=n).filter(|&m| n % m == 0).count() == 2)
.collect();
println!("素数リスト: {:?}", prime);
// Rn = (10^n - 1) / 9
let repunit: Vec< BigInt > = prime
.iter()
.map(|&x| (ten.pow(x as u32) - &one) / &nine)
.collect();
println!("レピュニット数リスト（Rn）: {:?}", repunit);
let repunit_prime: Vec< BigInt > = repunit
.into_iter()
.filter(|rn| is_prime_miller_rabin(rn, 10)) // 10回テスト
.collect();
println!("レピュニット素数: {:?}", repunit_prime);
}
fn is_prime_miller_rabin(n: &BigInt, repeat: usize) -> bool {
let zero = BigInt::zero();
let one = BigInt::one();
let two = 2.to_bigint().unwrap();
if n <= &one { return false; }
if n <= &3.to_bigint().unwrap() { return true; }
if n % &two == zero { return false; }
// n - 1 = 2^k * q を計算して、kとqを絞り出す。
let n_minus_one = n - &one;
let mut q = n_minus_one.clone();
let mut k = 0u32;
while &q % &two == zero {
q /= &two;
k += 1;
}
//乱択の準備
let mut rng = thread_rng();
//repeat回試行する。
'outer: for _ in 0..repeat {
// 2 <= a < n-1 の範囲でランダムな整数 a を選択
let a = rng.gen_bigint_range(&two, &n_minus_one);
// v = a^q mod n を求めてテストの準備
let mut v = a.modpow(&q, n);
// a^q ≡ 1 (mod n) なら次の試行へ
if v == one || v == n_minus_one {
continue 'outer;
}
// a^{ (2^i) * q } ≡ -1 (mod n) になる i があるかチェック
for _ in 0..k - 1 {
v = v.modpow(&two, n); // 順次 2乗していく
if v == n_minus_one {
continue 'outer;
}
}
// どの条件も満たさなければ、間違いなくnは合成数であり、aはその証人となる。
return false;
}
//たぶん素数
true
}
[OUT]cargo run --quiet
素数リスト: [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317]
レピュニット数リスト（Rn）: [11, 111, 11111, 1111111, 11111111111, 1111111111111, 11111111111111111, 1111111111111111111, 11111111111111111111111, 11111111111111111111111111111, 1111111111111111111111111111111, 1111111111111111111111111111111111111, 11111111111111111111111111111111111111111, 1111111111111111111111111111111111111111111, 11111111111111111111111111111111111111111111111, 11111111111111111111111111111111111111111111111111111, 11111111111111111111111111111111111111111111111111111111111, 1111111111111111111111111111111111111111111111111111111111111, 1111111111111111111111111111111111111111111111111111111111111111111, 11111111111111111111111111111111111111111111111111111111111111111111111, 1111111111111111111111111111111111111111111111111111111111111111111111111, 1111111111111111111111111111111111111111111111111111111111111111111111111111111, 11111111111111111111111111111111111111111111111111111111111111111111111111111111111, 11111111111111111111111111111111111111111111111111111111111111111111111111111111111111111, 1111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111, 11111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111, 1111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111, 11111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111, 1111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111, 11111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111, 1111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111, 11111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111, 11111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111, 1111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111, 11111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111, 1111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111, 1111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111, 1111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111, 11111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111, 11111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111, 11111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111, 1111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111, 11111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111, 1111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111, 11111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111, 1111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111, 1111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111, 1111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111, 11111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111, 1111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111, 11111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111, 11111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111, 1111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111, 11111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111, 11111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111, 11111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111, 11111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111, 1111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111, 1111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111, 11111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111, 1111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111, 11111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111, 1111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111, 11111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111, 1111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111, 11111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111]
レピュニット素数: [11, 1111111111111111111, 11111111111111111111111, 11111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111]