素数と戯れる いやもつれる？ いやいや組んず解れつ？
発端はこちら。
このクイズに答えようとしたわけである。
パターンはわかった。
だけど、そんな巨大な数値を計算できるのか？
というわけで試してみたんである。
Python Wiki 20 lines: Prime numbers
まずは安易にこちらのPythonコードを利用した。
「20 lines: Prime numbers sieve w/fancy generators」のコードである。このコード自体は1000未満の素数を画面に出力する。
いちいちリンク先に飛ぶのも面倒だろうから、ソースコードをちょいと転記する。
import itertools
def iter_primes():
# an iterator of all numbers between 2 and +infinity
numbers = itertools.count(2)
# generate primes forever
while True:
# get the first number from the iterator (always a prime)
prime = next(numbers)
yield prime
# this code iteratively builds up a chain of
# filters...slightly tricky, but ponder it a bit
numbers = filter(prime.__rmod__, numbers)
for p in iter_primes():
if p > 1000:
break
print (p)「yield」などを使って、Python独特のようなコードである。このプログラムについてはこちらで解説さした。
これだと3桁の素数しか得られない。
なので、どーんと！
1000 を
100000000000 に
変更してみた。ゼロを8個追加。
すなわち1億倍！（笑）。
それをダラダラと標準出力にアウトプットされても困るのでファイルにリダイレクトした。
そうすると…………。
2時間45分くらいで強制終了した……orz。
ちなみに、できたファイルの最終行は次の通りである。
130013行目：1727101
全然、桁数足りない！
「filter」って、なんだか重そうだもんね（ただの偏見）。
検出した最大の素数は1727101
検出した素数の数は130013個
2〜1727101の自然数に
素数が130013個ある
130013÷1727101=0.07527816844527332
約7.5%が素数か
11桁の自然数であれば
100,000,000,000×0.075
=7,500,000,000
=約7G個
げ😱
1個の素数を64bit（8byte）で表現すると
7,500,000,000*8[byte]=60,000,000,000[byte]
=約55Gbyte
いや、あかん、あかん。
うちのスマホの空き容量足りまへん。いや、足りんではないけど、7G個の素数も魅力的ではあるけど、55Gは、ちょっと…………。そもそも、それを文字列にしたら更に倍！
う〜ん。
素数連鎖
と名付けていいのかどうか。
素数を全部リストアップするつもりやったけど、素数が予想以上に膨大で、スマホのストレージを圧迫しかねない。そうなると仕方がないよね。冒頭のクイズの条件に合う素数だけをピックアップするしかない。名付けて
素数連鎖
そんなネーミングでええのか？
コードはガラッと変更した。
import time
def is_prime(n):
if (n == 2):
return True;
if ((n & 0x01) == 0):
return False;
for i in range(3, n//2):
if (n % i) == 0:
return False;
return True;
def prime_n(min_n, max_n):
pn_next = min_n
pn = 0
for i in range(2, max_n):
isp = is_prime(i)
if isp == True:
pn = pn + 1
#print(pn, i)
if (pn == pn_next):
print(pn, i)
pn_next = i
print(time.localtime())
prime_n(6, 100000000000)
print(time.localtime())2，3，4，5，6，…………
と順に自然数を素数かどうかチェックしていく。
素数だったら「pn」をインクリメントする。
「pn」により、何番目の素数かがわかる。
「pn_next」番目の素数をピックアップする。
「pn_next」番目が見つかったら、その素数が次の「pn_next」になる。
実行してみたんだが。
time.struct_time(tm_year=2024, tm_mon=10, tm_mday=31, tm_hour=10, tm_min=9, tm_sec=42, tm_wday=3, tm_yday=305, tm_isdst=0)
6 13
13 41
41 179
179 1063
1063 8527
8527 87803待てど暮せど終わらないorz。
記録を取らなかったんだけど、たしか仕事中ずっと実行し続けてバッテリーだけをひたすら消費して帰宅するまでバッテリーが持つのかどうかが危ぶまれてきてCtrl-Cを押したんだったか。
エラトステネスの篩
素数の抽出って、もしかして意外にも最もベーシックなこの方法がいいんじゃね？
と、なんの根拠もなく自らの直感だけを頼りに第三の方法に打って出た（というほどに大袈裟なことでもないけど）。
もう、Pythonはいい。
高速だのみであれば君だ、君。
C言語
ソースコードはこれ。
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
void prime_chain(
const unsigned long long n_start,
const unsigned long long n_end,
unsigned long long*const p_next,
unsigned long long*const p_number
)
{
char*const prime_map = malloc(n_end - n_start + 1);
memset(prime_map, 1, n_end - n_start + 1);
unsigned long long n = 0;
unsigned long long n_half = n_end >> 2;
for (n = 2; n <= n_half; n++)
{
unsigned long long m = n_start / n;
if (m < 2) { m = 2; }
for (;;)
{
const unsigned long long nm = n * m;
if (nm < n_start) { m++; continue; }
if (n_end < nm) { break; }
// composite number
prime_map[nm - n_start] = 0;
m++;
}
}
#if 0
// show prime map
for (n = 0; n < 0x100; n++)
{
printf("%d", prime_map[n]);
if ((n & 0x0f) == 0x0f) {printf("\n");}
}
#endif
for (n = n_start; n < n_end; n++)
{
if (prime_map[n - n_start] != 0)
{
(*p_number)++;
if (*p_number == *p_next)
{
time_t now = time(NULL);
printf("%ld %llu %llu\n", now, *p_number, n);
*p_next = n;
}
}
}
free(prime_map);
return;
}
int main()
{
unsigned long long p_next = 6;
unsigned long long p_number = 0;
prime_chain(2, 1000000, &p_next, &p_number);
for (unsigned long long start = 1000001; start < 100000000000; start += 1000000)
{
if ((start % 1000000000) == 1)
{
time_t now = time(NULL);
printf("%ld start = %llu p_number = %llu\n", now, start, p_number);
}
prime_chain(start, start + 999999, &p_next, &p_number);
}
return 0;
}えーっと、いっちゃん最初のループでエラトステネスの篩を実行する。
２つ目の「#if 0 〜 #endif」まではデバッグ用。エラトステネスの篩の結果を表示している。
３つ目のループで素数連鎖を作り出して出力している。この関数「prime_chain」はちょっと汚いんやけど、気が向いたらそのうち綺麗にする。
ただ、関数「prime_chain」で一気に11桁の素数を対象にするとまたアプリがとんでいくので、ちょっとずつ実行することにする。それが「main」関数。1回の関数「prime_chain」呼び出しで実行するのは６桁だけ。これをあと5桁分繰り返すことになる。
結果は…………。
1729908585 6 13
1729908585 13 41
1729908585 41 179
1729908585 179 1063
1729908585 1063 8527
1729908585 8527 87803
1729908585 87803 1128889
1729908586 1128889 17624813
1729908668 17624813 326851121
1729909980 start = 1000000001 p_order = 50847534
1729912713 start = 2000000001 p_order = 98222287
1729917997 start = 3000000001 p_order = 144449537
1729923270 start = 4000000001 p_order = 189961812
1729933029 start = 5000000001 p_order = 234954223
1729955507 start = 6000000001 p_order = 279545368
1729964406 start = 7000000001 p_order = 323804352
1729965037 326851121 7069067389
1729974963 start = 8000000001 p_order = 367783654
1730010551 start = 9000000001 p_order = 411523195
1730038415 start = 10000000001 p_order = 455052511
1730064994 start = 11000000001 p_order = 498388617おしい！（のか？）。
あと一歩やねんけど。
でも、この「あと一歩」が凄まじい。
一番左の列は時間をナノ秒で表示している。
最後のナノ秒から、最初のナノ秒を差し引くと…。
1730064994[ns] − 1729908585[ns]
=156,409[sec]
=43.446944444444[Hr]
😱43時間！
まぁ、週末に実行してることも忘れて放置してたから(~_~;)
43時間計算し続けてみたけど、498,388,617番目までしか計算できなかった。そういうことになる。
9桁まで…………。
326,851,121番目の素数：7,069,067,389
自然数10,000,000,001〜11,000,000,001（10億個の自然数）の素数判定に7時間を超えていて、このペースで計算し続けたとして、11桁の素数を列挙し尽くすのにどれだけかかるのかというと…………。
1120時間 = 約46日
4、46日！🤯
こんなにもかかってしまふのか(~_~;)。
自然数の数値が大きくなると更に時間がかかるだろうしなぁ。
ああ、ぁぁ、…………。