|
|
!90n+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=24
DIM B(QP)
DATA 1,7,11,13,17,19,23,29,31,37,41,43,47,49,53,59,61,67,71,73,77,79,83,89
MAT READ B
LET Q=90 !1E8=5761455(0.33秒)
LET k1=1E8 !1E9=50847534(4.87秒) 1E6=78498
LET ka=IP(k1/Q) !1E10=455052511(70.92秒)
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
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
--------------------------------------------------------------------
100 ! 素数積n +kの値と個数を求めるprogram
110 LET W=90 !素数積を入力(素数階乗 6,30,210,2310,30030,510510,,,)
120 LET p=w
130 LET z=w
140 OPTION BASE 0 ! DIM文より手前の行にOPTION BASE 0を追加する。
150 DIM A(p)
160 DIM B(z)
170 LET t0=TIME
180 LET p=0 ! p の初期値を0とすると、初回 1になる。140の!を削除
190 LET z=1
191 PRINT"DATA ";
192 LET cc=0
193 LET ct=20
200 FOR k=1 TO w STEP 2 !素数倍を代入
210 !
220 !
230 FOR n=1 TO 20 !素数倍周期 何周でも可能
240 LET m=w*n+k !素数倍を代入
250 FOR i=3 TO SQR(m) STEP 2 !篩
260 IF MOD(m,i)=0 THEN 300
270 NEXT i
280 !
290 LET B(z)=m
300 NEXT n
310 IF B(z)=0 THEN GOTO 370
320 LET A(p)=p
330 LET p=p+1
340 LET z=z+1
350 !
355 PRINT STR$(k);
356 LET cc=cc+1
357 IF cc=ct THEN
358 PRINT
359 ELSE
360 PRINT ",";
361 END IF
362 IF cc=ct THEN
363 !PRINT
364 PRINT"DATA ";
365 LET ct=ct+20
366 END IF
370 NEXT k
380 !
390 PRINT "k=1から数えて";p;"個"
400 PRINT TIME-t0;"秒で計算しました"
410 END
|
|