十進BASIC第2掲示板過去ログ

素数関連の関数

投稿者：山中和義投稿日：2016年 3月12日(土)20時29分50秒

素数判定の関数PrimeQをベースにした関連の関数を実装してみました。 !n番目の素数　※UBASIC prm(n)関数 FOR N=2 TO 168 PRINT N; Prime(N); N*LOG(N)+N*LOG(LOG(N)) !n番目の素数の近似 NEXT N PRINT !nより大きい素数の最小のもの　※UBASIC nxtprm(x)関数 FOR N=0 TO 100 PRINT N; NextPrime(N) NEXT N PRINT !n以下の素数の数 LET N=1000 PRINT "素数は";PrimePi(N);"個です。" PRINT N/LOG(N) !ガウス・アダマール・ピュイサンの素数定理による近似 PRINT N/(LOG(N)-1.08366) !ルジャンドルの近似 END !試行割算法（6k±1） 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まで検証する !+1,+3,+5は２の倍数（偶数）、+1,+4は３の倍数、+5は５の倍数 IF MOD(n,k)=0 THEN EXIT FUNCTION !5,11,17,23,29,35,41,47,… IF MOD(n,k+2)=0 THEN EXIT FUNCTION !7,13,19,25,31,37,43,49,… LET k=k+6 LOOP LET PrimeQ=1 !最後まで割り切れなければ、素数である END IF END FUNCTION EXTERNAL FUNCTION Prime(n) !n番目の素数 IF n<1 OR n<>INT(n) THEN EXIT FUNCTION !引数を確認する IF n=1 THEN LET Prime=2 !素数2 ELSEIF n=2 THEN LET Prime=3 !素数3 ELSE LET c=2 !見つけた個数 LET k=5 DO !+1,+3,+5は２の倍数（偶数）、+1,+4は３の倍数、+5は５の倍数 IF PrimeQ(k)<>0 THEN LET Prime=k LET c=c+1 !5,11,17,23,29,… が素数なら IF c=n THEN EXIT DO END IF IF PrimeQ(k+2)<>0 THEN LET Prime=k+2 LET c=c+1 !7,13,19,25,31,… IF c=n THEN EXIT DO END IF LET k=k+6 !次へ LOOP END IF END FUNCTION EXTERNAL SUB PrimeList(n,p()) !n個の素数列を返す IF n<1 OR n<>INT(n) THEN EXIT SUB !引数を確認する LET c=1 !見つけた個数 LET p(c)=2 !2は素数 IF n=1 THEN EXIT SUB LET c=2 LET p(c)=3 !3は素数 IF n=2 THEN EXIT SUB LET k=5 DO IF PrimeQ(k)<>0 THEN !5,11,17,23,29,… が素数なら LET c=c+1 LET p(c)=k IF c=n THEN EXIT DO END IF IF PrimeQ(k+2)<>0 THEN !7,13,19,25,31,… LET c=c+1 LET p(c)=k+2 IF c=n THEN EXIT DO END IF !+1,+3,+5は２の倍数（偶数）、+1,+4は３の倍数、+5は５の倍数 LET k=k+6 !次へ LOOP END SUB EXTERNAL FUNCTION PrimePi(n) !n以下の素数の数 IF n<2 THEN LET PrimePi=0 !見つけた個数 ELSEIF n<3 THEN LET PrimePi=1 !見つけた個数 ELSEIF n<5 THEN LET PrimePi=2 !見つけた個数 ELSE LET c=2 !見つけた個数 LET k=5 DO !+1,+3,+5は２の倍数（偶数）、+1,+4は３の倍数、+5は５の倍数 IF k>n THEN EXIT DO IF PrimeQ(k)<>0 THEN LET c=c+1 !5,11,17,23,29,… が素数なら IF k+2>n THEN EXIT DO IF PrimeQ(k+2)<>0 THEN LET c=c+1 !7,13,19,25,31,… LET k=k+6 !次へ LOOP LET PrimePi=c END IF END FUNCTION EXTERNAL FUNCTION NextPrime(n) !nより大きい素数の最小のもの IF n<2 THEN LET NextPrime=2 ELSEIF n<3 THEN LET NextPrime=3 ELSE LET k=CEIL(n/6)*6 -1 !6m±1 DO IF k>n THEN !見つかるまで LET NextPrime=k IF PrimeQ(k)<>0 THEN EXIT DO END IF IF k+2>n THEN !見つかるまで LET NextPrime=k+2 IF PrimeQ(k+2)<>0 THEN EXIT DO END IF LET k=k+6 !次へ LOOP END IF END FUNCTION 以下、サンプルです。ただし、サブルーチン部分は省略します。 ●素数階乗 FOR N=2 TO 30 !n# LET S=1 LET i=2 !素数2 DO WHILE i<=N LET S=S*i LET i=NextPrime(i) !次の素数 LOOP PRINT N;S NEXT N END ●区間[a,b]の素数の個数 LET A=1*10^7 ![a,b] LET B=2*10^7 LET C=0 IF PrimeQ(A)<>0 THEN LET C=C+1 !aが素数 LET i=NextPrime(A) DO WHILE i<=B LET C=C+1 LET i=NextPrime(i) LOOP PRINT C;"個" END ●区間篩い法 LET L=1*10^7 ![L,H]、1≦L≦H ※10^7以下 LET H=2*10^7 DIM F(L TO H) !篩い MAT F=ZER IF L=1 THEN LET F(1)=1 !1は素数でない FOR J=2^2 TO H STEP 2 !素数2 IF J>=L THEN LET F(J)=1 NEXT J LET i=3 !奇素数 DO WHILE i^2<=H !√hまで FOR J=i^2 TO H STEP i*2 IF J>=L THEN LET F(J)=1 NEXT J LET i=NextPrime(i) LOOP LET C=0 FOR i=L TO H IF F(i)=0 THEN LET C=C+1 !!!PRINT i !debug END IF NEXT i PRINT C;"個" END ●包除原理（個数定理） LET N=10^6 !1からnまで LET K=PrimePi(SQR(N)) !1から√nまでの素数を用意する DIM P(K) !k個の素数列 CALL PrimeList(K,P) LET NC=0 CALL try(1,1,1, N,K,P, NC) PRINT NC+K-1;"個" END EXTERNAL SUB try(i,S,B, N,K,P(), C) !m個の素数の積による[n/p] IF i=K+1 THEN LET C=C+B*INT(N/S) ![n/s] ELSE CALL try(i+1,S,B, N,K,P, C) !ビットパターンが0のとき LET T=S*P(i) !ビットパターンが1のとき IF T<=N THEN CALL try(i+1,T,-B, N,K,P, C) !s＞nのとき、[n/s]=0 と昇順列なので、これ以降は0となる END IF END SUB

エラトステネスの篩い法

投稿者：山中和義投稿日：2016年 3月20日(日)00時25分20秒

> No.3999[元記事へ] 技術メモ i:1,2,3,4, 5, 6, 7, 8, 9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,… t:3,5,7,9,11,13,15,17,19,21,23,25,27,29,31,33,35,37,39,41,43,45,47,49,51,53,55,57,59,…　奇数 3 * * * * * * * * * ↑(t+1)*i から　t ずつ 5 * * * * 7 * t=2i+1で素数のとき、t^2=(2i+1)^2=4i^2+4i+1=2{((2i+1)+1)i}+1　∴(t+1)iの位置となる LET t0=TIME LET N=10^7 !１からＮまでの素数 PUBLIC NUMERIC B_bits !C言語 bit F(B); の意 LET B_bits=32 LET B=INT((N-1)/2) !ビット数 LET M=CEIL(B/B_bits) !ビットを正の整数へ PRINT B; M !debug DIM F(M) !素数候補として、奇数 3,5,7,9,11,13,15,… を用意する MAT F=ZER LET c=0 IF N>1 THEN LET c=1 !２は素数 FOR i=1 TO B !３以上 CALL B_bit(i,F, rc) !i番目が素数なら IF rc=0 THEN LET c=c+1 LET t=2*i+1 !自然数に換算する !!PRINT c; t FOR k=(t+1)*i TO B STEP t !その２乗以上の奇数倍を削除する　※偶数倍は偶数となる CALL B_bitset(k,F) !F(k)=1 NEXT k END IF NEXT i PRINT c !個数 PRINT "計算時間=";TIME-t0 END !ビット配列　３２ビットを正の整数へ EXTERNAL SUB B_bit(n,B(), rc) !n番目のビット値　※n,B()は整数 LET x=B(INT(n/B_bits)+1) LET d=2^MOD(n,B_bits) LET rc=MOD(INT(x/d),2) END SUB EXTERNAL SUB B_bitset(n,B()) !n番目のビットを１にする　※n,B()は非負整数 LET t=INT(n/B_bits)+1 LET x=B(t) !n桁 LET d=2^MOD(n,B_bits) LET B(t)=(INT(x/d/2)*2+1)*d+MOD(x,d) !大きい桁＋１＋小さい桁 END SUB

Re: 素数関連の関数

投稿者：山中和義投稿日：2016年 3月23日(水)10時54分44秒

> No.3999[元記事へ] 問題素数しりとり 100までの素数（25個）を文字と見て、しりとりをすることを考える。例えば 2 → 23 → 31 → 11 → 19 → 97 → 73 → 37 → 79 なら、9連のしりとりとなる。さて、最長何連のしりとりが可能か？考察次につなげられる素数は、 FOR M=1 TO 9 PRINT M;"→"; LET P=2 DO WHILE P<=100 !100以下の素数 ※ LET T=P DO WHILE T>=10 LET T=INT(T/10) LOOP IF T=M THEN PRINT P; LET P=NextPrime(P) LOOP PRINT NEXT M END ※サブルーチン部分は省略する。(他も同様) として、 1 → 11 13 17 19 2 → 2 23 29 3 → 3 31 37 4 → 41 43 47 5 → 5 53 59 6 → 61 67 7 → 7 71 73 79 8 → 83 89 9 → 97 を得る。状態遷移表を作成すると、 DIM A(9,9) !A(m,n): mで始まってnで終わるものの個数 MAT A=ZER LET P=2 DO WHILE P<=100 !100以下の素数 ※ LET T=P DO WHILE T>=10 LET T=INT(T/10) LOOP LET W=MOD(P,10) LET A(T,W)=A(T,W)+1 LET P=NextPrime(P) LOOP MAT PRINT A; !debug END として、 1 0 1 0 0 0 1 0 1 0 1 1 0 0 0 0 0 1 1 0 1 0 0 0 1 0 0 1 0 1 0 0 0 1 0 0 0 0 1 0 1 0 0 0 1 1 0 0 0 0 0 1 0 0 1 0 1 0 0 0 1 0 1 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 を得る。異なる素数列に注意して、この表を使ってシミュレーションする。実行結果 1 1 3 1 7 1 9 7 3 3 ( 9 ) は、 1-1 1-3 3-1 1-7 7-1 1-9 9-7 7-3 3-3 を表す。具体的には、1-1に属する素数11を選択する。複数存在する場合は任意に1つ選べばよい。他も同様に、 11 13 31 17 71 19 97 73 3 となる。（終わり） LET t0=TIME DIM A(9,9) !A(m,n): mで始まってnで終わるものの個数 MAT A=ZER LET P=2 DO WHILE P<=100 !100以下の素数 ※ LET T=P DO WHILE T>=10 LET T=INT(T/10) LOOP LET W=MOD(P,10) LET A(T,W)=A(T,W)+1 LET P=NextPrime(P) LOOP MAT PRINT A; !debug DIM X(0 TO 100) !※ FOR M=1 TO 9 FOR N=1 TO 9 PRINT "M=";M;"N=";N LET W=A(M,N) !save it IF W>0 THEN LET A(M,N)=A(M,N)-1 !1つ使う LET X(0)=M !m-n LET X(1)=N CALL try(2,N,A,X) LET A(M,N)=W !restore it END IF NEXT N NEXT M PRINT TIME-t0 END EXTERNAL SUB try(P,M,A(,),X()) !バックトラック法でシミュレーションする LET FLG=-1 FOR N=1 TO 9 !遷移先を選ぶ LET W=A(M,N) !save it IF W>0 THEN LET A(M,N)=A(M,N)-1 LET X(P)=N CALL try(P+1,N,A,X) !次へ LET A(M,N)=W !restore it LET FLG=0 END IF NEXT N IF FLG<>0 THEN !mからの遷移がすべて不可能なとき FOR i=0 TO P-1 !結果を表示する PRINT X(i); NEXT i PRINT " (";P-1;")" END IF END SUB 実行結果 1 0 1 0 0 0 1 0 1 0 1 1 0 0 0 0 0 1 1 0 1 0 0 0 1 0 0 1 0 1 0 0 0 1 0 0 0 0 1 0 1 0 0 0 1 1 0 0 0 0 0 1 0 0 1 0 1 0 0 0 1 0 1 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 M= 1 N= 1 1 1 3 1 7 1 9 7 3 3 7 7 9 ( 12 ) 1 1 3 1 7 1 9 7 3 3 7 9 ( 11 ) 1 1 3 1 7 1 9 7 3 7 7 9 ( 11 ) 1 1 3 1 7 1 9 7 3 7 9 ( 10 ) 1 1 3 1 7 1 9 7 7 3 3 7 9 ( 12 ) 1 1 3 1 7 1 9 7 7 3 7 9 ( 11 ) 1 1 3 1 7 1 9 7 7 9 ( 9 ) 1 1 3 1 7 1 9 7 9 ( 8 ) 1 1 3 1 7 3 3 7 1 9 7 7 9 ( 12 ) 1 1 3 1 7 3 3 7 1 9 7 9 ( 11 ) 1 1 3 1 7 3 3 7 7 1 9 7 9 ( 12 ) 1 1 3 1 7 3 3 7 7 9 7 1 9 ( 12 ) 1 1 3 1 7 3 3 7 9 7 1 9 ( 11 ) 1 1 3 1 7 3 3 7 9 7 7 1 9 ( 12 ) 1 1 3 1 7 3 7 1 9 7 7 9 ( 11 ) 1 1 3 1 7 3 7 1 9 7 9 ( 10 ) 1 1 3 1 7 3 7 7 1 9 7 9 ( 11 ) 1 1 3 1 7 3 7 7 9 7 1 9 ( 11 ) : :

Re: 素数関連の関数

投稿者：山中和義投稿日：2016年 3月24日(木)07時57分53秒

> No.4011[元記事へ] 続き状態遷移表 1 0 1 0 0 0 1 0 1 0 1 1 0 0 0 0 0 1 1 0 1 0 0 0 1 0 0 1 0 1 0 0 0 1 0 0 0 0 1 0 1 0 0 0 1 1 0 0 0 0 0 1 0 0 1 0 1 0 0 0 1 0 1 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 を見ると、100までの素数の個数は1個ずつなので、 11 0 13 0 0 0 17 0 19 0 2 23 0 0 0 0 0 29 31 0 3 0 0 0 37 0 0 41 0 43 0 0 0 47 0 0 0 0 53 0 5 0 0 0 59 61 0 0 0 0 0 67 0 0 71 0 73 0 0 0 7 0 79 0 0 83 0 0 0 0 0 89 0 0 0 0 0 0 97 0 0 のように素数を記録すると、結果を表示するのに都合がよい。 13連の場合の数を求めると、 LET t0=TIME DIM A(9,9) !A(m,n): mで始まってnで終わるものの個数 MAT A=ZER LET P=2 DO WHILE P<=100 !100以下の素数 ※ LET T=P DO WHILE T>=10 LET T=INT(T/10) LOOP LET W=MOD(P,10) LET A(T,W)=P LET P=NextPrime(P) LOOP MAT PRINT A; !debug PUBLIC NUMERIC C !場合の数 LET C=0 DIM X(100) !※ FOR M=1 TO 9 !開始位置を決める FOR N=1 TO 9 PRINT "M=";M;"N=";N LET W=A(M,N) !save it IF W>0 THEN LET A(M,N)=0 !1つ使う LET X(1)=W CALL try(2,N,A,X) LET A(M,N)=W !restore it END IF NEXT N NEXT M PRINT C;"通り" PRINT TIME-t0 END EXTERNAL SUB try(P,M,A(,),X()) !バックトラック法でシミュレーションする FOR N=1 TO 9 !遷移先を選ぶ LET W=A(M,N) !save it IF W>0 THEN LET A(M,N)=0 LET X(P)=W CALL try(P+1,N,A,X) !次へ LET A(M,N)=W !restore it END IF NEXT N IF P-1>=13 THEN FOR i=1 TO P-1 !結果を表示する PRINT X(i); NEXT i PRINT " (";P-1;")" LET C=C+1 END IF END SUB 実行結果 11 0 13 0 0 0 17 0 19 0 2 23 0 0 0 0 0 29 31 0 3 0 0 0 37 0 0 41 0 43 0 0 0 47 0 0 0 0 53 0 5 0 0 0 59 61 0 0 0 0 0 67 0 0 71 0 73 0 0 0 7 0 79 0 0 83 0 0 0 0 0 89 0 0 0 0 0 0 97 0 0 M= 1 N= 1 M= 1 N= 2 M= 1 N= 3 M= 1 N= 4 M= 1 N= 5 M= 1 N= 6 M= 1 N= 7 M= 1 N= 8 M= 1 N= 9 M= 2 N= 1 M= 2 N= 2 2 23 31 11 13 3 37 71 17 7 79 97 73 ( 13 ) 2 23 31 11 13 3 37 71 17 79 97 7 73 ( 13 ) 2 23 31 11 13 3 37 7 71 17 79 97 73 ( 13 ) 2 23 31 11 13 3 37 7 79 97 71 17 73 ( 13 ) 2 23 31 11 13 3 37 79 97 71 17 7 73 ( 13 ) 2 23 31 11 13 3 37 79 97 7 71 17 73 ( 13 ) : : M= 8 N= 9 M= 9 N= 1 M= 9 N= 2 M= 9 N= 3 M= 9 N= 4 M= 9 N= 5 M= 9 N= 6 M= 9 N= 7 M= 9 N= 8 M= 9 N= 9 2064 通り

Re: 素数関連の関数

投稿者：山中和義投稿日：2016年 3月26日(土)11時48分37秒

> No.3999[元記事へ] ●ゴールドバッハの予想：「4以上の全ての偶数は2つの素数の和であらわせる」 FOR N=4 TO 100 STEP 2 PRINT "n=";N LET P=2 DO WHILE P<=N/2 IF PrimeQ(N-P)<>0 THEN PRINT P;N-P LET P=NextPrime(P) LOOP NEXT N END ※サブルーチン部分は省略する ●素数ものさし +-----------------------------------+ | | | | | | | | | | | +---+-+---+---+-------+---+-------+-+ 0 2 3 5 7 11 13 17 18(=p+1) LET K=7 !個数 DIM P(0 TO K+1) !目盛り CALL PrimeList(K,P) LET P(0)=0 LET P(K+1)=P(K)+1 MAT PRINT P; DIM F(P(K+1)) !フラグ MAT F=ZER FOR i=0 TO K !2つの素数の差 FOR J=i+1 TO K+1 LET S=P(J)-P(i) LET F(S)=F(S)+1 NEXT J NEXT i MAT PRINT F; !場合の数 END ※サブルーチン部分は省略する

Re: 素数関連の関数

投稿者：山中和義投稿日：2016年 3月27日(日)16時53分53秒

> No.3999[元記事へ] ●Project Euler 35 問題　巡回素数 197は巡回素数と呼ばれる。桁を回転させたときに得られる数 197,971,719 が全て素数だからである。 100未満には巡回素数が13個ある。（具体的には、2,3,5,7,11,13,17,31,37,71,73,79,および97である） 100万未満の巡回素数はいくつあるか? LET t0=TIME LET C=4 !個数 1桁の素数は2,3,5,7 LET P=11 DO WHILE P<=10^6 !1から10^6まで LET T=P LET K=0 !k桁 DO WHILE T>0 LET W=MOD(T,10) !一の位 IF MOD(W,2)=0 OR W=5 THEN EXIT DO !偶数または5を含むなら LET T=INT(T/10) LET K=K+1 LOOP IF T=0 THEN !2桁以上は、一の位は偶数または5になることはない LET K=K-1 !（5以外の奇数のみで構成される） LET T=P !巡回させる FOR i=1 TO K LET T=MOD(T,10)*10^K+INT(T/10) !!!PRINT T; !debug IF PrimeQ(T)=0 THEN EXIT FOR NEXT i !!!PRINT "←" !debug IF i>K THEN !題意を満たす LET C=C+1 PRINT P END IF END IF LET P=NextPrime(P) !次の素数へ LOOP PRINT C;"個" PRINT TIME-t0 END ※サブルーチン部分は省略する

Re: 素数関連の関数

投稿者：山中和義投稿日：2016年 3月28日(月)23時13分55秒

> No.3999[元記事へ] ●Project Euler 37 問題　切り詰め可能素数 3797は面白い性質を持っている。まずそれ自身が素数であり、左から右に桁を除いたときに全て素数になっている (3797,797,97,7)。同様に右から左に桁を除いたときも全て素数である (3797,379,37,3)。右から切り詰めても左から切り詰めても素数になるような素数は11個しかない。総和を求めよ。注. 2,3,5,7を切り詰め可能な素数とは考えない。 LET t0=TIME LET S=0 !和 LET C=0 !個数 LET P=11 DO WHILE P<=10^6 !10^6以下 LET L=P LET K=0 !桁数 DO LET K=K+1 LET L=INT(L/10) !左 LET R=MOD(P,10^K) !右 LOOP WHILE PrimeQ(L)<>0 AND PrimeQ(R)<>0 IF L=0 THEN !題意を満たす LET C=C+1 LET S=S+P PRINT P END IF LET P=NextPrime(P) LOOP PRINT C;"個" PRINT S PRINT TIME-t0 END ※サブルーチン部分は省略するその２ LET t0=TIME PUBLIC NUMERIC C !個数 LET C=0 PUBLIC NUMERIC S !和 LET S=0 FOR K=1 TO 13 !桁数 PRINT "K=";K CALL try(1,K,i) NEXT K PRINT S PRINT C;"個" PRINT S-(2+3+5+7) !748317 PRINT TIME-t0 END EXTERNAL SUB try(M,K,P) FOR i=1 TO 9 !左を確認する（右に付加する） ※2桁以上では一の位は、1,3,7,9のみ LET T=P*10+i IF PrimeQ(T)<>0 THEN IF M=K THEN !k桁なら FOR J=1 TO K-1 !右を確認する IF PrimeQ(MOD(T,10^J))=0 THEN EXIT FOR NEXT J IF J>K-1 THEN !題意を満たす PRINT T LET C=C+1 LET S=S+T END IF ELSE CALL try(M+1,K,T) !次の桁へ END IF END IF NEXT i END SUB ※サブルーチン部分は省略する

Re: エラトステネスの篩い法

投稿者：たろさ投稿日：2016年 6月24日(金)02時03分25秒

> No.4006[元記事へ] 山中和義さんへのお返事です。お世話になっております。32Bt変換を理解したくてしばらくの間プログラムを眺めて来ました。動作報告です。 -------------------------------- OPTION ARITHMETIC NATIVE !2進モード LET t0=TIME OPEN #1:NAME "prime_10_9xv.txt",RECTYPE INTERNAL ERASE #1 WRITE #1:2 LET N=10^9 !１からＮまでの素数 PUBLIC NUMERIC B_bits !C言語 bit F(B); の意 LET B_bits=52 LET B=INT((N-1)/2) !ビット数 LET M=CEIL(B/B_bits) !ビットを正の整数へ PRINT B; M !debug DIM F(M) !素数候補として、奇数 3,5,7,9,11,13,15,… を用意する MAT F=ZER LET c=0 IF N>1 THEN LET c=1 !２は素数 FOR i=1 TO B !３以上 CALL B_bit(i,F, rc) !i番目が素数なら IF rc=0 THEN LET c=c+1 LET t=2*i+1 !自然数に換算する !PRINT c; t WRITE #1:t ! FOR k=(t+1)*i TO B STEP t !その２乗以上の奇数倍を削除する　※偶数倍は偶数となる CALL B_bitset(k,F) !F(k)=1 NEXT k END IF NEXT i CLOSE #1 PRINT c !個数 PRINT "計算時間=";TIME-t0;"秒" END !ビット配列　３２ビットを正の整数へ EXTERNAL SUB B_bit(n,B(), rc) !n番目のビット値　※n,B()は整数 OPTION ARITHMETIC NATIVE !2進モード LET x=B(INT(n/B_bits)+1) LET d=2^MOD(n,B_bits) LET rc=MOD(INT(x/d),2) END SUB ! EXTERNAL SUB B_bitset(n,B()) !n番目のビットを１にする　※n,B()は非負整数 OPTION ARITHMETIC NATIVE !2進モード LET t=INT(n/B_bits)+1 LET x=B(t) !n桁 LET d=2^MOD(n,B_bits) LET B(t)=(INT(x/d/2)*2+1)*d+MOD(x,d) !大きい桁＋１＋小さい桁 END SUB ------------------------------------ 52Btでしょうか? よく分かりません。 LET N=10^10 !１からＮまでの素数 BASIC Acc 成功しました。 N=10^11は駄目でした。 2から9999999967　ファイルサイズは5.87GB 4999999999 96153847 455052511 計算時間=-82507.465 秒ついでに、このエディター凄いと思いました。しばっちさんのエラトステネスの篩も ----------------------------------- !エラトステネスの篩 OPTION ARITHMETIC NATIVE !OPTION ARITHMETIC DECIMAL_HIGH !1000桁モード OPTION BASE 0 LET t0=TIME LET z=52 LET CC=1 OPEN #1:NAME "E:\prime\prime_1E4_xv.txt",RECTYPE INTERNAL ERASE #1 WRITE #1:2 LET N=2+1 !'3以上の奇数のみ LET M=1.E4-1 !'3以上の奇数のみ DIM X(INT(1/2*(M-N)/z)+1) FOR I=3 TO SQR(M) STEP 2 IF MOD(N,I)=0 THEN LET P=N ELSE LET P=INT(N/I+1)*I !'N以上の最小のIの倍数 END IF LET NN=(N-1)/2 IF N<=I THEN LET PP=(P-1)/2 LET L=INT((PP-NN)/z) LET K=MOD((PP-NN),z) IF BITAND(X(L),2^K)=0 THEN LET P=P+I END IF FOR J=P TO M STEP I IF MOD(J,2)=1 THEN !'奇数のみ LET JJ=(J-1)/2 LET L=INT((JJ-NN)/z) LET K=MOD((JJ-NN),z) LET X(L)=BITOR(X(L),2^K) !'印をつける END IF NEXT J NEXT I FOR J=N TO M STEP 2 LET JJ=(J-1)/2 LET L=INT((JJ-NN)/z) LET K=MOD((JJ-NN),z) IF BITAND(X(L),2^K)=0 THEN LET C=C+1 !PRINT j;":";c+cc;":";c WRITE #1:j END IF NEXT J CLOSE #1 PRINT j;":";c+cc;":";c PRINT c;"個" PRINT TIME-t0;"秒で計算しました" END --------------------------------------- 計算結果は合っているでしょうか? WRITE #1:j ファイルサイズを適当な大きさに分けて出す方法を模索しています。 64Btにするにはどうすればいいのか分かりません。

http://blogs.yahoo.co.jp/donald_stinger

Re: エラトステネスの篩い法

投稿者：しばっち投稿日：2016年 6月24日(金)20時55分5秒

> No.4091[元記事へ] たろささんへのお返事です。 > WRITE #1:j ファイルサイズを適当な大きさに分けて出す方法を模索しています。下記のようにすると、とりあえず幾つかのファイルに分割できます。 LET J=1 OPEN #1:NAME "SAMPLE_"&STR$(J)&".txt" FOR I=1 TO 10000 IF I>=1000*J THEN CLOSE #1 LET J=J+1 OPEN #1:NAME "SAMPLE_"&STR$(J)&".txt" END IF WRITE #1:I NEXT I CLOSE #1 END > 64Btにするにはどうすればいいのか分かりません。 1000桁モード、有理数モードなら可能です (1000桁モードでは3000bit以上の精度があるようです PRINT 2^3000) 2進モードや1000桁モードで下記を実行してみてください。 FOR I=1 TO 64 PRINT I;":";2^I NEXT I END 2進モードでも 2^63まで(オプション-数値の「表示桁数を多く」にチェック)は正確に表現できるようですので、63bitまでは可能です。 (但し、BITANDやBITORは53bitまでの対応のようです)

Re: エラトステネスの篩い法

投稿者：白石　和夫投稿日：2016年 6月27日(月)08時20分6秒

> No.4092[元記事へ] 2進モードの数値変数は53ビットの精度です。次のプログラムを実行してみてください。 OPTION ARITHMETIC NATIVE FOR I=1 TO 64 LET x=2^i LET y=x+1 PRINT i,y-x NEXT I END

Re: エラトステネスの篩い法

投稿者：たろさ投稿日：2016年 6月30日(木)11時54分20秒

> No.4092[元記事へ] しばっちさんへのお返事です。白石　和夫様しばっち様大変参考になりました。ありがとうございます。私事ですが、自分の子供が全国中学校体育大会の地区予選に参加するため私はサポートに追われています。遅くなってごめんなさい。動作報告です。 ----------------------------------------------- !エラトステネスの篩 !OPTION ARITHMETIC NATIVE OPTION ARITHMETIC DECIMAL_HIGH !1000桁モード OPTION BASE 0 LET t0=TIME LET z=52 LET CC=0 OPEN #1:NAME "C:\prime\prime_1E20_xv.txt",RECTYPE INTERNAL ERASE #1 !WRITE #1:2 LET N=1E20-1E6+1 !'3以上の奇数のみ LET M=1E20-1 !'3以上の奇数のみ DIM X(INT(1/2*(M-N)/z)+1) FOR I=3 TO SQR(M) STEP 2 IF MOD(N,I)=0 THEN LET P=N ELSE LET P=INT(N/I+1)*I !'N以上の最小のIの倍数 END IF LET NN=(N-1)/2 IF N<=I THEN LET PP=(P-1)/2 LET L=INT((PP-NN)/z) LET K=MOD((PP-NN),z) IF BITAND(X(L),2^K)=0 THEN LET P=P+I END IF FOR J=P TO M STEP I IF MOD(J,2)=1 THEN !'奇数のみ LET JJ=(J-1)/2 LET L=INT((JJ-NN)/z) LET K=MOD((JJ-NN),z) LET X(L)=BITOR(X(L),2^K) !'印をつける END IF NEXT J NEXT I FOR J=N TO M STEP 2 LET JJ=(J-1)/2 LET L=INT((JJ-NN)/z) LET K=MOD((JJ-NN),z) IF BITAND(X(L),2^K)=0 THEN LET C=C+1 !PRINT j;":";c+cc;":";c WRITE #1:j END IF NEXT J CLOSE #1 !PRINT j;":";c+cc;":";c PRINT c;"個" PRINT TIME-t0;"秒で計算しました" END --------------------------------------------- 計算結果 21858 個 58275.64 秒で計算しました 16時間11分15.64秒(Intel Core i7-4790K@4.5GHz) 99999999999999000017 : 2220819602560896983 : 1 99999999999999999989 : 2220819602560918840 : 21858個以前のMiller-Rabin法の計算結果と一致しています。

http://blogs.yahoo.co.jp/donald_stinger

戻る

十進BASIC 第2掲示板 過去ログ

十進BASIC 第2掲示板過去ログ