|
|
素数階乗の探究
ウィルソンの定理を元にしています。
!素数階乗の探究
OPTION ARITHMETIC RATIONAL !有理数モード
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)
OPEN #1: TextWindow1
ERASE #1
PRINT #1:STR$(2)&",";
PRINT #1:STR$(3)&",";
LET cc=2
FOR r=1 TO 100
LET s=A(r)^2
LET t=A(r+1)^2-1
LET x=1
FOR i=1 TO r
LET x=x*A(i)
NEXT i
FOR n=s TO t
LET p1=n/x
IF DENOM(p1)=x THEN
LET cc=cc+1
PRINT #1:STR$(n)&",";
END IF
NEXT n
NEXT r
CLOSE #1
PRINT n;":";cc
LET TM=TIME-t0
PRINT USING"####." & REPEAT$("#",2):TM;
PRINT "秒"
END
--------------------------------------------
ウィルソンの定理は当たり前だと言う方にはつまらないプログラムですが
、このプログラムからスタートしてます。
!k=4 素数 2*3*5*7=210
OPTION ARITHMETIC RATIONAL !有理数モード
FOR n=7 TO 11^2-1
PRINT "N=";STR$(n)&" :";n/210
LET P=N*8/35
LET p1=+MOD(N,2)/2+MOD(N,3)/3+MOD(N,5)/5+MOD(N,7)/7
LET p2=-MOD(N,6)/6-MOD(N,10)/10-MOD(N,15)/15-MOD(N,14)/14-MOD(N,21)/21-MOD(N,35)/35
LET p3=+MOD(N,30)/30+MOD(N,105)/105+MOD(N,42)/42+MOD(N,70)/70
LET p4=-MOD(N,210)/210
LET Pn=P+p1+p2+p3+p4
PRINT "P=";STR$(pn)&"+3"
PRINT n;":";p1+p3+p
PRINT n;":";p2+p4
PRINT
NEXT n
END
-------------------------------------------
!素数階乗の素数篩 MOD
OPTION ARITHMETIC RATIONAL !有理数モード
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)
OPEN #1: TextWindow1
ERASE #1
PRINT #1:STR$(2)&",";
PRINT #1:STR$(3)&",";
LET cc=2
FOR r=1 TO 50
LET s=A(r)^2
LET t=A(r+1)^2-1
LET x=1
FOR i=1 TO r
LET x=x*A(i)
NEXT i
FOR n=s TO t
LET p1=0
FOR i=1 TO r
LET p1=p1+MOD(N,A(i))/A(i)
NEXT i
IF DENOM(p1)=x THEN
LET cc=cc+1
PRINT #1:STR$(n)&",";
END IF
NEXT n
NEXT r
CLOSE #1
PRINT n;":";cc
LET TM=TIME-t0
PRINT USING"####." & REPEAT$("#",2):TM;
PRINT "秒"
END
-----------------------------------------
やたらと計算時間がかかります。
ここで、最大公約数を使用しました。
!素数階乗の探究 GCD
DECLARE EXTERNAL FUNCTION GCD
OPTION ARITHMETIC DECIMAL_HIGH !1000桁モード
!OPTION ARITHMETIC RATIONAL !有理数モード
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)
OPEN #1: TextWindow1
ERASE #1
PRINT #1:STR$(2)&",";
PRINT #1:STR$(3)&",";
LET cc=2
FOR r=1 TO 100
LET s=A(r)^2
LET t=A(r+1)^2-1
LET x=1
FOR i=1 TO r
LET x=x*A(i)
NEXT i
FOR n=s TO t
LET p1=GCD(x,n)
IF p1=1 THEN
LET cc=cc+1
PRINT #1:STR$(n)&",";
END IF
NEXT n
NEXT r
CLOSE #1
PRINT n;":";cc
LET TM=TIME-t0
PRINT USING"####." & REPEAT$("#",2):TM;
PRINT "秒"
END
EXTERNAL FUNCTION gcd(a,b)
OPTION ARITHMETIC DECIMAL_HIGH !1000桁モード
DO
LET r=MOD(a,b)
IF r=0 THEN EXIT DO
LET a=b
LET b=r
LOOP
LET gcd=b
END FUNCTION
------------------------------------
数学的知識が乏しいので、素数階乗を使うとGCDで素数判定ができることに驚きました。
まだ30万程度の素数を確かめただけです。先の事は分かりません。
!素数階乗の探究 GCD
OPTION ARITHMETIC RATIONAL !有理数モード
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)
OPEN #1: TextWindow1
ERASE #1
PRINT #1:STR$(2)&",";
PRINT #1:STR$(3)&",";
LET cc=2
FOR r=1 TO 100
LET s=A(r)^2
LET t=A(r+1)^2-1
LET x=1
FOR i=1 TO r
LET x=x*A(i)
NEXT i
FOR n=s TO t
LET p1=GCD(x,n)
IF p1=1 THEN
LET cc=cc+1
PRINT #1:STR$(n)&",";
END IF
NEXT n
NEXT r
CLOSE #1
PRINT n;":";cc
LET TM=TIME-t0
PRINT USING"####." & REPEAT$("#",2):TM;
PRINT "秒"
END
http://blogs.yahoo.co.jp/donald_stinger
|
|