|
|
『オイラーの無限解析』レオンハルト・オイラー著 海鳴社 参照
----------------------------------------------------------
!4n±1篩 で円周率を求める。
OPTION ARITHMETIC NATIVE !2進モード
LET t0=TIME
LET k=11839
LET k2=1421
DIM P(k)
DIM A(k2) !素数
SUB prime(v) !エラトステネスの奇数列篩
LET k9=v
LET h1=1
LET A(h1)=2
LET h1=2
FOR n1=3 TO k9 STEP 2
IF P(n1)=0 THEN
LET A(h1)=n1
LET h1=h1+1
IF h1>k9+1 THEN GOTO 20
END IF
FOR k1=n1 TO k9 STEP 2
LET m1=n1*k1
IF m1>k9 THEN GOTO 10
LET P(m1)=1
NEXT k1
10 NEXT n1
20
END SUB
CALL prime(k)
LET fnpi=1-1/2
LET Q=4
LET k6=1.4E8 !篩の計算範囲
LET k5=k6/Q+1
DIM Au(k5),A7(k6)
MAT Au = ZER !(4*n-1)
FOR n=1 TO k2-1
LET Pu=A(n+1)
IF MOD(Pu+1,Q)=0 THEN
LET ru=(Pu+1)/Q
FOR i=1 TO k5
IF Pu*i+ru>k5 THEN EXIT FOR
LET Au(Pu*i+ru)=1
NEXT i
END IF
IF MOD(Pu-1,Q)=0 THEN
LET ru=(Pu-1)/Q
FOR i=1 TO k5
IF Pu*i-ru>k5 THEN EXIT FOR
LET Au(Pu*i-ru)=1
NEXT i
END IF
NEXT n
FOR n=1 TO k5
IF Au(n)=1 THEN GOTO 100
LET fnpi=fnpi*(1+1/(Q*n-1))
100 NEXT n
MAT Au = ZER !(4*n+1)
FOR n=1 TO k2-1
LET Pu=A(n+1)
IF MOD(Pu+1,Q)=0 THEN
LET ru=(Pu+1)/Q
FOR i=1 TO k5
IF Pu*i-ru>k5 THEN EXIT FOR
LET Au(Pu*i-ru)=1
NEXT i
END IF
IF MOD(Pu-1,Q)=0 THEN
LET ru=(Pu-1)/Q
FOR i=1 TO k5
IF Pu*i+ru>k5 THEN EXIT FOR
LET Au(Pu*i+ru)=1
NEXT i
END IF
NEXT n
FOR n=1 TO k5
IF Au(n)=1 THEN GOTO 200
LET fnpi=fnpi*(1-1/(Q*n+1))
200 NEXT n
PRINT 1/fnpi
PRINT PI/2
LET TM=TIME-t0
PRINT USING"####." & REPEAT$("#",2):TM;
PRINT "秒"
END
---------------------------------------------
計算結果
BASIC Acc 1億5千万までの素数
fnpi= 1.57080705636575
pi/2= 1.5707963267949
1億辺りから1.570変わらずです。
---------------------------------------------
素数生成プログラム
!4n±1篩
OPTION ARITHMETIC NATIVE !2進モード
LET t0=TIME
LET k=10039
LET k2=1233
DIM P(k)
DIM A(k2) !素数
SUB prime(v) !エラトステネスの奇数列篩
LET k9=v
LET h1=1
LET A(h1)=2
LET h1=2
FOR n1=3 TO k9 STEP 2
IF P(n1)=0 THEN
LET A(h1)=n1
LET h1=h1+1
IF h1>k9+1 THEN GOTO 20
END IF
FOR k1=n1 TO k9 STEP 2
LET m1=n1*k1
IF m1>k9 THEN GOTO 10
LET P(m1)=1
NEXT k1
10 NEXT n1
20
END SUB
CALL prime(k)
LET Q=4
LET k6=1E8 !篩の計算範囲
LET k5=k6/Q+1
DIM Au(k5),A7(k6)
MAT Au = ZER !(4*n-1)
FOR n=1 TO k2-4
LET Pu=A(n+1)
IF MOD(Pu+1,Q)=0 THEN
LET ru=(Pu+1)/Q
FOR i=1 TO k5
IF Pu*i+ru>k5 THEN EXIT FOR
LET Au(Pu*i+ru)=1
NEXT i
END IF
IF MOD(Pu-1,Q)=0 THEN
LET ru=(Pu-1)/Q
FOR i=1 TO k5
IF Pu*i-ru>k5 THEN EXIT FOR
LET Au(Pu*i-ru)=1
NEXT i
END IF
NEXT n
FOR n=1 TO k5
IF Au(n)=1 THEN GOTO 100
LET A7(Q*n-1)=1
100 NEXT n
MAT Au = ZER !(4*n+1)
FOR n=1 TO k2-4
LET Pu=A(n+1)
IF MOD(Pu+1,Q)=0 THEN
LET ru=(Pu+1)/Q
FOR i=1 TO k5
IF Pu*i-ru>k5 THEN EXIT FOR
LET Au(Pu*i-ru)=1
NEXT i
END IF
IF MOD(Pu-1,Q)=0 THEN
LET ru=(Pu-1)/Q
FOR i=1 TO k5
IF Pu*i+ru>k5 THEN EXIT FOR
LET Au(Pu*i+ru)=1
NEXT i
END IF
NEXT n
FOR n=1 TO k5
IF Au(n)=1 THEN GOTO 200
LET A7(Q*n+1)=1
200 NEXT n
OPEN #1: TextWindow1
ERASE #1
PRINT #1:STR$(2)&",";
FOR n=3 TO k6 STEP 2
IF A7(n)=1 THEN
LET c1=c1+1
PRINT #1:STR$(n)&",";
END IF
NEXT n
CLOSE #1
PRINT c1+1
LET TM=TIME-t0
PRINT USING"####." & REPEAT$("#",2):TM;
PRINT "秒"
END
-------------------------------------------------
素数→円周率→ゼータ関数
『オイラーの無限解析』レオンハルト・オイラー著 海鳴社
やっと理解出来ました。十進BASICがなかったら生涯無理でした。
パソコンが無い時代にオイラーは、凄いと思いました。
素数周期の最小は4n±1と言えそうですが、双子素数が無限に続く予想
には、無力です。4n±1と6n±1の合わせ業とか? 模索してます。
メルセンヌ数が無限に続く予想には?
http://blogs.yahoo.co.jp/donald_stinger
|
|