素数階乗と有理数モードを使用した素数判定

 投稿者：たろさ  投稿日：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

  

Re: 素数階乗と有理数モードを使用した素数判定

 投稿者：しばっち  投稿日：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


 		 

Re: 素数階乗と有理数モードを使用した素数判定

 投稿者：たろさ  投稿日：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

  

戻る