|
投稿者：たろさ
投稿日：2016年 5月18日(水)16時37分3秒
|
|
|
|
素数判定の興味からプログラムしました。
!π(x)素数判定
OPTION ARITHMETIC RATIONAL !有理数モード
LET t0=TIME
LET K9=2000000
DIM A(k9)
OPEN #1:NAME"prime_200.txt",ACCESS INPUT
FOR i=1 TO K9
INPUT #1: A(i)
NEXT i
CLOSE #1
PRINT TIME-t0;"秒で読み込み完了"
PRINT "素数階乗計算開始です。"
LET m1=1
LET m2=1
LET m3=1
LET m4=1
LET m5=1
LET m6=1
LET m7=1
LET m8=1
LET m9=1
LET m10=1
LET m11=1
LET m12=1
LET m13=1
LET m14=1
LET m15=1
LET m16=1
LET m17=1
LET m18=1
LET m19=1
LET m20=1
LET m21=1
LET m22=1
LET m23=1
LET m24=1
LET m25=1
LET m26=1
LET m27=1
LET m28=1
LET m29=1
LET m30=1
LET m31=1
LET m32=1
LET m33=1
LET m34=1
LET m35=1
LET m36=1
LET m37=1
LET m38=1
LET m39=1
LET m40=1
FOR n=2 TO K9
SELECT CASE n
CASE 2 TO 50000
LET m1=m1*A(n)
CASE 50001 TO 100000
LET m2=m2*A(n)
CASE 100001 TO 150000
LET m3=m3*A(n)
CASE 150001 TO 200000
LET m4=m4*A(n)
CASE 200001 TO 250000
LET m5=m5*A(n)
CASE 250001 TO 300000
LET m6=m6*A(n)
CASE 300001 TO 350000
LET m7=m7*A(n)
CASE 350001 TO 400000
LET m8=m8*A(n)
CASE 400001 TO 450000
LET m9=m9*A(n)
CASE 450001 TO 500000
LET m10=m10*A(n)
CASE 500001 TO 550000
LET m11=m11*A(n)
CASE 550001 TO 600000
LET m12=m12*A(n)
CASE 600001 TO 650000
LET m13=m13*A(n)
CASE 650001 TO 700000
LET m14=m14*A(n)
CASE 700001 TO 750000
LET m15=m15*A(n)
CASE 750001 TO 800000
LET m16=m16*A(n)
CASE 800001 TO 850000
LET m17=m17*A(n)
CASE 850001 TO 900000
LET m18=m18*A(n)
CASE 900001 TO 950000
LET m19=m19*A(n)
CASE 950001 TO 1000000
LET m20=m20*A(n)
CASE 1000001 TO 1050000
LET m21=m21*A(n)
CASE 1050001 TO 1100000
LET m22=m22*A(n)
CASE 1100001 TO 1150000
LET m23=m23*A(n)
CASE 1150001 TO 1200000
LET m24=m24*A(n)
CASE 1200001 TO 1250000
LET m25=m25*A(n)
CASE 1250001 TO 1300000
LET m26=m26*A(n)
CASE 1300001 TO 1350000
LET m27=m27*A(n)
CASE 1350001 TO 1400000
LET m28=m28*A(n)
CASE 1400001 TO 1450000
LET m29=m29*A(n)
CASE 1450001 TO 1500000
LET m30=m30*A(n)
CASE 1500001 TO 1550000
LET m31=m31*A(n)
CASE 1550001 TO 1600000
LET m32=m32*A(n)
CASE 1600001 TO 1650000
LET m33=m33*A(n)
CASE 1650001 TO 1700000
LET m34=m34*A(n)
CASE 1700001 TO 1750000
LET m35=m35*A(n)
CASE 1750001 TO 1800000
LET m36=m36*A(n)
CASE 1800001 TO 1850000
LET m37=m37*A(n)
CASE 1850001 TO 1900000
LET m38=m38*A(n)
CASE 1900001 TO 1950000
LET m39=m39*A(n)
CASE 1950001 TO 2000000
LET m40=m40*A(n)
END SELECT
NEXT n
PRINT "素数階乗計算完了です。"
LET C1=29844570422669
FOR n=10^15+1 TO 10^15+10^5-1 STEP 2
LET f=NUMER(n/m1)
IF n=f THEN ELSE GOTO 10
LET f1=NUMER(n/m2)
IF n=f1 THEN ELSE GOTO 10
LET f2=NUMER(n/m3)
IF n=f2 THEN ELSE GOTO 10
LET f3=NUMER(n/m4)
IF n=f3 THEN ELSE GOTO 10
LET f4=NUMER(n/m5)
IF n=f4 THEN ELSE GOTO 10
LET f5=NUMER(n/m6)
IF n=f5 THEN ELSE GOTO 10
LET f6=NUMER(n/m7)
IF n=f6 THEN ELSE GOTO 10
LET f7=NUMER(n/m8)
IF n=f7 THEN ELSE GOTO 10
LET f8=NUMER(n/m9)
IF n=f8 THEN ELSE GOTO 10
LET f9=NUMER(n/m10)
IF n=f9 THEN ELSE GOTO 10
LET f10=NUMER(n/m11)
IF n=f10 THEN ELSE GOTO 10
LET f11=NUMER(n/m12)
IF n=f11 THEN ELSE GOTO 10
LET f12=NUMER(n/m13)
IF n=f12 THEN ELSE GOTO 10
LET f13=NUMER(n/m14)
IF n=f13 THEN ELSE GOTO 10
LET f14=NUMER(n/m15)
IF n=f14 THEN ELSE GOTO 10
LET f15=NUMER(n/m16)
IF n=f15 THEN ELSE GOTO 10
LET f16=NUMER(n/m17)
IF n=f16 THEN ELSE GOTO 10
LET f17=NUMER(n/m18)
IF n=f17 THEN ELSE GOTO 10
LET f18=NUMER(n/m19)
IF n=f18 THEN ELSE GOTO 10
LET f19=NUMER(n/m20)
IF n=f19 THEN ELSE GOTO 10
LET f20=NUMER(n/m21)
IF n=f20 THEN ELSE GOTO 10
LET f21=NUMER(n/m22)
IF n=f21 THEN ELSE GOTO 10
LET f22=NUMER(n/m23)
IF n=f22 THEN ELSE GOTO 10
LET f23=NUMER(n/m24)
IF n=f23 THEN ELSE GOTO 10
LET f24=NUMER(n/m25)
IF n=f24 THEN ELSE GOTO 10
LET f25=NUMER(n/m26)
IF n=f25 THEN ELSE GOTO 10
LET f26=NUMER(n/m27)
IF n=f26 THEN ELSE GOTO 10
LET f27=NUMER(n/m28)
IF n=f27 THEN ELSE GOTO 10
LET f28=NUMER(n/m29)
IF n=f28 THEN ELSE GOTO 10
LET f29=NUMER(n/m30)
IF n=f29 THEN ELSE GOTO 10
LET f30=NUMER(n/m31)
IF n=f30 THEN ELSE GOTO 10
LET f31=NUMER(n/m32)
IF n=f31 THEN ELSE GOTO 10
LET f32=NUMER(n/m33)
IF n=f32 THEN ELSE GOTO 10
LET f33=NUMER(n/m34)
IF n=f33 THEN ELSE GOTO 10
LET f34=NUMER(n/m35)
IF n=f34 THEN ELSE GOTO 10
LET f35=NUMER(n/m36)
IF n=f35 THEN ELSE GOTO 10
LET f36=NUMER(n/m37)
IF n=f36 THEN ELSE GOTO 10
LET f37=NUMER(n/m38)
IF n=f37 THEN ELSE GOTO 10
LET f38=NUMER(n/m39)
IF n=f38 THEN ELSE GOTO 10
LET f39=NUMER(n/m40)
IF n=f39 THEN ELSE GOTO 10
IF n=f AND n=f1 AND n=f2 AND n=f3 AND n=f4 AND n=f5 AND n=f6 AND n=f7 AND n=f8 AND n=f9 AND n=f10 AND n=f11 AND n=f12 AND n=f13 AND n=f14 AND n=f15 AND n=f16 AND n=f17 AND n=f18 AND n=f19 AND n=f20 AND n=f21 AND n=f22 AND n=f23 AND n=f24 AND n=f25 AND n=f26 AND n=f27 AND n=f28 AND n=f29 AND n=f30 AND n=f31 AND n=f32 AND n=f33 AND n=f34 AND n=f35 AND n=f36 AND n=f37 AND n=f38 AND n=f39 THEN
LET C=C+1
PRINT n;":";c+c1;":";c
END IF
10 NEXT n
PRINT c;"個"
LET t1=TIME-t0
PRINT USING "----%.####":t1;
PRINT "秒"
PRINT TIME-t0;"秒で計算しました"
END
------------------------------
素数リスト作成プログラム
------------------------------
!サンダラムの篩
OPTION ARITHMETIC NATIVE !2進モード
LET t0=TIME
!PRINT 2
OPEN #1:NAME "prime_200.txt",RECTYPE INTERNAL
ERASE #1
WRITE #1:2
LET L=32452868/2 !計算結果は2倍の数値になります。
DIM P(L)
LET b=1
LET k=0.7
FOR a=3 TO INT(L*k) STEP 2
LET f=f+1
FOR N=1 TO f STEP 1
LET z=a*n+b
! PRINT STR$(z)&","
IF Z>L THEN GOTO 100
LET P(z)=1
NEXT n
100 LET b=b+1
! PRINT "a=";a;"b=";b-1
NEXT a
FOR i=1 TO L-1
IF P(i)=0 THEN
LET z=I*2+1
WRITE #1:z
END IF
!PRINT i*2+1
NEXT i
PRINT TIME-t0;"秒で計算しました"
END
-------------------------------------
32452843(2000000th prime)
32452867(2000001st prime)
素数判定範囲 32452843+1 TO 32452867^2-1
試し割りよりは効率がよいが、素数が大きくなるに従って、急速に速度が低下するため、実用的ではない。
------------------------------
計算方法の仕組み
------------------------------
OPTION ARITHMETIC RATIONAL !有理数モード
LET c=12
LET m1=3*5*7*11*13*17
LET m2=19*23*29*31*37
FOR n=37 TO 41^2-1 STEP 2
LET x=n/m1
LET f=NUMER(x)
LET y=n/m2
LET f1=NUMER(y)
IF n=f AND n=f1 THEN
LET C=C+1
PRINT n;":";c
END IF
NEXT n
END
----------------------------
プログラムの書き方が不慣れですので、もう少し改良できたらと思い投稿しました。
よろしくお願いします。
http://blogs.yahoo.co.jp/donald_stinger
|
|
|
|
投稿者：しばっち
投稿日：2016年 5月20日(金)20時52分16秒
|
|
|
|
> No.4053[元記事へ]
たろささんへのお返事です。
素数階乗を使うとループ1つ分をはぶくことができますが
素数階乗での計算が多倍長精度なので計算に時間がかかります。
OPTION ARITHMETIC RATIONAL
LET KMAX=500
DIM SOSU(KMAX)
LET SOSU(1)=3
LET SOSU(2)=5
LET SOSU(3)=7
LET K=3
LET P=1
LET NUM=4
FOR I=1 TO KMAX-1
LET P=P*SOSU(I)
FOR J=SOSU(I)^2+2 TO SOSU(I+1)^2-2 STEP 2
IF MOD(J,10)<>5 THEN
IF NUMER(J/P)=J THEN
LET NUM=NUM+1
PRINT NUM;":";J
IF K<KMAX THEN
LET K=K+1
LET SOSU(K)=J
END IF
END IF
END IF
NEXT J
NEXT I
END
ループを使って割り算していくほうが、(素数で)割り切れた時に
(ループを)打ち切りできるのでムダな計算を減らせると思います。
OPTION ARITHMETIC NATIVE
LET KMAX=500
DIM SOSU(KMAX)
LET SOSU(1)=3
LET SOSU(2)=5
LET SOSU(3)=7
LET K=3
LET P=1
LET NUM=4
!'OPEN #1:NAME "sosu.txt"
FOR I=1 TO KMAX-1
FOR J=SOSU(I)^2+2 TO SOSU(I+1)^2-2 STEP 2
IF MOD(J,10)<>5 THEN
LET FL=0
FOR P=1 TO I !'ループ
IF MOD(J,SOSU(P))=0 THEN
LET FL=1
EXIT FOR !'打ち切り
END IF
NEXT P
IF FL=0 THEN
LET NUM=NUM+1
PRINT NUM;":";J
IF K<KMAX THEN
LET K=K+1
LET SOSU(K)=J
!' WRITE #1:J
!'ELSE
!' CLOSE #1
!' STOP
END IF
END IF
END IF
NEXT J
NEXT I
END
|
|
|
|
投稿者：たろさ
投稿日：2016年 5月25日(水)18時10分27秒
|
|
|
|
> No.4054[元記事へ]
しばっちさんへのお返事です。
ありがとうございます。プログラムを理解するまでは至りません(頭が篩=古い)が、自分なりに分析してみました。
テスト機AMD Athlon 64 3800+ (2.4GHz) Win8.1 32bt
BASIC Acc 計算結果
----------------------------
エラトステネスの奇数& 6n±1 篩
1026167 : 80426
249.841 秒
----------------------------
試割&6n±1 primenumber
1026167 : 80426
252.731000000007 秒
----------------------------
Sieve of Sundaram
1026167 : 80426
251.246999999996 秒
----------------------------
!エラトステネスの奇数& 6n±1 篩
OPTION ARITHMETIC NATIVE
LET t0=TIME
PRINT " 2 : 1"
LET Kmin=1 !99991(9592nd prime) 999983(78498th prime)
LET KMAX=Kmin-1+1.E6+26167 !99999989(5761455th prime)
LET C=1 !Kmin以下の素数の個数
DIM s(Kmin TO KMAX)
MAT s=ZER
FOR n=3 TO KMAX STEP 2
IF N>3 AND MOD(n,3)=0 OR N>5 AND MOD(n,10)=5 THEN GOTO 10
IF n < Kmin THEN GOTO 10
IF s(n)=0 THEN
LET C=C+1
PRINT n;":";c
END IF
10 FOR k=n^2 TO KMAX STEP n
IF K < Kmin THEN GOTO 20
LET s(k)=1
20 NEXT k
NEXT n
PRINT TIME-t0;"秒"
END
----------------------------------
!試割&6n±1 primenumber
OPTION ARITHMETIC NATIVE
LET t0=TIME
LET KMAX=169
DIM SOSU(KMAX)
LET SOSU(1)=3
LET SOSU(2)=5
LET SOSU(3)=7
LET NUM=4
LET K=3
FOR I=1 TO KMAX-1
FOR J=SOSU(I)^2+2 TO SOSU(I+1)^2-2 STEP 2
IF MOD(J,3)=0 OR MOD(J,10)=5 THEN GOTO 10
FOR P=3 TO I
IF MOD(J,SOSU(P))=0 THEN GOTO 10
NEXT p
LET NUM=NUM+1
PRINT J;":";NUM
IF K<KMAX THEN
LET K=K+1
LET SOSU(K)=J
END IF
10 NEXT j
NEXT i
PRINT TIME-t0;"秒"
END
----------------------------------
!Sieve of Sundaram
OPTION ARITHMETIC NATIVE
LET t0=TIME
LET L=513084 !計算結果は2倍の数値になります。
LET C=1
DIM P(L)
LET b=1
LET k=0.7
FOR a=3 TO INT(L*k) STEP 2
LET f=f+1
FOR N=1 TO f STEP 1
LET z=a*n+b
IF Z>L THEN GOTO 100
LET P(z)=1
NEXT n
100 LET b=b+1
NEXT a
FOR i=1 TO L-1
IF P(i)=0 THEN
LET z=I*2+1
LET c=c+1
PRINT z;":";c
END IF
NEXT i
PRINT TIME-t0;"秒"
END
-------------------------------------
!エラトステネスの奇数& 6n±1 篩 は1億以上での動作も確認しました。
!エラトステネスの奇数& 6n±1 篩
OPTION ARITHMETIC NATIVE
LET t0=TIME
!PRINT " 2 : 1"
LET Kmin=1.E9 !99991(9592nd prime) 999983(78498th prime)
LET KMAX=Kmin+1.E6 !99999989(5761455th prime) 999999937(50847534th prime)
LET C=50847534 !Kmin以下の素数の個数
DIM s(Kmin TO KMAX)
MAT s=ZER
FOR n=3 TO KMAX STEP 2
IF N>3 AND MOD(n,3)=0 OR N>5 AND MOD(n,10)=5 THEN GOTO 10
IF n < Kmin THEN GOTO 10
IF s(n)=0 THEN
LET C=C+1
PRINT n;":";c
END IF
10 FOR k=n^2 TO KMAX STEP n
IF K < Kmin THEN GOTO 20
LET s(k)=1
20 NEXT k
NEXT n
PRINT TIME-t0;"秒"
END
-------------------------------------
1000999949 : 50895689
276.904999999999 秒
素数階乗 #3996
重宝しています。重ねてお礼申し上げます。
http://blogs.yahoo.co.jp/donald_stinger
|
|
|
戻る