素数積篩 360n+k

 投稿者:たろさ  投稿日:2021年11月 8日(月)03時16分56秒
  !360n+k篩
OPTION ARITHMETIC NATIVE
LET t0=TIME

LET k=316241 !MAX 100,008,370,081
LET k2=27294
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 QP=96
   DIM B(QP)
   DATA 1,7,11,13,17,19,23,29,31,37,41,43,47,49,53,59,61,67,71,73
   DATA 77,79,83,89,91,97,101,103,107,109,113,119,121,127,131,133,137,139,143,149
   DATA 151,157,161,163,167,169,173,179,181,187,191,193,197,199,203,209,211,217,221,223
   DATA 227,229,233,239,241,247,251,253,257,259,263,269,271,277,281,283,287,289,293,299
   DATA 301,307,311,313,317,319,323,329,331,337,341,343,347,349,353,359


   MAT READ B

   LET Q=360       !1E8=5761455(0.33秒)
   LET k1=1E8      !1E9=50847534(4.38秒) 1E6=78498
   LET ka=IP(k1/Q) !1E10=455052511(55.36秒)
   LET kb=k1
   LET kc=IP(SQR(k1))
   LET kd=IP(ka/7+7)
   DIM D(ka)


   LET k3=0
   FOR n=1 TO k2
      LET PP=A(n)
      IF PP > kc THEN EXIT FOR
      LET k3=k3+1
   NEXT n
   ! PRINT "k3=";k3;PP;kc

   LET cj =QP-24

   FOR r=1 TO QP
      LET rr=B(r)
      MAT D = ZER

      FOR t=4 TO k3
         LET x=A(t)   !素数篩
         IF x^2 >  k1 THEN EXIT FOR

         IF MOD(x-rr,Q)=0 THEN
            LET y=(x-rr)/Q
            FOR f=1 TO kd
               LET ii=x*f+y
               IF ii>ka THEN EXIT FOR
               LET D(ii)=1
            NEXT f

         ELSE

            FOR i=1 TO QP
               LET rv=B(i)
               IF MOD(x*rv+rr,Q)=0 THEN
                  LET y=-(x*rv+rr)/Q
                  EXIT FOR
               END IF
            NEXT i

            FOR f=1 TO kd
               LET ii=x*f+y
               IF ii > ka THEN EXIT FOR
               IF ii <= 0 THEN GOTO 100
               LET D(ii)=1
100          NEXT f
          END IF
          !PRINT x;kd;f;ii;ka
       NEXT t

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

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

戻る