|
投稿者：山中和義
投稿日：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
|
|
|
|
投稿者：山中和義
投稿日：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 )
:
:
|
|
|
|
投稿者：山中和義
投稿日：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 通り
|
|
|
|
投稿者：山中和義
投稿日：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
※サブルーチン部分は省略する
|
|
|
|
投稿者：山中和義
投稿日：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
※サブルーチン部分は省略する
|
|
|
|
投稿者：山中和義
投稿日：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
※サブルーチン部分は省略する
|
|
|
|
投稿者：たろさ
投稿日：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
|
|
|
|
投稿者：しばっち
投稿日：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までの対応のようです)
|
|
|
|
投稿者：白石 和夫
投稿日：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
|
|
|
|
投稿者：たろさ
投稿日：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
|
|
|
戻る
