果てしない素数
このワークシートはMath by Codeの一部です。
果てしない夢
無限の可能性
素数を直線上に並べていけば、終わりがないですね。
今回は、素数が無限にあることにかかわることをさぐっていこう。
１．ユークリッドの証明
ユークリッドの証明
たとえば、
２の倍数に１をたした数は２の倍数ではない。
３の倍数に１をたした数は３の倍数ではない。
５の倍数に１をたした数は５の倍数ではない。
これはだれもが認める事実です。
＜ユークリッドさん＞
ユークリッドはこれらの事実を利用します。
2×3×5は２、３、５の倍数ですから、
2×3×5＋１は、２の倍数でも、３の倍数でも、５の倍数でもないことになります。
ユークリッド「原論」の論法はこうでした。
素数が有限個しかなく、そのリストが２，３、...、ｐのように、ｐで終わるとする。
N＝２×３×…×ｐ＋１は、すべての素数２からｐで割り切れません。
つまりNは素数を使って合成できないので、Nは素数です。
しかし、Nは素数リストの最大数ｐより大きいので、矛盾します。
だから、素数は有限個とする仮定が誤りだから、素数は無限にある。
このユークリッドさんの論法は現実離れしてますね。
そもそも、すべての素数という物自体が確定できないのに、それを仮定する思考だからです。
プログラミングしようとしても、すべての素数の集合は作れないので計算不可能ですね。
本当に神の視点に立った論理ですね。
たとえば、素数＝ [2,3,5,7,11,13,17,19,23]。
これがすべての素数だとして、Nを計算するとNは合成数になってしまいます。
「すべての素数」といっておきながら、途中の素数をすっとばして、
でかいNを作っているわけですから、Nが素数にならない可能性が出てきますね。
だから、この「すべての素数」という、現実的には、できもしないことを前提にした論理ですから、
現実の有限数でたしかめようとしても破綻しますね。
無限個なものを有限とあえて偽っているのですから。
[IN]Python
#Euclid
from sympy import isprime
from operator import mul
from functools import reduce
def prime_euc():
Ps = [2,3,5,7,11,13,17,19,23]
N = reduce(mul,Ps)+1
print(isprime(N),end=":")
print(get_divisors(N))
Ps = Ps + [N]
return Ps
print(prime_euc())
[OUT]
False:[1, 317, 703763, 223092871]
[2, 3, 5, 7, 11, 13, 17, 19, 23, 223092871]
223092871は素数ではないが、最小の素因数317はPsの最大の素数23より大きい。
ここに着目すれば、少数の素数リストからスタートして、素数を追加することはできそうだ。
＜ピタゴラスとマリンさん＞
ユークリッドの証明法のルーツはピタゴラスの論法だという説があります。
その論法を明確化して、おまけをつけたマリンという人がいます。
素数リストをPsとしましょう。
（ベース）Psに2を入れます。
（反復）N=（Psにあるすべての素数の積）＋１。Nの１より大きい最小約数MをPsに加える。
MがPsの素数とかぶらないので、この反復は無限に続きますから、素数は無限リストになります。
「順番はともかくこのリストにすべて素数が入る」という仮説があります。
課題：ピタゴラス・マリンの素数列をつくるにはどうしたらよいですか。
ジェネレータや関数プログラミングツールを使うことで、日本語をただ英語に直す感じで
コードをかいてみよう。
[IN]Python
from operator import mul
from functools import reduce
def prime_seq():
Ps = [2]
while True:
yield Ps
N = reduce(mul,Ps)+1
for i in range(2, N+1):
if N % i == 0:
M = i # Mは1より大きいNの最小約数で、素数
break
Ps = Ps + [M]
generator = prime_seq()
for _ in range(8):
print(next(generator))
[OUT]
[2]
[2, 3]
[2, 3, 7]
[2, 3, 7, 43]
[2, 3, 7, 43, 13]
[2, 3, 7, 43, 13, 53]
[2, 3, 7, 43, 13, 53, 5]
[2, 3, 7, 43, 13, 53, 5, 6221671]
９ステップ目になかなかいかないので、sympyにたよってみよう。
[IN]Python
from functools import reduce
from operator import mul
from sympy import factorint
def prime_seq_fast():
Ps = [2]
while True:
yield Ps
N = reduce(mul, Ps) + 1
# factorint(N) は {素因数: 個数} の辞書を高速に返します
factors = factorint(N)
# 最小の素因数（最小の約数 M）を取得
M = min(factors.keys())
Ps = Ps + [M]
generator = prime_seq_fast()
# 少し多めに回してみます
for i in range(17):
print(f"Step {i+1}: {next(generator)}")
[OUT]
Step 1: [2]
Step 2: [2, 3]
Step 3: [2, 3, 7]
Step 4: [2, 3, 7, 43]
Step 5: [2, 3, 7, 43, 13]
Step 6: [2, 3, 7, 43, 13, 53]
Step 7: [2, 3, 7, 43, 13, 53, 5]
Step 8: [2, 3, 7, 43, 13, 53, 5, 6221671]
Step 9: [2, 3, 7, 43, 13, 53, 5, 6221671, 38709183810571]
Step 10: [2, 3, 7, 43, 13, 53, 5, 6221671, 38709183810571, 139]
Step 11: [2, 3, 7, 43, 13, 53, 5, 6221671, 38709183810571, 139, 2801]
Step 12: [2, 3, 7, 43, 13, 53, 5, 6221671, 38709183810571, 139, 2801, 11]
Step 13: [2, 3, 7, 43, 13, 53, 5, 6221671, 38709183810571, 139, 2801, 11, 17]
Step 14: [2, 3, 7, 43, 13, 53, 5, 6221671, 38709183810571, 139, 2801, 11, 17, 5471]
Step 15: [2, 3, 7, 43, 13, 53, 5, 6221671, 38709183810571, 139, 2801, 11, 17, 5471, 52662739]
Step 16: [2, 3, 7, 43, 13, 53, 5, 6221671, 38709183810571, 139, 2801, 11, 17, 5471, 52662739, 23003]
Step 17: [2, 3, 7, 43, 13, 53, 5, 6221671, 38709183810571, 139, 2801, 11, 17, 5471, 52662739, 23003, 30693651606209]
すざまじい素数の波ですね。津波がきたとおもったら、やさしいさざ波になったり、また大波だ。
２．素数の隙間
素数の隙間と言えば
双子素数が有名だ。
３と５
１１と１３
１７と１９
のように、素数が最低の隙間、１つの整数があく２つの素数だ。
ではその反対、素数の空白地帯を考えてみよう。
２×３×５+2=32
２×３×５+3=33
２×３×５+4=34
２×３×５+5=35
２×３×５+6=36
２×３×５×７+2 =212
２×３×５×７+3 =213
２×３×５×７+4 =214
２×３×５×７+5 =215
２×３×５×７+6 =216
２×３×５×７+7 =217
２×３×５×７+8 =218
２×３×５×７+9 =219
２×３×５×７+10=220
課題：素数の隙間をもっと長く作るにはどうしたらよいでしょうか。
たとえばｐを1009にすると、2から1001までの整数は素数リスト（2から1009）のどれかで割り切れますね。だから、N=product(２から1009)として、iをrange(2,1002)で動かしたN+iは素数になりません。
でも、これだと1000個の非素数列ができますが、長すぎます。
100個を出すコードを作ってみよう。
[IN]Python
from sympy import isprime
from operator import mul
from functools import reduce
def primefree():
N = reduce(mul,Target)
for i in range(2,101):
nonP = N + i
judge = isprime(nonP)
print(nonP, judge)
return res
Target=[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47]
print(Target)
primefree()
[OUT]
#[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47]
614889782588491412 False
614889782588491413 False
614889782588491414 False
614889782588491415 False
614889782588491416 False
614889782588491417 False
614889782588491418 False
614889782588491419 False
614889782588491420 False
614889782588491421 False
614889782588491422 False
614889782588491423 False
614889782588491424 False
614889782588491425 False
614889782588491426 False
614889782588491427 False
614889782588491428 False
614889782588491429 False
614889782588491430 False
614889782588491431 False
614889782588491432 False
614889782588491433 False
614889782588491434 False
614889782588491435 False
614889782588491436 False
614889782588491437 False
614889782588491438 False
614889782588491439 False
614889782588491440 False
614889782588491441 False
614889782588491442 False
614889782588491443 False
614889782588491444 False
614889782588491445 False
614889782588491446 False
614889782588491447 False
614889782588491448 False
614889782588491449 False
614889782588491450 False
614889782588491451 False
614889782588491452 False
614889782588491453 False
614889782588491454 False
614889782588491455 False
614889782588491456 False
614889782588491457 False
614889782588491458 False
614889782588491459 False
614889782588491460 False
614889782588491461 False
614889782588491462 False
614889782588491463 False
614889782588491464 False
614889782588491465 False
614889782588491466 False
614889782588491467 False
614889782588491468 False
614889782588491469 False
614889782588491470 False
614889782588491471 False
614889782588491472 False
614889782588491473 False
614889782588491474 False
614889782588491475 False
614889782588491476 False
614889782588491477 False
614889782588491478 False
614889782588491479 False
614889782588491480 False
614889782588491481 False
614889782588491482 False
614889782588491483 False
614889782588491484 False
614889782588491485 False
614889782588491486 False
614889782588491487 False
614889782588491488 False
614889782588491489 False
614889782588491490 False
614889782588491491 False
614889782588491492 False
614889782588491493 False
614889782588491494 False
614889782588491495 False
614889782588491496 False
614889782588491497 False
614889782588491498 False
614889782588491499 False
614889782588491500 False
614889782588491501 False
614889782588491502 False
614889782588491503 False
614889782588491504 False
614889782588491505 False
614889782588491506 False
614889782588491507 False
614889782588491508 False
614889782588491509 False
614889782588491510 False
＜振り返り＞
これって、むりやり空白地帯を作るやり方でした。
では地道に、素数探検隊に実地調査してもらったらどうなるでしょうか。
平方数の発想で、1000個の隙間ができるためには多めに見積もって、
1000×1000＝1000000までの範囲で100個の素数空白地帯を探してみよう。
[IN]Python
##NonPrimeArray2
from sympy import primerange
def find_first(goal_length):
# 探索範囲の素数を少し大きめに範囲をとる
primes = list(primerange(2, 1000000))
for i in range(len(primes) - 1):
p1 = primes[i]
p2 = primes[i+1]
gap = p2 - p1
if gap >= goal_length:
return p1, p2, gap
p1, p2, gap = find_first(101)
print(f"素数探検隊の報告：素数{p1}と素数{p2}の間に、",end="")
print(f"{gap-1} 個の空白が発見されました！")
[OUT]
素数探検隊の報告：素数370261と素数370373の間に、111 個の空白が発見されました！
こんな身近に連続する100個以上の素数空白地帯があったんだね。
3。等差数列の中にも無限の素数
等差数列といえば、
初項と公差が同じである「倍数の数列」からみれば、
等差数列は初項と公差が一致するとはかぎらない。
言い換えると、倍数の数列をずらしたものだ。
たとえば、「３の倍数」の数列をエラトステネスのふるいと比べてみよう。
３，６，９，１２，１５、……………
エラトステネスのふるいでは「３だけ残して、あと数列の残党はバッサリ」捨てられている。
素数の倍数列では素数は先頭しかない。
合成数の倍数列ではどうだろう。「まるごと捨て」られる。
４，８，１２，１６，２０、…………
それにくらべて等差数列はどうだろう。
うまく、エラトステネスのふるいから「丸ごと捨てられるような無残なこと」にはならない。
それどころか、けっこう素数が出現する。ところどころ捨てられても生き残りが多数ありそうだ。
けっこうどころか無限に素数があるというのだ。
＜３で割った余りが２の数列＞
３で割って１余る数列をS1と書きしょう。
S1＝｛１，4，7，１0，１3，１6、１９，２２，２５，２８，３１，３４、………｝
素数は７，１３，１９，３１，…と出現しますね。
３で割って２余る数列をS2と書きしょう。
S2=｛２，５，８，１１，１４，１７、２０，２３，２６，２９，３２，３５、………｝
素数は２，５，１１，１７，２３，２９、…とわりと多く出現しますね。
合成数に目をつけましょう。
S1の４＝2^2、10＝２×5、16＝2^4、22＝2×11、25＝5^2、28＝2^2*7, 34=２×17,....
S2の８＝2^3、14＝2×7、20＝2^2*5, 26=2×13、32＝2^5, 35＝５×7
合同式で考えると、S2の要素=-1(mod 3)、S1の要素=1(mod 3)です。
１は何個かけても１と合同ですが、－１を作るにはー１を奇数個かけないといけません。
１をー１と１をまぜて作るにはー１を偶数個かけないといけません。
だから、
S１の合成数を見ると、S２の素数を偶数個かけてます。
S２の合成数を見ると、S２の素数を1個以上の奇数個かけてます。
さて、ここで、２からｐまでの素数の積をNとしましょう。
Nは３の倍数ですね。
だからM１＝N+1、M2=N-1とすると、M1はS1に登場し、M2はS2に登場する巨大数ですね。
ユークリッド、ピタゴラスの素数無限論で登場するのはM1型でしたね。
Ps＝[2, 3, 5, 7, 11, 13, 17, 19, 23]と最大をp=23としてN＝product(Ps)+1としたとき、
N=223092871は素数ではないが、最小の素因数317はPsの最大の素数23より大きかった。
ここに着目すれば、少数の素数リストからスタートして、素数を追加することはできた。
それが、マリンのアルゴリズムだった。
M1型の発想で作ったNはS1に属するけれど、Nの素因数がS1に続するかどうかは不明だね。
M２型ならどうでしょう。M2が素数ならｐより大きい素数があるということになります。
M２が合成数ならどうでしょうか。NはS2型だからNの素因数にS２型が1個はあるはずです。
つまり、S2には素数は無限に入るということだね。
課題：S２型の素数が無限に作られるコードにとりくんでみよう。
土台になる素数リストはS2型に限定しなくてもいいです。
新しく追加された素数をS2型としましょう。
#S2type-mullin
from functools import reduce
from operator import mul
from sympy import factorint
def s2_prime_generator():
Ps = [2, 3, 5] # S1,S2に関係なく土台として使う素数リスト
Start = Ps
S2 =[]
while True:
yield S2
# M2型: これまでの積 - 1
M2 = reduce(mul, Ps) - 1
factors = factorint(M2).keys()
M = min([f for f in factors if f % 3 == 2]) #S2型の最小素因数
Ps = Ps + [M]
S2 = [x for x in Ps if x not in Start]
# 探検隊、M2型のルートで8ステップ追跡
generator = s2_prime_generator()
for i in range(8):
current_list = next_generator = next(generator)
print(f"Step {i+1}: {current_list}")
[OUT]
Step 1: []
Step 2: [29]
Step 3: [29, 11]
Step 4: [29, 11, 1367]
Step 5: [29, 11, 1367, 13082189]
Step 6: [29, 11, 1367, 13082189, 89]
Step 7: [29, 11, 1367, 13082189, 89, 59]
Step 8: [29, 11, 1367, 13082189, 89, 59, 29819952677]
＜３で割った余りが１の数列＞
では、３で割った余りが１の数列、S1に無限の素数があるだろうか。
最初に書き出した、S1とS2の数列をみてみよう。
S1＝｛１，4，7，１0，１3，１6、１９，２２，２５，２８，３１，３４、………｝
S2=｛２，５，８，１１，１４，１７、２０，２３，２６，２９，３２，３５、………｝
素数がS1とS2をほぼ、交互にいったりきたりしているようだから、S２だけに素数が無限にあるのも
違う気がする。
ユークリッド、ピタゴラスの素数無限論で登場するM1型の実験を思い出そう。
Ps＝[2, 3, 5, 7, 11, 13, 17, 19, 23]と最大をp=23としてN＝product(Ps)+1としたとき、
N=223092871は素数ではないが、最小の素因数317はPsの最大の素数23より大きかった。
それで、317を素数リストPsに追加したね。
しかし、S1とS2に分けていなかったから、よかったけれど、317は3で割った余りは２だから使えない。
703763もよく見ると３で割った余りは２だ。223092871は素数でない。
ということは、Psを土台にして、S1の先を作ろうとしたのだが、いきなり破綻した。
そこで、Nをユークリッド式の1次式にするのをやめる。
整数を３で割った余りで分類しているから、群としては｛0，1，2｝という剰余類になる。
剰余群は巡回群になるから、複素平面の円周の3等分点と同型になる。だから、z^3=1となる。
ｚ^3-1=(z-1)(z^2+z+1)と分解できて、ｚ＝１では動きがない。
だから、z^2+z+1と同型にしたらどうだろうか。
たとえば、t=product(Ps)とするとき、N=ｔ＋１ではなく、N=t^2+t+1にするのだ。
検証しやすく土台はPs=[2,3,5]にしてみよう。
tは３の倍数だから、Nは素数かどうかは別として余り１だ。Nの素因数で余り１を追加していこう。
最初はｔ＝２×3×5＝30、N＝931。よさそうだね。
課題：３で割って１余る数列で無限に素数を作るコードにとりくもう。
[IN]Python
#S1type-mullin
from functools import reduce
from operator import mul
from sympy import factorint
def s1_prime_generator():
Ps = [2, 3, 5] # S1,S2に関係なく土台として使う素数リスト
Start = Ps
S1 =[]
while True:
yield S1
# 2次の生成式
t = reduce(mul, Ps)
N = t**2 + t + 1
factors = factorint(N).keys()
M = min([f for f in factors if f % 3 == 1]) #S1型の最小素因数
Ps = Ps + [M]
S1 = [x for x in Ps if x not in Start]
# 探検隊、S1型のルートで8ステップ追跡
generator = s1_prime_generator()
for i in range(10):
current_list = next_generator = next(generator)
print(f"Step {i+1}: {current_list}")
[OUT]
Step 1: []
Step 2: [7]
Step 3: [7, 73]
Step 4: [7, 73, 13]
Step 5: [7, 73, 13, 163]
Step 6: [7, 73, 13, 163, 43]
Step 7: [7, 73, 13, 163, 43, 1006867933]
Step 8: [7, 73, 13, 163, 43, 1006867933, 607]
Step 9: [7, 73, 13, 163, 43, 1006867933, 607, 1117]
Step 10: [7, 73, 13, 163, 43, 1006867933, 607, 1117, 909309821128128549785509674722967283825203785371]
よい感じで波打って出力できました。
＜振り返り＞
なぜ、N = t^2 + t + 1がS1型の素数ｐしか生まないのだろうか。
ある素数ｐがNを割り切るとしたら、N≡０（modp)とかけるね。
t^2 + t + 1≡０（modp)の両辺にt-1をかけると、t^3-1≡0(modp)となる。
これは、ｔの3乗がmod pでは1と合同だから、
ｔは3乗して1になる周期（位数3）を持つということだね。
フェルマーの小定理（ｔp-1 ≡1 (mod p)）から、p-1 が 3 の倍数でなければならないと言える。
これから素数pは３の倍数＋１なので、S1に追加できたわけだ。
課題：geogabraを使って、マリン数列が作れるでしょうか。
Step 9: [2, 3, 7, 43, 13, 53, 5, 6221671, 38709183810571]
iterationListがうまく動くかわからないのと、
もともと6221671のあとにけた数が爆発するので、手動でいいでしょう。
a=Element(divisors(product({2,3})+1),2)
b=Element(divisors(product({2,3,7})+1),2)
c=Element(divisors(product({2,3,7,43})+1),2)
d=Element(divisors(product({2,3,7,43,13})+1),2)
e=Element(divisors(product({2,3,7,43,13,53})+1),2)
f=Element(divisors(product({2,3,7,43,13,53,5})+1),2)