1だけならぶ素数
このワークシートはMath by Codeの一部です。
<レプ・ユニット素数ってなんだっけ?>
レプ・ユニット数は、反復する1、つまり、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の5個だけだといわれています。
実験してみましょう。
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が3の倍数個あると3の倍数になるからとばそう。
1が2より大きい偶数個(2n)あると、左右にわけた1がn個ならんだ数の倍数になるからとばそう。
いやいや、
それだけではない。
また、m,nに対するレプ・ユユニット数Rm,Rnがあるとき、mがnの倍数ならば、RmはRnの倍数になる。
すると、mが合成数ならば、Rmはmの因子nのRnを因数にもつから、合成数になる。
レプユニット素数Rpを探すには、pが素数であることは必要だが、十分でなない。
これって、メルセンヌ素数と同様に、mが素数であることは必要だが十分ではないときと同じだね。
質問:m,nに対するレプ・ユユニット数Rm,Rnがあるとき、mがnの倍数ならば、RmはRnの倍数になる。これはどうやって証明しますか。
レプユニット素数の候補数をならべよう。
Pythonの場合はJuliaより無名関数やmap,filterでリストを作ると、記述が長めになる。
だから、約数個数を求める関数divCountを作って、1行に詰め込みすぎないようにしよう。
2から100までの整数のリストMで、dicCountが2個になる整数Pが素数のリストとなる。
素数個の1を並べることで、レプ・ユニット素数の候補ができる。
#[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でも同様なコードをかいて実行しても、時短にはならなかった。残念。
質問:juliaでn=19までのレプ・ユニット素数を検索するにはどうかいたらよいでしょうか。