毎日コツコツ素数づくり

 投稿者:名和長泰メール  投稿日:2010年 5月19日(水)17時56分22秒
  2010年1月はじめから、十進BASIC(2進およそ16桁モード)で素数を作っています。
ファイルに書き出し、適当な大きさで止めて、新たなファイルに書き出しています。
PC(Pentium D 2.80GHz, 1.2GB RAM)を走らせ続け、
5ヶ月かかってようやく10億個めの素数をこえました。

ご参考まで、9桁の素数(10億+1個目からの100個)までの一部です。
 

Re: 毎日コツコツ素数づくり

 投稿者:山中和義  投稿日:2010年 5月21日(金)15時37分12秒
  > No.1237[元記事へ]

名和長泰さんへのお返事です。

> 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
 

戻る