名和長泰さんへのお返事です。
> 2010年1月はじめから、十進BASIC(2進およそ16桁モード)で素数を作っています。
> ファイルに書き出し、適当な大きさで止めて、新たなファイルに書き出しています。
私も作ってみました。
概要
5000個の素数を予め求めておく。(試行割算法の実用範囲)
この素数を使って、√nまでの篩いを行う。
「篩い」表の全体をメモリに実装できないので、分割して部分的に定義する。
1~約23億までの範囲の素数を求めることができる。
WindowsMe、PentiumⅢ700MHz、192MB 2進モードにて計測すると、
10万個なら 24 秒、100万個なら 315 秒(5.2分)となる。
※手続き PrimeList の場合、100万個なら20分程度となる。
ファイル出力なし、n番目の素数のみ表示では、100万個なら 17 秒、1千万個なら 305 秒(5分)となる。
出力を1つのテキストファイルにすると、10万個なら1.62MB、100万個なら18.2MBとなる。
1億個なら容量は2GB程度が予想される。いくつかに分割する必要もある。
!nまでの素数
LET t0=TIME
LET CP=5000 !※範囲を広げる場合は増やす 10,000で100億
DIM P(CP) !「篩い」の素数列
CALL PrimeList(CP,P) !2,3,5,7,11,13,17,19,23,…
!!!MAT PRINT P; !debug
LET N=P(CP)^2
PRINT USING "1 ~ ###,###,###,### までが対象範囲です。": N
!OPEN #1: NAME "TEST.TXT" !素数表
!!!OPEN #1: TextWindow1 !debug
!ERASE #1
!FOR i=1 TO CP !1~√n
!PRINT #1: i; P(i)
!NEXT i
LET SZ=4096
DIM F(SZ) !「篩い」表(奇数列 1,3,5,7,9,11,13,15,17,19,21,…) 0:素数でない、1:素数
LET c=CP !!100000001 !100,000,001番目 2,038,074,751
LET M=P(c) !!2038074751
LET top=M+1 !1~Nを範囲[top,btm]に分割する
IF MOD(top,2)=0 THEN LET top=top+1
IF top<3 THEN
LET c=1 !2は素数
!PRINT #1: 1; 2
END IF
DO
LET btm=top+2*(SZ-1)
!!!PRINT top; btm !debug
MAT F=CON !エラトステネスの篩い
IF top=1 THEN LET F(1)=0 !1を除く
FOR x=2 TO CP !3~√nの素数で
LET t=P(x)
IF t*t>btm THEN EXIT FOR
LET tt=MAX( INT((top-1)/t)+1, t )
IF MOD(tt,2)=0 THEN LET tt=tt+1
FOR k=t*tt TO btm STEP 2*t !2乗以上の倍数を削除する
LET F(MOD((k-M)/2-1,SZ)+1)=0
NEXT k
NEXT x
FOR k=1 TO SZ !この範囲での素数を表示する
IF F(k)=1 THEN
LET c=c+1 !c番目の素数
!PRINT #1: c; (top-1)+2*k-1
!!!PRINT c; (top-1)+2*k-1 !debug
IF c=1000001 THEN !指定の個数なら終了!
!!IF c=100000100 THEN !指定の個数なら終了!
PRINT c;"番目"; (top-1)+2*k-1 !!!
!CLOSE #1
PRINT "計算時間=";TIME-t0
STOP
END IF
END IF
NEXT k
LET top=top+2*SZ !次へ
LOOP
END
!試行割算法
EXTERNAL FUNCTION PrimeQ(n) !素数判定 1:素数、0:素数でない
LET PrimeQ=0
IF n<2 OR n<>INT(n) THEN EXIT FUNCTION !引数を確認する
!2以上の自然数なら
IF MOD(n,2)=0 THEN !2の倍数
IF n=2 THEN LET PrimeQ=1 !2は素数
ELSEIF MOD(n,3)=0 THEN !3の倍数
IF n=3 THEN LET PrimeQ=1 !3は素数
ELSE
LET k=5
DO WHILE k*k<=n !√nまで検証する
IF MOD(n,k)=0 THEN !5,11,17,23,29,…
EXIT FUNCTION
ELSEIF MOD(n,k+2)=0 THEN !7,13,19,25,31,…
EXIT FUNCTION
END IF !+1,+3,+5は2の倍数(偶数)、+1,+4は3の倍数、+5は5の倍数
LET k=k+6
LOOP
LET PrimeQ=1 !最後まで割り切れなければ、素数である
END IF
END FUNCTION
EXTERNAL SUB PrimeList(n,p()) !n個の素数列を返す
IF n<1 THEN EXIT SUB !引数を確認する
LET c=1 !見つけた個数
LET p(c)=2 !2は素数
IF n=1 THEN EXIT SUB
LET c=2
LET p(c)=3
IF n=2 THEN EXIT SUB
LET k=5
DO
IF PrimeQ(k)<>0 THEN !5,11,17,23,29,… が素数なら
LET c=c+1
LET p(c)=k
IF c=n THEN EXIT DO
END IF
IF PrimeQ(k+2)<>0 THEN !7,13,19,25,31,…
LET c=c+1
LET p(c)=k+2
IF c=n THEN EXIT DO
END IF
!+1,+3,+5は2の倍数(偶数)、+1,+4は3の倍数、+5は5の倍数
LET k=k+6 !次へ
LOOP
END SUB