|
|
!180n+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=48
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
MAT READ B
LET Q=180 !1E8=5761455(0.33秒)
LET k1=1E8 !1E9=50847534(4.72秒) 1E6=78498
LET ka=IP(k1/Q) !1E10=455052511(66.10秒)
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-7
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
|
|