|
投稿者：たろさ
投稿日：2016年12月22日(木)05時53分9秒
|
|
|
|
!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=1E7 !篩の計算範囲
DIM A6(k6),A7(k6)
MAT A6 = ZER !(6*n-1)
FOR n=1 TO k2-2
LET P6=n+2
IF MOD(A(P6)+1,6)=0 THEN
LET r6=(A(P6)+1)/6
FOR i=1 TO k6
IF A(P6)*i+r6>k6/6 THEN EXIT FOR
LET A6(A(P6)*i+r6)=1
NEXT i
END IF
IF MOD(A(P6)-1,6)=0 THEN
LET r6=(A(P6)-1)/6
FOR i=1 TO k6
IF A(P6)*i-r6>k6/6 THEN EXIT FOR
LET A6(A(P6)*i-r6)=1
NEXT i
END IF
NEXT n
FOR n=1 TO k6/6
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=n+2
IF MOD(A(P6)+1,6)=0 THEN
LET r6=(A(P6)+1)/6
FOR i=1 TO k6
IF A(P6)*i-r6>k6/6 THEN EXIT FOR
LET A6(A(P6)*i-r6)=1
NEXT i
END IF
IF MOD(A(P6)-1,6)=0 THEN
LET r6=(A(P6)-1)/6
FOR i=1 TO k6
IF A(P6)*i+r6>k6/6 THEN EXIT FOR
LET A6(A(P6)*i+r6)=1
NEXT i
END IF
NEXT n
FOR n=1 TO k6/6
IF A6(n)=1 THEN GOTO 200
LET A7(6*n+1)=1
200 NEXT n
OPEN #1: TextWindow1
ERASE #1
PRINT #1:2;3
FOR n=1 TO k6
IF A7(n)=1 THEN
LET c1=c1+1
PRINT #1:STR$(n)&",";
END IF
NEXT n
PRINT c1+2
LET TM=TIME-t0
PRINT USING"####." & REPEAT$("#",2):TM;
PRINT "秒"
END
------------------------------------------
数値分析用プログラム
!6, 11, 13, 16, 20, 21, 24, 26, 27, 31, 34, 35, 36, 37, 41, 46, 48, 50
OPTION ARITHMETIC NATIVE !2進モード
LET k=2012
LET k2=305
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=1E4
DIM A6(k6)
FOR n=1 TO k2-2
LET P6=n+2
IF MOD(A(P6)+1,6)=0 THEN
LET r6=(A(P6)+1)/6
FOR i=1 TO k6
IF A(P6)*i+r6>k6 THEN EXIT FOR
LET A6(A(P6)*i+r6)=1
NEXT i
END IF
IF MOD(A(P6)-1,6)=0 THEN
LET r6=(A(P6)-1)/6
FOR i=1 TO k6
IF A(P6)*i-r6>k6 THEN EXIT FOR
LET A6(A(P6)*i-r6)=1
NEXT i
END IF
NEXT n
FOR n=1 TO k6
LET m=6*n-1
IF ISPRIME(m)=0 THEN PRINT TAB(6);n;
IF A6(n)=1 THEN PRINT n!;
NEXT n
END
EXTERNAL FUNCTION ISPRIME(N)
OPTION ARITHMETIC NATIVE
IF N=2 OR N=3 OR N=5 OR N=7 OR N=11 OR N=13 OR N=17 THEN
LET ISPRIME=1
EXIT FUNCTION
END IF
IF N=1 OR MOD(N,2)=0 THEN
LET ISPRIME=0
EXIT FUNCTION
END IF
LET D=(N-1)/2
LET S=0
DO WHILE MOD(D, 2)=0
LET D=INT(D/2)
LET S=S+1
LOOP
FOR I=1 TO 7
LET ISP=0
READ A !' n < 341550071728321 なら a=2,3,5,7,11,13,17
DATA 2,3,5,7,11,13,17
LET ISP=0
LET R=POWMOD(A,D,N)
IF R=1 OR R=N-1 THEN
LET ISP=1
END IF
LET R=POWMOD(R,2,N)
FOR J=0 TO S-1
IF R=N-1 THEN
LET ISP=1
END IF
LET R=POWMOD(R,2,N)
NEXT J
IF ISP=0 THEN
LET ISPRIME=0
EXIT FUNCTION
END IF
NEXT I
LET ISPRIME=1
END FUNCTION
EXTERNAL FUNCTION POWMOD(B,P,M)
OPTION ARITHMETIC NATIVE
LET RESULT=1
DO WHILE P>0
IF MOD(P,2)=1 THEN
LET RESULT=MOD(RESULT*B,M)
END IF
LET B=MOD(B*B,M)
LET P=INT(P/2)
LOOP
LET POWMOD=RESULT
END FUNCTION
------------------------------------------
数値分析用プログラム
!4, 8, 9, 14, 15, 19, 20, 22, 24, 28, 29, 31, 34, 36, 39, 41, 42, 43, 44, 48, 49, 50
OPTION ARITHMETIC NATIVE !2進モード
LET k=2012
LET k2=305
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=1E4
DIM A6(k6)
FOR n=1 TO k2-2
LET P6=n+2
IF MOD(A(P6)+1,6)=0 THEN
LET r6=(A(P6)+1)/6
FOR i=1 TO k6
IF A(P6)*i-r6>k6 THEN EXIT FOR
LET A6(A(P6)*i-r6)=1
NEXT i
END IF
IF MOD(A(P6)-1,6)=0 THEN
LET r6=(A(P6)-1)/6
FOR i=1 TO k6
IF A(P6)*i+r6>k6 THEN EXIT FOR
LET A6(A(P6)*i+r6)=1
NEXT i
END IF
NEXT n
FOR n=1 TO k6
LET m=6*n+1
IF ISPRIME(m)=0 THEN PRINT TAB(6);n;
IF A6(n)=1 THEN PRINT n!;
NEXT n
END
EXTERNAL FUNCTION ISPRIME(N)
OPTION ARITHMETIC NATIVE
IF N=2 OR N=3 OR N=5 OR N=7 OR N=11 OR N=13 OR N=17 THEN
LET ISPRIME=1
EXIT FUNCTION
END IF
IF N=1 OR MOD(N,2)=0 THEN
LET ISPRIME=0
EXIT FUNCTION
END IF
LET D=(N-1)/2
LET S=0
DO WHILE MOD(D, 2)=0
LET D=INT(D/2)
LET S=S+1
LOOP
FOR I=1 TO 7
LET ISP=0
READ A !' n < 341550071728321 なら a=2,3,5,7,11,13,17
DATA 2,3,5,7,11,13,17
LET ISP=0
LET R=POWMOD(A,D,N)
IF R=1 OR R=N-1 THEN
LET ISP=1
END IF
LET R=POWMOD(R,2,N)
FOR J=0 TO S-1
IF R=N-1 THEN
LET ISP=1
END IF
LET R=POWMOD(R,2,N)
NEXT J
IF ISP=0 THEN
LET ISPRIME=0
EXIT FUNCTION
END IF
NEXT I
LET ISPRIME=1
END FUNCTION
EXTERNAL FUNCTION POWMOD(B,P,M)
OPTION ARITHMETIC NATIVE
LET RESULT=1
DO WHILE P>0
IF MOD(P,2)=1 THEN
LET RESULT=MOD(RESULT*B,M)
END IF
LET B=MOD(B*B,M)
LET P=INT(P/2)
LOOP
LET POWMOD=RESULT
END FUNCTION
http://blogs.yahoo.co.jp/donald_stinger
|
|
|
|
投稿者：たろさ
投稿日：2017年 2月11日(土)23時23分20秒
|
|
|
|
> No.4207[元記事へ]
たろささんへのお返事です。
> !6n±1篩
!6n±1篩 1億までの素数をカウント
OPTION ARITHMETIC NATIVE !2進モード
LET t0=TIME
LET k=1E4 !SQR(1E8)
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 EXIT SUB
END IF
FOR k1=n1 TO k9 STEP 2
LET m1=n1*k1
IF m1>k9 THEN EXIT FOR
LET P(m1)=1
NEXT k1
NEXT n1
END SUB
CALL prime(k)
LET k6=1E8 !篩の計算範囲
LET k5=IP(k6/6)+1
DIM A6(k5),A7(k5)
MAT A6 = ZER !(6*n-1)
FOR n=3 TO k2
LET P6=A(n)
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
MAT A7 = ZER !(6*n+1)
FOR n=3 TO k2
LET P6=A(n)
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 A7(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 A7(P6*i+r6)=1
NEXT i
END IF
NEXT n
FOR n=1 TO k5
IF A6(n)=1 THEN GOTO 100
LET c1=c1+1
100 IF A7(n)=1 THEN GOTO 200
LET c1=c1+1
200 NEXT n
PRINT c1+2
LET TM=TIME-t0
PRINT USING"####." & REPEAT$("#",2):TM;
PRINT "秒"
END
-----------------------------------------
素数Counter 1億π(x)5761455
AMD 10進BASIC 28.01秒
AMD BASIC Acc 8.49秒
Intel i7 4.5GHz 1.67秒
-----------------------------------------
■ パソコン環境
ウィンドウズ : Microsoft Windows 8.1
システムの種類 : 32 ビット
プロセッサー : AMD Athlon(tm) 64 Processor 3800+
周波数 : 2.40 GHz
メモリー : 3.00 GB
-----------------------------------------
■ パソコン環境
ウィンドウズ : Microsoft Windows 7
システムの種類 : 32 ビット
メモリー : 3.00 GB
プロセッサー :Intel Core i7 4790K
-----------------------------------------
http://blogs.yahoo.co.jp/donald_stinger
|
|
|
|
投稿者：たろさ
投稿日：2017年 2月11日(土)23時32分23秒
|
|
|
|
> No.4262[元記事へ]
たろささんへのお返事です。
DECLARE EXTERNAL FUNCTION cprime
OPTION ARITHMETIC NATIVE !2進モード
LET t0=TIME
LET k6=1E8
DATA 14143,17321,20011,22367,24499,26459,28289,30011,31627,33179
DATA 34649,36061,37423,38729,40009,41231,42433,43591,44729,45827
DATA 46919,47963,48989,50021,50993,51971,52919,53857,54773,55681
DATA 56569,57457,58309,59167,60013,60859,61651,62459,63247,64033
DATA 64811,65579,66337,67103,67829,68567,69313,70001,70717,71419
DATA 72139,72817,73517,74161,74843,75503,76157,76819,77471,78121
DATA 78779,79379,80021,80627,81281,81853,82463,83071,83689,84263
DATA 84857,85447,86027,86627,87179,87751,88321,88883,89443,90001
DATA 90583,91121,91673,92203,92737,93281,93809,94343,94873,95393
DATA 95917,96443,96953,97499,97987,98491,98999,99523,100003,100501
DIM Ba(100)
MAT READ Ba
DATA 1664,1991,2263,2503,2718,2907,3081,3246,3402,3555
DATA 3700,3830,3962,4082,4204,4317,4436,4541,4649,4747
DATA 4847,4943,5034,5134,5222,5316,5400,5486,5572,5651
DATA 5737,5824,5903,5983,6058,6132,6203,6270,6338,6417
DATA 6478,6549,6614,6685,6757,6818,6881,6936,7005,7074
DATA 7141,7199,7258,7316,7380,7436,7497,7552,7609,7671
DATA 7725,7779,7838,7889,7953,8004,8059,8113,8167,8215
DATA 8266,8317,8365,8422,8465,8523,8564,8611,8661,8714
DATA 8768,8810,8860,8904,8959,9010,9051,9099,9148,9195
DATA 9244,9287,9331,9375,9418,9459,9505,9548,9593,9633
DIM Bb(100)
MAT READ Bb
DATA 1177,1922,2497,2979,3400,3778,4122,4441,4741,5022
DATA 5289,5544,5788,6023,6247,6466,6676,6880,7078,7272
DATA 7458,7643,7823,7995,8169,8336,8502,8663,8823,8978
DATA 9132,9281,9431,9576,9718,9857,10001,10138,10276,10410
DATA 10542,10672,10802,10927,11056,11180,11300,11427,11547,11667
DATA 11781,11900,12014,12134,12246,12360,12473,12582,12691,12799
DATA 12914,13021,13124,13232,13337,13438,13552,13642,13741,13848
DATA 13944,14043,14142,14237,14338,14434,14530,14625,14720,14814
DATA 14903,14998,15089,15183,15276,15365,15456,15545,15634,15724
DATA 15812,15916,15986,16068,16158,16244,16329,16411,16411,16499
DIM Bc(100)
MAT READ Bc
LET c1=5761455
PRINT " 1";":";c1
FOR n=1 TO 100
LET k=Ba(n) !SQR(1E8)
LET k2=Bb(n) !π(x)
LET S=Bc(n)
LET k6=k6+1E8
LET C=cprime(k,k2,k6,S)
LET c1=c1+c
PRINT n+1;":";c1;":";c
NEXT n
LET TM=TIME-t0
PRINT USING"####." & REPEAT$("#",2):TM;
PRINT "秒"
END
EXTERNAL FUNCTION cprime (k,k2,k6,S)
OPTION ARITHMETIC NATIVE
LET k4=k6-1E8
LET U=IP(k6/6)+1
LET W=IP(k4/6)-1
DIM x(w TO u)
DIM y(w TO u)
!エラトステネスの篩
DIM P(k)
DIM A(k2) !素数
LET h1=1
LET A(h1)=2
LET h1=2
FOR n1=3 TO k STEP 2
IF P(n1)=0 THEN
LET A(h1)=n1
LET h1=h1+1
IF h1>k+1 THEN GOTO 20
END IF
FOR k1=n1 TO k STEP 2
LET m1=n1*k1
IF m1>k THEN GOTO 10
LET P(m1)=1
NEXT k1
10 NEXT n1
20
FOR n=3 TO k2
LET P6=A(n)
IF MOD(P6+1,6)=0 THEN !(6*n-1)
LET r=(P6+1)/6
FOR i=S TO u
IF P6*i+r>u THEN EXIT FOR
IF P6*i+r<w THEN GOTO 30
LET x(P6*i+r)=1
30 NEXT i
END IF
IF MOD(P6-1,6)=0 THEN
LET r=(P6-1)/6
FOR i=S TO u
IF P6*i-r>u THEN EXIT FOR
IF P6*i-r<w THEN GOTO 40
LET x(P6*i-r)=1
40 NEXT i
END IF
IF MOD(P6+1,6)=0 THEN !(6*n+1)
LET r=(P6+1)/6
FOR i=S TO u
IF P6*i-r>u THEN EXIT FOR
IF P6*i-r<w THEN GOTO 50
LET y(P6*i-r)=1
50 NEXT i
END IF
IF MOD(P6-1,6)=0 THEN
LET r=(P6-1)/6
FOR i=S TO u
IF P6*i+r>u THEN EXIT FOR
IF P6*i+r<w THEN GOTO 60
LET y(P6*i+r)=1
60 NEXT i
END IF
NEXT n
FOR n=w TO u
IF x(n)=1 THEN
GOTO 100
ELSE
IF k6-1E8>n*6-1 OR k6<n*6-1 THEN GOTO 100
LET cc=cc+1
END IF
100 IF y(n)=1 THEN
GOTO 200
ELSE
IF k6-1E8>n*6+1 OR k6<n*6+1 THEN GOTO 200
LET cc=cc+1
END IF
200 NEXT n
LET cprime=cc
END FUNCTION
-----------------------------------------
10億までは高速です。
-----------------------------------------
6n±1篩を途中から計算するプログラムです。
10億までは、しばっち様の
Re: エラトステネスの篩
#4101
プログラムとBASIC Accの計算結果ならいい勝負です。
素数リスト作成プログラムは、下記の方が速い。
1000億までの素数リスト作成
Re: Sieve of Sundaram サンダラムの篩(ふるい)
6n±1篩を途中から計算する方法
教えてください。
http://blogs.yahoo.co.jp/donald_stinger
|
|
|
|
投稿者：しばっち
投稿日：2017年 2月12日(日)14時33分32秒
|
|
|
|
> No.4263[元記事へ]
たろささんへのお返事です。
6n±1篩というのがどういうものかよくわかりませんが
たろささんのエラトステネスの篩はムダな処理が多いようです
LET T=TIME
LET k=1e+8
LET K2=1E+7
DIM P(k)
DIM A(k2) !素数
LET h1=1
LET A(h1)=2
LET h1=2
FOR n1=3 TO k STEP 2
IF P(n1)=0 THEN
LET A(h1)=n1
LET h1=h1+1
IF h1>k+1 THEN GOTO 20
END IF
FOR k1=n1 TO k STEP 2
LET m1=n1*k1
IF m1>k THEN GOTO 10
LET P(m1)=1
NEXT k1
10 NEXT n1
20 PRINT TIME-T
PRINT H1
MAT P=ZER
MAT A=ZER
LET T=TIME
LET H1=1
LET A(H1)=2
LET H1=2
FOR I=3 TO K STEP 2
IF P(I)=0 THEN
LET A(H1)=I
LET H1=H1+1
FOR J=I*I TO K STEP I
LET P(J)=1
NEXT J
END IF
NEXT I
PRINT TIME-T
PRINT H1
MAT P=ZER
MAT A=ZER
LET T=TIME
LET H1=1
LET A(H1)=2
LET H1=2
FOR I=3 TO SQR(K) STEP 2
IF P(I)=0 THEN
FOR J=I*I TO K STEP I
LET P(J)=1
NEXT J
END IF
NEXT I
FOR I=3 TO K STEP 2
IF P(I)=0 THEN
LET A(H1)=I
LET H1=H1+1
END IF
NEXT I
PRINT TIME-T
PRINT H1
END
|
|
|
|
投稿者：たろさ
投稿日：2017年 2月12日(日)19時01分51秒
|
|
|
|
しばっちさんへのお返事です。
お世話になります。エラトステネスの篩を組み込みました。6m±1篩を探究しています。
!6n±1篩 1億までの素数をカウント
OPTION ARITHMETIC NATIVE !2進モード
LET t0=TIME
LET k=1E4 !SQR(1E8)
LET k2=1229
DIM P(k)
DIM A(k2) !素数
SUB prime(k) !エラトステネスの篩
MAT P=ZER
MAT A=ZER
LET H1=1
LET A(H1)=2
LET H1=2
FOR I=3 TO SQR(K) STEP 2
IF P(I)=0 THEN
FOR J=I*I TO K STEP I
LET P(J)=1
NEXT J
END IF
NEXT I
FOR I=3 TO K STEP 2
IF P(I)=0 THEN
LET A(H1)=I
LET H1=H1+1
END IF
NEXT I
END SUB
CALL prime(k)
LET k6=1E8 !篩の計算範囲
LET k5=IP(k6/6)+1
DIM A6(k5),A7(k5)
MAT A6 = ZER !(6*n-1)
FOR n=3 TO k2
LET P6=A(n)
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
MAT A7 = ZER !(6*n+1)
FOR n=3 TO k2
LET P6=A(n)
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 A7(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 A7(P6*i+r6)=1
NEXT i
END IF
NEXT n
FOR n=1 TO k5
IF A6(n)=1 THEN GOTO 100
LET c1=c1+1
100 IF A7(n)=1 THEN GOTO 200
LET c1=c1+1
200 NEXT n
PRINT c1+2
LET TM=TIME-t0
PRINT USING"####." & REPEAT$("#",2):TM;
PRINT "秒"
END
-------------------------------------------------
6n±1篩は、サンダラムの篩の出力部 2n+1 の例に倣って
6n±1 のnを篩っています。
EXTERNAL FUNCTION PrimeQ(n)のネットのダウンロードサイトlink切れです。
5+6n = 6n-1 , 7+6n = 6n+1 とても効率的な方法です。
DECLARE EXTERNAL FUNCTION count
DECLARE EXTERNAL FUNCTION PrimeQ
FOR n=1 TO 100
LET S=n*100-100
IF S=0 THEN LET S=1
LET T=n*100
PRINT S;"TO";T
PRINT count(S,T);"個"
PRINT
NEXT n
END
EXTERNAL FUNCTION count(S,T)
DIM A(S TO T)
MAT A=ZER
FOR n=S TO T
LET v=n+1
IF MOD(v,6)=0 AND PrimeQ(n)=1 THEN
LET C=C+1
LET A(n)=1
PRINT c;":";n
END IF
NEXT n
PRINT
LET count=C
END FUNCTION
EXTERNAL FUNCTION PrimeQ(n) !素数判定 1:素数、0:素数でない
LET PrimeQ=0
IF n<2 OR n<>INT(n) THEN EXIT FUNCTION !引数を確認する
!２以上の自然数なら
IF MOD(n,2)=0 THEN !２の倍数
IF n=2 THEN LET PrimeQ=1 !２は素数
ELSEIF MOD(n,3)=0 THEN !３の倍数
IF n=3 THEN LET PrimeQ=1 !３は素数
ELSE
LET k=5
DO WHILE k*k<=n !√nまで検証する
IF MOD(n,k)=0 THEN !5,11,17,23,29,…
EXIT FUNCTION
ELSEIF MOD(n,k+2)=0 THEN !7,13,19,25,31,…
EXIT FUNCTION
END IF !+1,+3,+5は２の倍数（偶数）、+1,+4は３の倍数、+5は５の倍数
LET k=k+6
LOOP
LET PrimeQ=1 !最後まで割り切れなければ、素数である
END IF
END FUNCTION
配列の再定義の方法については、外部関数定義だけの知識です。
配列を共有する方法があったらと思います。
http://blogs.yahoo.co.jp/donald_stinger
|
|
|
|
投稿者：たろさ
投稿日：2017年 2月13日(月)01時40分48秒
|
|
|
|
> No.4262[元記事へ]
たろささんへのお返事です。
1億では変わらずなので3億までの素数をカウントしました。
BASIC Acc 専用です。
------------------------------------------------
!6n±1篩 3億までの素数をカウント
OPTION ARITHMETIC NATIVE !2進モード
LET t0=TIME
LET k=17321 !SQR(3E8)
LET k2=1991
DIM P(k)
DIM A(k2) !素数
SUB prime(k) !エラトステネスの篩
MAT P=ZER
MAT A=ZER
LET H1=1
LET A(H1)=2
LET H1=2
FOR I=3 TO SQR(K) STEP 2
IF P(I)=0 THEN
FOR J=I*I TO K STEP I
LET P(J)=1
NEXT J
END IF
NEXT I
FOR I=3 TO K STEP 2
IF P(I)=0 THEN
LET A(H1)=I
LET H1=H1+1
END IF
NEXT I
END SUB
CALL prime(k)
LET k6=3E8 !篩の計算範囲
LET B6=IP(k6/30)+1
LET k5=IP(k6/6)+1
DIM A6(k5),A7(k5)
MAT A6 = ZER !(6*n-1)
MAT A7 = ZER !(6*n+1)
FOR n=3 TO k2
LET P6=A(n)
IF MOD(P6+1,6)=0 THEN
LET r6=(P6+1)/6
CALL nap(P6,r6)
CALL napi(P6,r6*(-1))
END IF
IF MOD(P6-1,6)=0 THEN
LET r6=(P6-1)/6
CALL nap(P6,r6*(-1))
CALL napi(P6,r6)
END IF
NEXT n
SUB nap(z,r6)
FOR i=1 TO B6
IF z*i+r6>k5 THEN EXIT SUB
LET A6(z*i+r6)=1
NEXT i
END SUB
SUB napi(z,r6)
FOR i=1 TO B6
IF z*i+r6>k5 THEN EXIT SUB
LET A7(z*i+r6)=1
NEXT i
END SUB
LET C1=2
FOR n=1 TO k5
IF A6(n)=1 THEN GOTO 100
LET c1=c1+1
IF c1>=16252325 THEN PRINT n*6-1
100 IF A7(n)=1 THEN GOTO 200
LET c1=c1+1
IF c1>=16252325 THEN PRINT n*6+1
200 NEXT n
PRINT c1
LET TM=TIME-t0
PRINT USING"####." & REPEAT$("#",2):TM;
PRINT "秒"
END
------------------------------------------------
計算結果
299999977
300000007
16252326 一つ余計に出ました。ので念のため素数も出力しました。
素数Counter 3億π(x)299999977(16252325th prime)
Intel i7 4.5GHz 5.02秒
-----------------------------------------
■ パソコン環境
ウィンドウズ : Microsoft Windows 7
システムの種類 : 32 ビット
メモリー : 3.00 GB
プロセッサー :Intel Core i7 4790K
-----------------------------------------
!6n±1篩 6億までの素数をカウント
OPTION ARITHMETIC NATIVE !2進モード
LET t0=TIME
LET k=24499 !SQR(3E8)
LET k2=2718
DIM P(k)
DIM A(k2) !素数
SUB prime(k) !エラトステネスの篩
MAT P=ZER
MAT A=ZER
LET H1=1
LET A(H1)=2
LET H1=2
FOR I=3 TO SQR(K) STEP 2
IF P(I)=0 THEN
FOR J=I*I TO K STEP I
LET P(J)=1
NEXT J
END IF
NEXT I
FOR I=3 TO K STEP 2
IF P(I)=0 THEN
LET A(H1)=I
LET H1=H1+1
END IF
NEXT I
END SUB
CALL prime(k)
LET k6=6E8 !篩の計算範囲
LET B6=IP(k6/30)+1
LET k5=IP(k6/6)+1
DIM A6(k5),A7(k5)
MAT A6 = ZER !(6*n-1)
MAT A7 = ZER !(6*n+1)
FOR n=3 TO k2
LET P6=A(n)
IF MOD(P6+1,6)=0 THEN
LET r6=(P6+1)/6
CALL nap(P6,r6)
CALL napi(P6,r6*(-1))
END IF
IF MOD(P6-1,6)=0 THEN
LET r6=(P6-1)/6
CALL nap(P6,r6*(-1))
CALL napi(P6,r6)
END IF
NEXT n
SUB nap(z,r6)
FOR i=1 TO B6
IF z*i+r6>k5 THEN EXIT SUB
LET A6(z*i+r6)=1
NEXT i
END SUB
SUB napi(z,r6)
FOR i=1 TO B6
IF z*i+r6>k5 THEN EXIT SUB
LET A7(z*i+r6)=1
NEXT i
END SUB
LET C1=2
FOR n=1 TO k5
IF A6(n)=1 THEN GOTO 100
LET c1=c1+1
IF c1>=31324703 THEN PRINT n*6-1
100 IF A7(n)=1 THEN GOTO 200
LET c1=c1+1
IF c1>=31324703 THEN PRINT n*6+1
200 NEXT n
PRINT c1
LET TM=TIME-t0
PRINT USING"####." & REPEAT$("#",2):TM;
PRINT "秒"
END
------------------------------------------
計算結果
599999971
600000001
600000007
31324705
10.22秒
http://blogs.yahoo.co.jp/donald_stinger
|
|
|
|
投稿者：たろさ
投稿日：2017年 2月16日(木)00時07分5秒
|
|
|
|
> No.4207[元記事へ]
仮称6n±1篩
プログラムの目的
6n-1,6n+1は、確率的素数です。素数に成らないnを篩、素数を生成します。
nの篩方は、私が偶然発見したものです。わたしは、一般化されたものを知りません。
6n-1, 合成数の数列
6, 11, 13, 16, 20, 21, 24, 26, 27, 31, 34, 35, 36, 37, 41, 46, 48, 50
6n+1, 合成数の数列
4, 8, 9, 14, 15, 19, 20, 22, 24, 28, 29, 31, 34, 36, 39, 41, 42, 43, 44, 48, 49, 50
上記の数列を生成したついでに、6n±1篩として素数生成プログラムを書いてきました。
上記の数列を生成時、エラトステネスの篩を使用して配列に素数を登録しています。
エラトステネスの篩は、しばっち様に教えて頂いたプログラムを実装しています。
Re: 6n±1篩 #4264
プログラム(6n±1篩)を書き始めてから一月が過ぎました。毎日プログラムを眺めているわけではないのですが、
LET S=6
DO WHILE S<IP(1E8/6)
LET S=S+5
LET C=C+1
LOOP
PRINT S;":";c
PRINT IP(1E8/6)/5
END
私の能力レベルでは、これは、ひらめき
BASIC Acc の配列上限の限界克服法は、サンダラムの篩と同じ方法です。
上記二点の変更で、飛躍的に処理速度が向上しました。
計算時間 BASIC Acc を使用
1億～100億
AMD 2.40 GHz 1005.68(16分45.68秒)
1億～1000億
AMD 2.40 GHz 10533.047(2時間55分33.047秒)10533/999=10.54
i7@4.5GHz 1912.56(31分52.56秒)1913/999=1.9149
1億～1兆
i7@4.5GHz 20411.903 (5時間40分12秒)
Intel Core i7 4790K(Win7 Pro 32bt)
AMD Athlon(tm) 64 Processor 3800+ (Win 8.1 32bt)
i7@4.5GHzの場合サンダラムの篩の約4倍速い。
-------------------------------------------------
DECLARE EXTERNAL FUNCTION cprime
DECLARE EXTERNAL SUB sqprime
OPTION ARITHMETIC NATIVE !2進モード
LET t0=TIME
LET k6=1E8
CALL sqprime(k6)
LET SA=2 !MIN 2 (1億to)
LET SB=1000 !MAX 10200 (1兆200億)
LET ki=SB
DIM Ba(ki)
OPEN #1:NAME "prime_n.txt",ACCESS INPUT
FOR i=1 TO ki
INPUT #1: BA(i)
NEXT i
CLOSE #1
DIM Bb(ki)
OPEN #2:NAME "prime_pi.txt",ACCESS INPUT
FOR i=1 TO ki
INPUT #2: BB(i)
NEXT i
CLOSE #2
!LET c1=4118054813 !1E11 (1000億)
!LET c1=450708777 !1E10 (100億)
!LET c1=50847534 !1E9 (10億)
LET c1=5761455 !1E8 (1億)
LET k6=SA*1E8-1E8
FOR n=SA TO SB
LET k=Ba(n) !SQR(1E8)
LET k2=Bb(n) !π(x)
LET k6=k6+1E8
LET C=cprime(k,k2,k6)
LET c1=c1+c
PRINT n;":";c1;":";c
NEXT n
LET TM=TIME-t0
PRINT USING"######." & REPEAT$("#",2):TM;
PRINT "秒"
END
EXTERNAL FUNCTION cprime (k,k2,k6)
OPTION ARITHMETIC NATIVE
!エラトステネスの篩
DIM P(k)
DIM A(k2) !素数
SUB prime(k)
MAT P=ZER
MAT A=ZER
LET H1=1
LET A(H1)=2
LET H1=2
FOR I=3 TO SQR(K) STEP 2
IF P(I)=0 THEN
FOR J=I*I TO K STEP I
LET P(J)=1
NEXT J
END IF
NEXT I
FOR I=3 TO K STEP 2
IF P(I)=0 THEN
LET A(H1)=I
LET H1=H1+1
END IF
NEXT I
END SUB
CALL prime(k)
LET k4=k6-1E8
LET B6=IP(k6/30)+1
LET U=IP(k6/6)+1 !k5=u
LET W=IP(k4/6)-1 !k3=w
LET M7=W-2
DIM x(W-M7 TO U-M7)
DIM y(W-M7 TO U-M7)
MAT x = ZER !(6*n-1)
MAT y = ZER !(6*n+1)
FOR n=3 TO k2
LET P6=A(n)
LET G1=IP(W/P6)-2
IF MOD(P6+1,6)=0 THEN !(6*n-1)
LET r6=(P6+1)/6
CALL nap(P6,G1,r6)
END IF
IF MOD(P6-1,6)=0 THEN
LET r6=(P6-1)/6
CALL nap(P6,G1,r6*(-1))
END IF
IF MOD(P6+1,6)=0 THEN !(6*n+1)
LET r6=(P6+1)/6
CALL napi(P6,G1,r6*(-1))
END IF
IF MOD(P6-1,6)=0 THEN
LET r6=(P6-1)/6
CALL napi(P6,G1,r6)
END IF
NEXT n
SUB nap(z,t1,r6) !(6*n-1)
FOR i=t1 TO B6
IF z*i+r6<W THEN GOTO 50
IF z*i+r6>U THEN EXIT SUB
LET x(z*i+r6-M7)=1
50 NEXT i
END SUB
SUB napi(z,t1,r6) !(6*n+1)
FOR i=t1 TO B6
IF z*i+r6<W THEN GOTO 60
IF z*i+r6>U THEN EXIT SUB
LET y(z*i+r6-M7)=1
60 NEXT i
END SUB
LET Cc=0
FOR n=W-M7 TO U-M7
LET ST=n+M7
IF x(n)=1 THEN GOTO 100
IF 6*ST-1>k4 AND 6*ST-1<k6 THEN
LET cc=cc+1
IF 1=cc THEN PRINT 6*ST-1;" :";cc;" 6*n-1"
END IF
100 IF y(n)=1 THEN GOTO 200
IF 6*ST+1>k4 AND 6*ST+1<k6 THEN
LET cc=cc+1
IF 1=cc THEN PRINT 6*ST+1;" :";cc;" 6*n+1"
END IF
200 NEXT n
LET cprime=cc
END FUNCTION
EXTERNAL SUB sqprime(x)
OPTION ARITHMETIC NATIVE
!エラトステネスの篩
LET k=1009997
DIM P(k)
DIM A(79251) !素数
MAT P=ZER
MAT A=ZER
LET H1=1
LET A(H1)=2
LET H1=2
FOR I=3 TO SQR(K) STEP 2
IF P(I)=0 THEN
FOR J=I*I TO K STEP I
LET P(J)=1
NEXT J
END IF
NEXT I
FOR I=3 TO K STEP 2
IF P(I)=0 THEN
LET A(H1)=I
LET H1=H1+1
END IF
NEXT I
OPEN #1:NAME "prime_n.txt",RECTYPE INTERNAL
ERASE #1
OPEN #2:NAME "prime_pi.txt",RECTYPE INTERNAL
ERASE #2
LET x1=SQR(x)
FOR i=1229 TO 79251
LET v=A(i)
IF x1=<v THEN
LET c=c+1
LET x=x+1E8
LET x1=SQR(x)
WRITE #1:v
WRITE #2:i
!PRINT c;":";v;":";i
END IF
NEXT i
CLOSE #1
CLOSE #2
END SUB
-------------------------------------
いつも分かりにくい。プログラム解説
基本的に1億間の素数をカウントします。
LET SA=2 !MIN 2 (1億to) 最小値
LET SB=1000 !MAX 10200 (1兆200億) 最大値
-------------------------------
LET k6=1E8
CALL sqprime(k6)まだ、未完成です。SQR(x)よりも余分に出ます。
1009997(79251) 1億単位の素数と個数リストを生成します。
6n±1篩のnを篩う素数は素数^2よりも少ない自然数になる?
5から100までの素数23個で、6n±1の合成数列1733まで篩えます。
1733*6-1=10397
SQR(10000)の範囲は暫定です。この6n±1篩の数値法則の検証目的でもあります。確認は1兆まで
現状でのプログラムは、手探り状態です。このSQR(10000)の範囲で重複する、6n±1の合成数列は約4倍見当でした。
ここの重複部分をどのようにするか? 今後の課題です。別の方法もある?
sqprime(k6)はリストを生成しますので、数値以下の場合は、何度も実行する必要なし。
1E14は、1E7^2なので、15桁の限界が100兆から1000兆の辺りにありますが、未確認です。
1000兆は、素数リスト生成可能領域なので、わたしの今後の課題です。
BASIC Accが2進モードで60桁になったらいいなー。
☆その他の手動入力
!LET c1=4118054813 !1E11 (1000億)
!LET c1=450708777 !1E10 (100億)
!LET c1=50847534 !1E9 (10億)
LET c1=5761455 !1E8 (1億)
プログラムを訂正しました。
100兆からの素数Counter 6n±1篩 BASIC Acc
計算結果
100000007 : 1 6*n-1 ☆ 1億から開始して最初の素数
2 : 11078937 : 5317482 ☆ 2億未満までの素数をカウント
200000033 : 1 6*n-1
3 : 16252325 : 5173388
300000007 : 1 6*n+1
プログラムに出力部を書く足せば素数リストも出力できます。精度は未確認です。
素数の個数は1兆まで確認しました。
最初のプログラム
6n±1篩 投稿者：たろさ 投稿日：2016年12月22日(木)05時53分9秒
#4207
http://blogs.yahoo.co.jp/donald_stinger
|
|
|
|
投稿者：たろさ
投稿日：2017年 2月19日(日)20時00分39秒
|
|
|
|
> No.4273[元記事へ]
追記 プログラムを修正しました。
1億～100億
AMD 2.40 GHz 998.74(16分38.74秒)BASIC Accを使用
990.81秒(16分30.81秒)BASICAcc0981 (win32 Lazarus + fpc 3.0.0)
精度確認中です。 現在10兆まで精度確認済
100000 : 346065536839 : 3339819
346065536839 !1E13 (10兆)
計算時間：2日と11時間0分39秒
-------------------------------------------
1 : 5761455 : 0
20000 : 73301896139 : 3531519
2000000000000 : 73301896139
11時間23分11秒 (Intel Core i7 4790K@4.5GHz)
-------------------------------------------
今回は、15桁までテストしました。
9999989E+8 TO 10000000E+8(1E+15)
9999989 : 0 : 0
9999990 : 2894210 : 2894210
9999991 : 5789728 : 2895518
9999992 : 8685927 : 2896199
9999993 : 11579320 : 2893393
9999994 : 14475753 : 2896433
9999995 : 17369832 : 2894079
9999996 : 20266976 : 2897144
9999997 : 23161706 : 2894730
9999998 : 26058360 : 2896654
9999999 : 28952608 : 2894248
10000000 : 31848339 : 2895731
AMD 2.40 GHz 229.88秒(3分49.88秒)旧プログラム
BASICAcc0981 (win32 Lazarus + fpc 2.6.4)176.88秒 (2分56.88秒)
BASICAcc0981 (win32 Lazarus + fpc 3.0.0)155.10秒 (2分35.1秒)
---------------------------------------------------------------
精度確認は、しばっち様のプログラムで確かめました。
エラトステネスの篩
ERATOS2_DLL(4スレッドで実行)68.47秒 (Intel Core i7 4790K@4.5GHz)
#4101
エラトステネスの篩 eratos_dll
#3985
#3986
#3987
-------------------------------------------------
!6n±1篩 prime number Counting !Ver.2.13
!素数計数関数（英: Prime-counting function）π(x)
DECLARE EXTERNAL FUNCTION cprime
DECLARE EXTERNAL SUB sqprime
OPTION ARITHMETIC NATIVE
!PRINT DATE$;"/"; TIME$
LET t0=TIME
LET k=31622809
LET k2=1951961
!エラトステネスの篩
LET Fu=5639
DIM P(Fu)
DIM A(740) !素数
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 GT(k2)
LET Q=6
LET k7=k !篩の計算範囲
LET k5=IP(k7/Q)+1
DIM Au(k5),Av(k5)
MAT Au = ZER !(6*n-1)
MAT Av = ZER !(6*n+1)
FOR n=3 TO 740
LET Pu=A(n)
IF MOD(Pu+1,Q)=0 THEN !(6*n-1)
LET ru=(Pu+1)/Q
FOR i=1 TO k5
IF Pu*i+ru>k5 THEN EXIT FOR
LET Au(Pu*i+ru)=1
NEXT i
END IF
IF MOD(Pu-1,Q)=0 THEN
LET ru=(Pu-1)/Q
FOR i=1 TO k5
IF Pu*i-ru>k5 THEN EXIT FOR
LET Au(Pu*i-ru)=1
NEXT i
END IF
IF MOD(Pu+1,Q)=0 THEN !(6*n+1)
LET ru=(Pu+1)/Q
FOR i=1 TO k5
IF Pu*i-ru>k5 THEN EXIT FOR
LET Av(Pu*i-ru)=1
NEXT i
END IF
IF MOD(Pu-1,Q)=0 THEN
LET ru=(Pu-1)/Q
FOR i=1 TO k5
IF Pu*i+ru>k5 THEN EXIT FOR
LET Av(Pu*i+ru)=1
NEXT i
END IF
NEXT n
LET GT(1)=2
LET GT(2)=3
LET cz=2
FOR n=1 TO k5
IF 6*n-1>k7 THEN GOTO 100
IF Au(n)=0 THEN
LET cz=cz+1
LET GT(cz)=6*n-1
END IF
100 IF 6*n+1>k7 THEN EXIT FOR
IF Av(n)=0 THEN
LET cz=cz+1
LET GT(cz)=6*n+1
END IF
NEXT n
LET SA=2 !11092 ,634831,190331, 70141,26451,754901,250000, 9999990,4000001,MIN 2 (1億to)
LET SB=100 !11093 ,634840,190340, 70150,26461,754910,250010,10000000,4100000,MAX 1E15-1
CALL sqprime(SA,SB,k2,GT)
DIM Ba(SA TO SB)
DIM Bb(SA TO SB)
OPEN #1:NAME "prime_n1.txt",ACCESS INPUT
OPEN #2:NAME "prime_pi1.txt",ACCESS INPUT
FOR i=SA TO SB
INPUT #1: BA(i)
INPUT #2: BB(i)
NEXT i
CLOSE #1
CLOSE #2
! !SA=数値入力+1
!LET c1=12273824155491 !400兆
!LET c1=12270876688507 !3999009億
!LET c1=11678533758926 !380兆
!LET c1=3204941750802 !1E14 (100兆)1000001億
!LET c1=2894232250783 !9E13 (90兆) 900001億
!LET c1=2441437469602 !754900 : 2441437469602
!LET c1=1638923764567 !5E13 (50兆) 500001億
!LET c1=2064698005481 !634830 : 2064698005481
!LET c1=838534798134 !249999 : 838534798134
!LET c1=644296877669 !190330 : 644296877669
!LET c1=346065536839 !1E13 (10兆) 100001億
!LET c1=245751018905 !70140 : 245751018905
!LET c1=95957135910 !26450 : 95957135910
!LET c1=41548852740 !11091 : 41548852740
!LET c1=37607912018 !1E12 (1兆) 10001億
!LET c1=4118054813 !1E11 (1000億) 1001億
!LET c1=455052511 !1E10 (100億) 101億
!LET c1=450708777 !9.9E9 (99億) 100億
!LET c1=50847534 !1E9 (10億) 11億
LET c1=5761455 !1E8 (1億) 2億
PRINT SA-1;":";c1;":";c
LET k6=SA*1E8-1E8
FOR n=SA TO SB
!LET t1=TIME
LET k=Ba(n) !SQR(1E8)
LET k2=Bb(n) !π(x)
LET k6=k6+1E8
LET C=cprime(k,k2,k6,GT)
LET c1=c1+c
PRINT n;":";c1;":";c
!LET TM=TIME-t1
!PRINT USING"######." & REPEAT$("#",2):TM;
!PRINT "秒"
NEXT n
LET TM=TIME-t0
PRINT USING"######." & REPEAT$("#",2):TM;
PRINT "秒"
!PRINT DATE$;"/"; TIME$
END
EXTERNAL FUNCTION cprime (k,k2,k6,GT())
OPTION ARITHMETIC NATIVE
LET k4=k6-1E8
LET B6=IP(k6/30)+1
LET U=IP(k6/6)+1
LET W=IP(k4/6)-1
LET M7=W-2
DIM x(W-M7 TO U-M7)
DIM y(W-M7 TO U-M7)
MAT x = ZER !(6*n-1)
MAT y = ZER !(6*n+1)
FOR n=3 TO k2
LET P6=GT(n)
LET G1=IP(W/P6)!-2
IF MOD(P6+1,6)=0 THEN !(6*n-1)
LET r6=(P6+1)/6
FOR i=G1 TO B6
IF P6*i+r6<W THEN GOTO 140
IF P6*i+r6>U THEN EXIT FOR
LET x(P6*i+r6-M7)=1
140 NEXT i
END IF
IF MOD(P6-1,6)=0 THEN
LET r6=(P6-1)/6
FOR i=G1 TO B6
IF P6*i-r6<W THEN GOTO 150
IF P6*i-r6>U THEN EXIT FOR
LET x(P6*i-r6-M7)=1
150 NEXT i
END IF
IF MOD(P6+1,6)=0 THEN !(6*n+1)
LET r6=(P6+1)/6
FOR i=G1 TO B6
IF P6*i-r6<W THEN GOTO 160
IF P6*i-r6>U THEN EXIT FOR
LET y(P6*i-r6-M7)=1
160 NEXT i
END IF
IF MOD(P6-1,6)=0 THEN
LET r6=(P6-1)/6
FOR i=G1 TO B6
IF P6*i+r6<W THEN GOTO 170
IF P6*i+r6>U THEN EXIT FOR
LET y(P6*i+r6-M7)=1
170 NEXT i
END IF
NEXT n
LET Cc=0
FOR n=W-M7 TO U-M7
LET ST=n+M7
IF x(n)=1 THEN GOTO 180
IF 6*ST-1>k4 AND 6*ST-1<k6 THEN
LET cc=cc+1
! IF 1=cc THEN PRINT 6*ST-1;" :";cc;" 6*n-1"
! IF 3102679=cc OR 3099882=cc THEN PRINT 6*ST-1;" :";cc;" 6*n-1"
END IF
180 IF y(n)=1 THEN GOTO 200
IF 6*ST+1>k4 AND 6*ST+1<k6 THEN
LET cc=cc+1
! IF 1=cc THEN PRINT 6*ST+1;" :";cc;" 6*n+1"
! IF 3102679=cc OR 3099882=cc THEN PRINT 6*ST+1;" :";cc;" 6*n+1"
END IF
200 NEXT n
LET cprime=cc
END FUNCTION
EXTERNAL SUB sqprime(SA,SB,k2,GT())
OPTION ARITHMETIC NATIVE
OPEN #1:NAME "prime_n1.txt",RECTYPE INTERNAL
ERASE #1
OPEN #2:NAME "prime_pi1.txt",RECTYPE INTERNAL
ERASE #2
LET x=1E8
FOR i=1 TO k2-1
LET v=GT(i)^2
LET vi=GT(i+1)^2
LET vv=vi-v
IF x =< v THEN
CALL sqrp(1)
IF vv>=1E8 THEN
LET vt=IP(vv/1E8)
CALL sqrp(vt)
END IF
END IF
NEXT i
SUB sqrp(vo)
FOR ns=1 TO vo
LET C=C+1
LET X=X+1E8
IF C>=SA AND SB>=C THEN
WRITE #1:GT(i)
WRITE #2:i
END IF
NEXT ns
END SUB
CLOSE #1
CLOSE #2
END SUB
-------------------------------------------------
変数
LET SA=2 !MIN 2 (1億to) 最小値 2から1億も分岐を付ければ計算可能です。6n±1なので5からです。
LET SB=20000!MAX 1E15-1 15桁まで
LET c1=5761455 !1E8 (1億)
----------------------------------------
EXTERNAL SUB sqprime は、計算の無駄を省きました。
素数= pn
(pn+1)^2-pn^2 の性質を探究中です。
http://blogs.yahoo.co.jp/donald_stinger
|
|
|
戻る