30n± k篩

 投稿者：たろさ  投稿日：2016年12月27日(火)17時26分46秒
  十進BASIC　1億まで
--------------------------------------------------
!30n± k篩
OPTION ARITHMETIC NATIVE
LET t0=TIME

LET k=104743
LET k2=10001
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
   FOR zx=1 TO 9
      PRINT #1:STR$(A(zx))&",";
   NEXT zx

   DIM B(8)
   DATA -1,1,7,11,13,17,19,23
   MAT READ B

   LET ka=3333334
   LET kb=ka*30
   DIM D(ka)
   DIM W(kb)

   FOR r=1 TO 8
      LET rr=B(r)
      MAT D = ZER
      FOR t=4 TO IP(SQR(ka*2))
         LET x=A(t)
         IF x>IP(SQR(kb))THEN EXIT FOR
         IF r=1 THEN
            IF MOD(x+1,30)=0 THEN
               LET y=(x+1)/30
               GOTO 80
            END IF

            IF MOD(x-rr,30)=0 THEN
               LET y=-(x-1)/30
               GOTO 80
            END IF
         ELSEIF r<>1 THEN
            IF MOD(x+rr,30)=0 THEN
               LET y=-(x+rr)/30
               GOTO 80
            END IF

            IF MOD(x-rr,30)=0 THEN
               LET y=(x-rr)/30
               GOTO 80
            END IF
         END IF

         SELECT CASE r
         CASE 1,5,8
            LET ss=3
         CASE 2
            LET ss=2
         CASE 3,4
            LET ss=4
         CASE 6
            LET ss=1
         CASE 7
            LET ss=-1
         END SELECT

         FOR i=-x TO ss !!
            IF MOD(30*i+rr,x)=0 THEN LET y=i
         NEXT i
80
         FOR f=1 TO ka
            IF x*f+y>ka THEN GOTO 90
            IF x*f+y=0  THEN GOTO 90
            LET D(x*f+y)=1
         NEXT f
90    NEXT t

      FOR n=1 TO ka
         IF D(n)=1 THEN GOTO 100
         IF n>ka-1 THEN GOTO 100
         !PRINT #1:STR$(n*30+rr)&",";
         LET W(n*30+rr)=1
         ! LET cj=cj+1
100    NEXT n
    NEXT r

    LET ci=9
    FOR v=29 TO kb STEP 2
       IF W(v)=1 THEN
          LET ci=ci+1
          PRINT #1:STR$(v)&",";
       END IF
    NEXT v

    CLOSE #1
    PRINT ci
    LET TM=TIME-t0
    PRINT USING"####." & REPEAT$("#",2):TM;
    PRINT "秒"
END


--------------------------------------------------
十進BASIC メモリーが足りないと出た場合は、BASIC Acc
--------------------------------------------------
!30n± k篩
OPTION ARITHMETIC NATIVE
LET t0=TIME
LET k=104743
LET k2=10001
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:NAME "prime_10_8v.txt",RECTYPE INTERNAL
   ERASE #1
   FOR zx=1 TO 9
   WRITE #1:A(zx)
   NEXT zx
   DIM B(8)
   DATA -1,1,7,11,13,17,19,23
   MAT READ B

   LET ka=3333334
   LET kb=ka*30
   DIM D(ka)
   DIM W(kb)

   FOR r=1 TO 8
      LET rr=B(r)
      MAT D = ZER
      FOR t=4 TO IP(SQR(ka*2))
         LET x=A(t)
         IF x>IP(SQR(kb))THEN EXIT FOR
         IF r=1 THEN
            IF MOD(x+1,30)=0 THEN
               LET y=(x+1)/30
               GOTO 80
            END IF

            IF MOD(x-rr,30)=0 THEN
               LET y=-(x-1)/30
               GOTO 80
            END IF
         ELSEIF r<>1 THEN
            IF MOD(x+rr,30)=0 THEN
               LET y=-(x+rr)/30
               GOTO 80
            END IF

            IF MOD(x-rr,30)=0 THEN
               LET y=(x-rr)/30
               GOTO 80
            END IF
         END IF

         SELECT CASE r
         CASE 1,5,8
            LET ss=3
         CASE 2
            LET ss=2
         CASE 3,4
            LET ss=4
         CASE 6
            LET ss=1
         CASE 7
            LET ss=-1
         END SELECT

         FOR i=-x TO ss !!
            IF MOD(30*i+rr,x)=0 THEN LET y=i
         NEXT i
80
         FOR f=1 TO ka
            IF x*f+y>ka THEN GOTO 90
            IF x*f+y=0  THEN GOTO 90
            LET D(x*f+y)=1
         NEXT f
90    NEXT t

      FOR n=1 TO ka
         IF D(n)=1 THEN GOTO 100
         IF n>ka-1 THEN GOTO 100
         !PRINT #1:STR$(n*30+rr)&",";
         LET W(n*30+rr)=1
         ! LET cj=cj+1
100    NEXT n
    NEXT r

    LET ci=9
    FOR v=29 TO kb STEP 2
       IF W(v)=1 THEN
          LET ci=ci+1
          WRITE #1:v
       END IF
    NEXT v

    CLOSE #1
    PRINT ci
    LET TM=TIME-t0
    PRINT USING"####." & REPEAT$("#",2):TM;
    PRINT "秒"
END

---------------------------------------------------
6n±1篩　1億まで　最速設定　30n± k篩より数秒速い。
---------------------------------------------------
!6n±1篩
!A046954 Numbers n such that 6n+1 is nonprime. http://oeis.org/A046954
OPTION ARITHMETIC NATIVE       !2進モード
LET t0=TIME
LET k=1E4
LET k2=1229
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 k6=1E8          !篩の計算範囲
   LET k5=k6/6+1
   DIM A6(k5),A7(k6)

   MAT A6 = ZER     !(6*n-1)
   FOR n=1 TO k2-2
      LET P6=A(n+2)
      IF MOD(P6+1,6)=0 THEN
         LET r6=(P6+1)/6
         FOR i=1 TO k5
            IF P6*i+r6>k5 THEN EXIT FOR
            LET A6(P6*i+r6)=1
         NEXT i
      END IF

      IF MOD(P6-1,6)=0 THEN
         LET r6=(P6-1)/6
         FOR i=1 TO k5
            IF P6*i-r6>k5 THEN EXIT FOR
            LET A6(P6*i-r6)=1
         NEXT i
      END IF
   NEXT n

   FOR n=1 TO k5
      IF A6(n)=1 THEN GOTO 100
      LET A7(6*n-1)=1
100 NEXT n

    MAT A6 = ZER !(6*n+1)
    FOR n=1 TO k2-2
       LET P6=A(n+2)
       IF MOD(P6+1,6)=0 THEN
          LET r6=(P6+1)/6
          FOR i=1 TO k5
             IF P6*i-r6>k5 THEN EXIT FOR
             LET A6(P6*i-r6)=1
          NEXT i
       END IF

       IF MOD(P6-1,6)=0 THEN
          LET r6=(P6-1)/6
          FOR i=1 TO k5
             IF P6*i+r6>k5 THEN EXIT FOR
             LET A6(P6*i+r6)=1
          NEXT i
       END IF
    NEXT n

    FOR n=1 TO k5
       IF A6(n)=1 THEN GOTO 200
       LET A7(6*n+1)=1
200 NEXT n

    OPEN #1: TextWindow1
    ERASE #1

    PRINT #1:STR$(2)&",";
    PRINT #1:STR$(3)&",";
    FOR n=1 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+2

    LET TM=TIME-t0
    PRINT USING"####." & REPEAT$("#",2):TM;
    PRINT "秒"
END

http://blogs.yahoo.co.jp/donald_stinger

  

Re: 30n± k篩

 投稿者：たろさ  投稿日：2017年 2月26日(日)22時18分54秒
  > No.4208[元記事へ]

素数計数関数を求めるプログラム

4n+k,6n+k,10n+k,16n+k,30n+k,60n+k,100n+k,210n+k

テストの結果(1億)　上記中一番計算時間が短い。

!30n+k篩  50億までの素数を数える。π(5E9)234954223個
OPTION ARITHMETIC NATIVE
LET t0=TIME

LET kb=5E9
LET k6=70709 !60億,77447 10億,31607 1億,9973
LET k2=7004  !60億,7608  10億,3401  1億,1229
!エラトステネスの篩
LET Fu=k6
DIM P(Fu)
DIM A(k2) !素数
MAT P=ZER
MAT A=ZER
LET A(1)=2
LET H1=1
FOR I=3 TO SQR(Fu) STEP 2
   IF P(I)=0 THEN
      FOR J=I*I TO Fu STEP I
         LET P(J)=1
      NEXT J
   END IF
NEXT I
FOR I=3 TO Fu STEP 2
   IF P(I)=0 THEN
      LET H1=H1+1
      LET A(H1)=I
   END IF
NEXT I

DIM B(8)
DATA 1,7,11,13,17,19,23,29
MAT READ B

LET Q=30
LET U=IP(kb/Q)+1
LET kd=IP(kb/207)

DIM D(U)

LET cj=10

FOR r=1 TO 8
   LET rr=B(r)
   MAT D = ZER
   FOR t=4 TO k2
      LET x=A(t)

      IF MOD(x+rr,Q)=0 THEN
         LET y=-(x+rr)/Q
         GOTO 80
      END IF

      IF MOD(x-rr,Q)=0 THEN
         LET y=(x-rr)/Q
         GOTO 80
      END IF

      FOR i=0 TO x
         IF MOD(Q*i-rr,x)=0 THEN
            LET y=-i
            GOTO 80
            EXIT FOR
         END IF
      NEXT i
80       FOR f=1 TO kd
            IF x*f+y>U THEN EXIT FOR
            LET D(x*f+y)=1
         NEXT f
      NEXT t

      FOR n=1 TO U
         IF D(n)=0 THEN
            IF n*Q+rr>kb THEN EXIT FOR
            LET cj=cj+1
         END IF
      NEXT n
   NEXT r

   PRINT cj
   LET TM=TIME-t0
   PRINT USING"####." & REPEAT$("#",2):TM;
   PRINT "秒"
END

----------------------------------
課題

不完全なループ

      FOR i=0 TO x
         IF MOD(Q*i-rr,x)=0 THEN
            LET y=-i
            GOTO 80
            EXIT FOR
         END IF
      NEXT i

素数増大と共に計算量が増える。
-----------------------------------
4n+k,6n+kには不用です。

http://blogs.yahoo.co.jp/donald_stinger

  
 

戻る