|
|
> No.4301[元記事へ]
しばっちさんへのお返事です。
素晴らしいプログラムをありがとうございました。
動作報告
prime quadruples (30n+11, 30n+13, 30n+17, 30n+19).
30n+11=1000000000067441 TO 1000000099262441
299 個
1000兆までは、まだ遠いです。下記のプログラムを頼りにしています。
エラトステネスの篩 #4264
素数判定 Miller-Rabin法 #3894
エラトステネスの篩 eratos_dll #3985
Re: エラトステネスの篩 #4101
460000000000000 : 14054559347931 : 0
470000000000000 : 14350653588527 : 2958176
計算時間：6日と4時間33分46秒
-----------------------------------------------
下記は数値実験用プログラムです。
A014561 Numbers n giving rise to prime quadruples (30n+11, 30n+13, 30n+17, 30n+19).
上記の数列を生成するプログラム
prime quadruples BASIC Acc program
-------------------------------
!30n± k篩
OPTION ARITHMETIC NATIVE
LET t0=TIME
LET k6=31627
LET k2=3402
!エラトステネスの篩 !#4264
LET Fu=k6
DIM P(Fu)
DIM A(k2) !素数
MAT P=ZER
MAT A=ZER
LET A(1)=2
LET H1=1
FOR I=3 TO SQR(Fu) STEP 2
IF P(I)=0 THEN
FOR J=I*I TO Fu STEP I
LET P(J)=1
NEXT J
END IF
NEXT I
FOR I=3 TO Fu STEP 2
IF P(I)=0 THEN
LET H1=H1+1
LET A(H1)=I
END IF
NEXT I
DIM B(4)
DATA 11,13,17,19
MAT READ B
LET Q=30
LET ka=INT(1E9/Q)+1
LET kb=1E9
LET kc=k2
LET kd=INT(kb/207)
DIM D(ka)
MAT D = ZER
LET cj=2 !(5, 7, 11, 13, 17), (11, 13, 17, 19, 23),
FOR r=1 TO 4
LET rr=B(r)
FOR t=4 TO kc
LET x=A(t)
IF MOD(x+rr,Q)=0 THEN
LET y=-(x+rr)/Q
GOTO 80
END IF
IF MOD(x-rr,Q)=0 THEN
LET y=(x-rr)/Q
GOTO 80
END IF
FOR i=2 TO x !*.86
IF MOD(Q*i-rr,x)=0 THEN
LET y=-i
EXIT FOR
END IF
NEXT i
80 FOR f=1 TO kd
IF x*f+y>ka THEN EXIT FOR
LET D(x*f+y)=1
NEXT f
NEXT t
NEXT r
OPEN #1:NAME "E:\pi4_1E_9.txt",RECTYPE INTERNAL
ERASE #1
WRITE #1:0
FOR n=1 TO ka
!IF n*Q+rr>kb THEN EXIT FOR
IF D(n)=0 THEN
LET cj=cj+1
WRITE #1:n
END IF
NEXT n
CLOSE #1
PRINT cj
LET TM=TIME-t0
PRINT USING"####." & REPEAT$("#",2):TM;
PRINT "秒"
END
-----------------------
計算結果
28388 37.07秒
30n+k 篩 の探究の結果 このアルゴリズムにたどり着きました。
精度確認中です。
FOR i=2 TO x !*.86
IF MOD(Q*i-rr,x)=0 THEN
LET y=-i
EXIT FOR
END IF
NEXT i
素数の係数を求めるループがネックで1000億位が限度です。
それで
6n±1 篩
十進BASIC 1000桁モードで 16桁も行ってみた。
!6n±1篩 Ver.1
!四つ子素数計数関数（π4(x) Prime Quadruplet-counting function）
DECLARE EXTERNAL FUNCTION cprime
DECLARE EXTERNAL SUB sqprime
OPTION ARITHMETIC DECIMAL_HIGH
PRINT DATE$;"/"; TIME$
LET t0=TIME
LET k=31716427
LET k2=1957369
DIM GT(k2)
OPEN #3:NAME "E:\prime_1E8xv.txt",ACCESS INPUT
FOR i=1 TO k2
INPUT #3: GT(i)
NEXT i
CLOSE #3
! MAX 1005931741646329=31716427^2
LET SA=100000001 !MIN 2 (1E+7 to)
LET SB=100001000 !MAX 1E15-1
CALL sqprime(SA,SB,k2,GT)
DIM Bb(SA TO SB)
OPEN #2:NAME "prime_pi1.txt",ACCESS INPUT
FOR i=SA TO SB
INPUT #2: BB(i)
NEXT i
CLOSE #2
! !SA=数値入力+1 (1E7)
LET c1=0 !1E7 (1E7) to 2E7
PRINT SA-1;":";c2;":";c1;":";c
LET k6=SA*1E7-1E7
FOR n=SA TO SB
LET t1=TIME
LET k2=Bb(n) !π(x)
LET k6=k6+1E7
LET C=cprime(k2,k6,GT)
LET c1=c1+c
PRINT n*1E7;":";c2+c1;":";c1;":";c;
LET TM=TIME-t1
PRINT USING"######." & REPEAT$("#",2):TM;
PRINT "秒"
NEXT n
LET TM=TIME-t0
PRINT USING"######." & REPEAT$("#",2):TM;
PRINT "秒"
PRINT DATE$;"/"; TIME$
END
EXTERNAL FUNCTION cprime (k2,k6,GT())
OPTION ARITHMETIC DECIMAL_HIGH
LET k4=k6-1E7
LET B6=INT(k6/30)
LET U=INT(k6/6)
LET W=INT(k4/6)
LET M7=W
DIM x(W-M7 TO U-M7)
MAT x = ZER !(6*n-1)
FOR n=3 TO k2
LET P6=GT(n)
LET G1=INT(W/P6)!-2
IF MOD(P6+1,6)=0 THEN !(6*n-1)
LET r6=(P6+1)/6
FOR i=G1 TO B6
IF P6*i+r6<W THEN GOTO 140
IF P6*i+r6>U THEN EXIT FOR
LET x(P6*i+r6-M7)=1
140 NEXT i
END IF
IF MOD(P6-1,6)=0 THEN
LET r6=(P6-1)/6
FOR i=G1 TO B6
IF P6*i-r6<W THEN GOTO 150
IF P6*i-r6>U THEN EXIT FOR
LET x(P6*i-r6-M7)=1
150 NEXT i
END IF
IF MOD(P6+1,6)=0 THEN !(6*n+1)
LET r6=(P6+1)/6
FOR i=G1 TO B6
IF P6*i-r6<W THEN GOTO 160
IF P6*i-r6>U THEN EXIT FOR
LET x(P6*i-r6-M7)=1
160 NEXT i
END IF
IF MOD(P6-1,6)=0 THEN
LET r6=(P6-1)/6
FOR i=G1 TO B6
IF P6*i+r6<W THEN GOTO 170
IF P6*i+r6>U THEN EXIT FOR
LET x(P6*i+r6-M7)=1
170 NEXT i
END IF
NEXT n
LET Cc=0
FOR n=W-M7 TO U-M7-1
LET ST=n+M7
IF x(n)=0 AND x(n+1)=0 AND 6*ST-1>k4 AND 6*ST-1<k6 THEN
LET cc=cc+1
END IF
NEXT n
LET cprime=cc
END FUNCTION
EXTERNAL SUB sqprime(SA,SB,k2,GT())
OPTION ARITHMETIC DECIMAL_HIGH
OPEN #2:NAME "prime_pi1.txt",RECTYPE INTERNAL
ERASE #2
LET x=1E7
FOR i=1 TO k2-1
LET v=GT(i)^2
LET vi=GT(i+1)^2
LET vv=vi-v
IF x =< v THEN
CALL sqrp(1)
IF vv>=1E7 THEN
LET vt=INT(vv/1E7)
CALL sqrp(vt)
END IF
END IF
NEXT i
SUB sqrp(vo)
FOR ns=1 TO vo
LET C=C+1
LET X=X+1E7
IF C>=SA AND SB>=C THEN
WRITE #2:i
END IF
NEXT ns
END SUB
CLOSE #2
END SUB
----------------------------------------------
双子素数-双子素数=6 配列にするとこのような状況
管理用メモリー限界のため一千万ループこれがネックです。
計算結果
100000000 : 3314576487 : 0 : 0
1000000010000000 : 3314576511 : 24 : 24
1000000020000000 : 3314576542 : 55 : 31
1000000030000000 : 3314576565 : 78 : 23
1000000040000000 : 3314576601 : 114 : 36
1000000050000000 : 3314576639 : 152 : 38
1000000060000000 : 3314576675 : 188 : 36
1000000070000000 : 3314576697 : 210 : 22
1000000080000000 : 3314576730 : 243 : 33
1000000090000000 : 3314576758 : 271 : 28
1000000100000000 : 3314576786 : 299 : 28
素数リストを読み込むので、遅いです。
計算時間 1千万あたり30秒*10 約300秒
Re: エラトステネスの篩 #4101
これでは、比較にならないです。
C++ のマルチプラットフォームも出来たらいいな～。
http://blogs.yahoo.co.jp/donald_stinger
|
|