|
|
> No.4186[元記事へ]
白石 和夫さんへのお返事です。
プログラムは、素数リストを読み込んで計算しています。
実験記録をまとめました。
実行プログラムの概要
計算開始 計算終了
K=300000 2016 11 20/17:55:42 から 2016 12 05/05:59:58
K=399998 Ver.2 2016 11 19/22:15:33 から 2016 11 22/11:55:55 エラー
k=500000 2016 11 20/06:51:14 から 2016 12 01/13:15:12
プログラム別の詳細
-----------------------------------------
K=300000
20161120/17:55:42
20161120/18:54:34出力
J1= 1 20161120/21:40:08
k_1 (出力ファイル)
18115519350289 , 614267441903
途中省略
18115520350288 , 614267474747
------------------------------
途中省略
------------------------------
J1= 136 20161205/05:30:28
k_137 (出力ファイル)
18115655350289 , 614271896866
途中省略
18115655550000 , 614271903479
計算終了
20161205/05:59:58
------------------------------
K=300000 計算時間
2016 11 20/17:55:42
2016 12 05/05:59:58
-----------------------------------------
k=500000
20161120/06:51:14 k=500000
20161120/09:52:50出力
10895.52秒
J1= 1 20161120/10:21:48
k_1 (出力ファイル)
54299021851369 , 1775034096424
途中省略
54299021851369 , 1775034096424
------------------------------
途中省略
------------------------------
J1= 589 20161201/13:01:43
k_500 (出力ファイル)
54299080751369 , 1775035959029
途中省略
54299080801680 , 1775035960614
計算終了
20161201/13:15:12
--------------------------------
k=500000 計算時間
2016 11 20/06:51:14
2016 12 01/13:15:12
-----------------------------------------
K=399998 Ver.2
k_1
33640568402401 , 1117217274291
k_159
33640584258593 , 1117217783657
計算開始
20161119/22:15:33
20161120/00:03:26出力
J1= 1 20161120/00:25:08
J1= 2 20161120/00:46:06
J1= 3 20161120/01:06:54
途中省略
J1= 156 20161122/11:10:26
J1= 157 20161122/11:33:14
J1= 158 20161122/11:55:55
エラー
-----------------------------------------
素数リスト 2～ 7368791 作成
! 素数リスト作成プログラム
OPTION ARITHMETIC NATIVE !2進モード
ASK DIRECTORY d$
LET f_name1$="prime_500001" ! ファイル名
LET f1$=d$&"\"&f_name1$&".txt"
PRINT "ファイル保存場所[1] = ";f1$
OPEN #1: NAME f1$
ERASE #1
LET t0=TIME
PRINT #1:2 ! 素数2をプリント
!エラトステネスの奇数列篩
LET k7=7368791
DIM P(k7)
LET h1=2
FOR n1=3 TO k7 STEP 2
IF P(n1)=0 THEN
PRINT #1:n1 ! 素数3からをプリント
LET h1=h1+1
END IF
FOR k1=n1 TO k7 STEP 2
LET m1=n1*k1
IF m1>k7 THEN GOTO 10
LET P(m1)=1
NEXT k1
10 NEXT n1
CLOSE #1
PRINT TIME-t0;"秒で計算しました"
END
-----------------------------------------
! K=300000 Ver.2.1
DECLARE EXTERNAL FUNCTION ISPRIME
DECLARE EXTERNAL FUNCTION POWMOD
OPTION ARITHMETIC RATIONAL
PRINT "K=300000"
PRINT DATE$;"/"; TIME$
LET k9=300000 !素数の個数
DIM A(k9+1)
OPEN #2: NAME "PRIME_300001.txt" ,ACCESS INPUT
FOR i=1 TO k9+1
INPUT #2: A(i)
NEXT i
CLOSE #2
LET X=1 !分子
LET Y=1 !分母
FOR i=1 TO K9
LET X=X*(A(i)-1) !1,2,4,6,10,12,16,18,…
LET Y=Y*A(i) !2,3,5,7,11,13,17,19,…
NEXT i
LET START=A(K9)^2 !4256233(300000th prime)
LET STOPT=A(K9+1)^2-1 !4256249(300001st prime)
LET prime=614267441903 !18115519350289 (614267441903rd prime)
LET SERIES=x/y
LET x1=NUMER(SERIES) !分子
LET y1=DENOM(SERIES) !分母
LET S=(prime-(START*x1/y1+k9-1))*y1 !DATA初項
LET g=y1-x1
OPEN #1:NAME "k_"&STR$(1)&".txt",RECTYPE INTERNAL
ERASE #1
LET J1=1
PRINT DATE$;"/"; TIME$;"出力"
FOR N=START+1 TO STOPT+1
IF cc>=1E6*J1 THEN
PRINT "J1=";J1;DATE$;"/"; TIME$
CLOSE #1
LET J1=J1+1
OPEN #1:NAME "k_"&STR$(J1)&".txt",RECTYPE INTERNAL
END IF
LET pn=((n-1)*x1+S)/y1+k9-1
WRITE #1:(n-1),pn
IF ISPRIME(N)=1 THEN
LET S=S+g
ELSE
LET S=S-x1
END IF
LET cc=cc+1
NEXT N
CLOSE #1
PRINT DATE$;"/"; TIME$
END
EXTERNAL FUNCTION ISPRIME(N)
OPTION ARITHMETIC RATIONAL
IF N=2 THEN
LET ISPRIME=1
EXIT FUNCTION
END IF
IF N = 1 OR MOD(N,2)=0 THEN
LET ISPRIME=0
EXIT FUNCTION
END IF
LET D=(N-1)/2
LET S=0
DO WHILE MOD(D,2)=0
LET D=INT(D/2)
LET S=S+1
LOOP
FOR I=1 TO 7
LET ISP=0
READ A !' n < 341550071728321 なら a = 2, 3, 5, 7, 11, 13, 17
DATA 2,3,5,7,11,13,17
LET ISP=0
LET R=POWMOD(A,D,N)
IF R=1 OR R=N-1 THEN
LET ISP=1
END IF
LET R=POWMOD(R,2,N)
FOR J=0 TO S-1
IF R=N-1 THEN
LET ISP=1
END IF
LET R=POWMOD(R,2,N)
NEXT J
IF ISP=0 THEN
LET ISPRIME=0
EXIT FUNCTION
END IF
NEXT I
LET ISPRIME=1
END FUNCTION
EXTERNAL FUNCTION POWMOD(B,P,M)
OPTION ARITHMETIC RATIONAL
LET RESULT=1
DO WHILE P>0
IF MOD(P,2)=1 THEN
LET RESULT=MOD(RESULT*B,M)
END IF
LET B=MOD(B*B,M)
LET P=INT(P/2)
LOOP
LET POWMOD=RESULT
END FUNCTION
------------------------------------
! K=500000
! k=x s=係数(DATA) Ver 3.1
! 素数の個数 π(x)=(n*x+S)/y+k-1
DECLARE EXTERNAL FUNCTION ISPRIME
DECLARE EXTERNAL FUNCTION POWMOD
OPTION ARITHMETIC RATIONAL
PRINT DATE$;"/"; TIME$;" ";"k=500000"
LET t0=TIME
LET k9=5E5 !素数の個数(1から)
DIM A(k9+1)
OPEN #2: NAME "c:\prime_1E8.txt" ,ACCESS INPUT
FOR i=1 TO k9+1
INPUT #2: A(i)
NEXT i
CLOSE #2
20 LET X=1 !分子
LET Y=1 !分母
FOR i=1 TO K9
LET X=X*(A(i)-1) !1,2,4,6,10,12,16,18,…
LET Y=Y*A(i) !2,3,5,7,11,13,17,19,…
! PRINT i;x/y
NEXT i
LET START=A(K9)^2 !7368787(500000th prime)
LET STOPT=A(K9+1)^2-1 !7368791(500001st prime)
LET prime=1775034096424 !54299021851369 (1775034096424)
LET SERIES=x/y
LET x1=NUMER(SERIES) !分子
LET y1=DENOM(SERIES) !分母
LET S=(prime-(START*x1/y1+k9-1))*y1 !DATA初項
OPEN #1:NAME "k_"&STR$(1)&".txt",RECTYPE INTERNAL
ERASE #1
LET J1=1
PRINT DATE$;"/"; TIME$;"出力"
LET TM=TIME-t0
PRINT USING"#########." & REPEAT$("#",2):TM;
PRINT "秒"
LET N=START
DO WHILE N < STOPT+1
IF cc>=1E5*J1 THEN
PRINT "J1=";J1;DATE$;"/"; TIME$
CLOSE #1
LET J1=J1+1
OPEN #1:NAME "k_"&STR$(J1)&".txt",RECTYPE INTERNAL
END IF
LET pn=(n*x1+S)/y1+k9-1
WRITE #1:n,pn
LET S=S-x1
IF ISPRIME(N+1)=1 THEN LET S=S+y1
LET cc=cc+1
LET N=N+1
LOOP
CLOSE #1
PRINT DATE$;"/"; TIME$
END
EXTERNAL FUNCTION ISPRIME(N)
OPTION ARITHMETIC RATIONAL
IF N=2 OR N=3 OR N=5 OR N=7 OR N=11 OR N=13 OR N=17 THEN
LET ISPRIME=1
EXIT FUNCTION
END IF
IF N=1 OR MOD(N,2)=0 THEN
LET ISPRIME=0
EXIT FUNCTION
END IF
LET D=(N-1)/2
LET S=0
DO WHILE MOD(D, 2)=0
LET D=INT(D/2)
LET S=S+1
LOOP
FOR I=1 TO 7
LET ISP=0
READ A !' n < 341550071728321 なら a=2,3,5,7,11,13,17
DATA 2,3,5,7,11,13,17
LET ISP=0
LET R=POWMOD(A,D,N)
IF R=1 OR R=N-1 THEN
LET ISP=1
END IF
LET R=POWMOD(R,2,N)
FOR J=0 TO S-1
IF R=N-1 THEN
LET ISP=1
END IF
LET R=POWMOD(R,2,N)
NEXT J
IF ISP=0 THEN
LET ISPRIME=0
EXIT FUNCTION
END IF
NEXT I
LET ISPRIME=1
END FUNCTION
EXTERNAL FUNCTION POWMOD(B,P,M)
OPTION ARITHMETIC RATIONAL
LET RESULT=1
DO WHILE P>0
IF MOD(P,2)=1 THEN
LET RESULT=MOD(RESULT*B,M)
END IF
LET B=MOD(B*B,M)
LET P=INT(P/2)
LOOP
LET POWMOD=RESULT
END FUNCTION
http://blogs.yahoo.co.jp/donald_stinger
|
|
