|
|
アトキンのふるい
https://en.wikipedia.org/wiki/Sieve_of_Atkin
https://tjkendev.github.io/procon-library/python/prime/sieve.html
http://fantom1x.blog130.fc2.com/blog-entry-212.html
ネットからの移植版です。エラトステネスの篩より速いとされるアトキンの篩です。
FOR I=1 TO 10
LET N=I*10^5
PRINT N;"th prime";NTHPRIME(N)
NEXT I
END
EXTERNAL FUNCTION NTHPRIME(M)
LET N=M*(LOG(M)+LOG(LOG(M)))
DIM ISPRIME(N)
FOR Z=1 TO 5 STEP 4
FOR Y=Z TO SQR(N) STEP 6
FOR X=1 TO SQR(N)
LET NN=4*X*X+Y*Y
IF NN>N THEN EXIT FOR
LET ISPRIME(NN)=BITNOT(ISPRIME(NN))
NEXT X
FOR X=Y+1 TO SQR(N) STEP 2
LET NN=3*X*X-Y*Y
IF NN>N THEN EXIT FOR
LET ISPRIME(NN)=BITNOT(ISPRIME(NN))
NEXT X
NEXT Y
NEXT Z
FOR Z=2 TO 4 STEP 2
FOR Y=Z TO SQR(N) STEP 6
FOR X=1 TO SQR(N) STEP 2
LET NN=3*X*X+Y*Y
IF NN>N THEN EXIT FOR
LET ISPRIME(NN)=BITNOT(ISPRIME(NN))
NEXT X
FOR X=Y+1 TO SQR(N) STEP 2
LET NN=3*X*X-Y*Y
IF NN>N THEN EXIT FOR
LET ISPRIME(NN)=BITNOT(ISPRIME(NN))
NEXT X
NEXT Y
NEXT Z
FOR Z=1 TO 2
FOR Y=3 TO SQR(N) STEP 6
FOR X=Z TO SQR(N) STEP 3
LET NN=4*X*X+Y*Y
IF NN>N THEN EXIT FOR
LET ISPRIME(NN)=BITNOT(ISPRIME(NN))
NEXT X
NEXT Y
NEXT Z
FOR P=5 TO SQR(N)
IF ISPRIME(P)=-1 THEN
FOR K=P*P TO N STEP P*P
LET ISPRIME(K)=0
NEXT K
END IF
NEXT P
LET ISPRIME(2)=-1
LET ISPRIME(3)=-1
LET COUNT=1
FOR I=3 TO N STEP 2
IF ISPRIME(I)=-1 THEN
LET COUNT=COUNT+1
IF COUNT=M THEN
LET NTHPRIME=I
EXIT FUNCTION
END IF
END IF
NEXT I
END FUNCTION
-------------------------------------------------------------------------------
LET LIMIT=1000000
DIM SIEVE(LIMIT)
LET FACTOR=INT(SQR(LIMIT))+1
FOR I=1 TO FACTOR
FOR J=1 TO FACTOR
LET N=4*I*I+J*J
IF N<=LIMIT AND (MOD(N,12)=1 OR MOD(N,12)=5) THEN LET SIEVE(N)=BITNOT(SIEVE(N))
LET N=3*I*I+J*J
IF N<=LIMIT AND MOD(N,12)=7 THEN LET SIEVE(N)=BITNOT(SIEVE(N))
IF I>J THEN
LET N=3*I*I-J*J
IF N<=LIMIT AND MOD(N,12)=11 THEN LET SIEVE(N)=BITNOT(SIEVE(N))
END IF
NEXT J
NEXT I
FOR I=5 TO FACTOR
IF SIEVE(I)=-1 THEN
FOR J=I*I TO LIMIT STEP I*I
LET SIEVE(J)=0
NEXT J
END IF
NEXT I
LET COUNT=3 ! [2 3 5]
FOR I=7 TO LIMIT STEP 2
IF SIEVE(I)=-1 THEN LET COUNT=COUNT+1
NEXT I
PRINT COUNT
END
-------------------------------------------------------------------------------
LET M=1000000
DIM P(M)
LET SM=INT(SQR(M))
FOR X=1 TO INT(SM/2)
LET V=4*X*X+1
LET Y=8
DO WHILE V<=M
IF MOD(V,12)<>9 THEN LET P(V)=BITXOR(P(V),1)
LET V=V+Y
LET Y=Y+8
LOOP
NEXT X
FOR X=1 TO INT(SM/3^.5) STEP 2
LET V=3*X*X+4
LET Y=12
DO WHILE V<=M
IF MOD(V,12)=7 THEN LET P(V)=BITXOR(P(V),1)
LET V=V+Y
LET Y=Y+8
LOOP
NEXT X
FOR X=2 TO INT(SM/2^.5)
LET V=2*X*(X+1)-1
LET Y=4*X-8
DO WHILE Y>=0 AND V<=M
IF MOD(V,12)=11 THEN LET P(V)=BITXOR(P(V),1)
LET V=V+Y
LET Y=Y-8
LOOP
NEXT X
FOR N=5 TO SM
IF P(N)>0 THEN
FOR Z=N*N TO M-1 STEP N*N
LET P(Z)=0
NEXT Z
END IF
NEXT N
LET P(2)=1
LET P(3)=1
FOR I=2 TO M-1
IF P(I)>0 THEN LET COUNT=COUNT+1
NEXT I
PRINT COUNT
END
|
|