最速の素数判定アルゴリズム
■ このスレッドは過去ログ倉庫に格納されています
0001 刺身定食
NGNG文献を読んでみたのですが、サッパリわかりません。
日本で完璧に理解している方、たまたまこのスレ
見ていたら、ついでに教えてください。
マジです。
0002デフォルトの名無しさん
NGNG理解できないのならあきらめるしかないでしょう。
他人をたよるのはよくありません。
0003デフォルトの名無しさん
NGNGその「さっぱり」って言葉がやるきねぇなぁと他人に感じさせるので、
そういう聞き方では知識ナルシスト以外教えてくれないと思われ。
0004デフォルトの名無しさん
NGNG0005デフォルトの名無しさん
NGNGデータとして保持しておくものだよ。
32bit整数MAXまでの素数リストがあれば実用上何も問題なし。
http://arch.arch.kumamoto-u.ac.jp/hagane/yamanari/joho2/kougi3.html
http://www.hotwired.co.jp/news/news/technology/story/2234.html
0006デフォルトの名無しさん
NGNG神父もびっくりだな
0007デフォルトの名無しさん
NGNG日本のチームはあるの？
京都産業大学なんか大学でチームつくってそうだな。
0008デフォルトの名無しさん
NGNG100万って、どのくらいなんだ？
100万文字ってなんKB？
0009デフォルトの名無しさん
NGNG0010デフォルトの名無しさん
NGNGBCDで持つとして、100万桁=50万バイト=5000キロバイト=5メガバイト
0011特殊早大生理論（広末涼子）
NGNG0013頭プアー
NGNG0014デフォルトの名無しさん
NGNG整数環の素元のこと
0015デフォルトの名無しさん
NGNG扱えるユニットがあるので
それで扱える桁数かな、とオモタわけ。
5MB程度ならいけるな、100万桁．．．．
ところで素数って正の数かしら？
0016 刺身定食
NGNG0017デフォルトの名無しさん
NGNGそのリストと照らし合わせる判定方法にすれば
それこそ速い遅いの問題ではないっていってるんだが？
通じないか？
0018デフォルトの名無しさん
NGNGリストってお前、500bitとか1000bitとかの素数まで全部保持するつもりか？
0019デフォルトの名無しさん
NGNG>32bit整数MAXまでの素数リストがあれば実用上何も問題なし。
って言ってんだろ?
0020デフォルトのやる気なし
NGNGHNにちょっとだけワラタ。 (^◇^;
0022デフォルトの名無しさん
NGNGメモリ管理とか考える方が先。
エラストテネスで32bitMAXまでの素数リストを作ろうと思っても
メモリがオーバー風呂して出来なかったYO。
あと、何万桁を計算する手法も考えるべき。
32bitとかそういうレベルじゃなく
それでいて素数判定が早いというのには意味があるが
実用的かどうかというと全然実用的じゃないと思う。
学術的には意味はありそうだな。
0023デフォルトの名無しさん
NGNGRabinの定理というものが知りたかったのか...
で、その定理の詳細ページどこかにない？
0024デフォルトの名無しさん
NGNGこれって適当に答え提出して「これは素数だ」って言い張ったらどうなるのだろう？
0025デフォルトの名無しさん
NGNGhttp://www.em.cache.waseda.ac.jp/research/complexity.html
確率性アルゴリズム
補題３
n≧2 が素数である必要十分条件は 0＜a＜n なる任意の a に対して
a^(n-1) ≡ 1 (mod n) が成り立つことである．
補題４
n が奇数の合成数ならば，
次の条件を満たす整数 0＜a＜n が少なくとも (n-1) / 2 個存在する
a^(n-1) ≡ 1 (mod n) が成り立たないか，または，
2i | (n - 1), 1＜gcd(a^((n-1)/2i -1), n)＜n を満たす整数 i が存在する．
これを使って 適当に乱数で検査して 素数らしさが十分だったら、それで
いいってのが１が欲しい乱数の作成方法でしょ？
0026デフォルトの名無しさん
NGNG素数が注目されているのは暗号化技術に必要だから、
例えば、公開鍵暗号には２つの素数 p,qが必要。
この素数p,qと別の数 e,dの関係が
e*d=n; n % (p-1)*(q-1) == 1 になるように e,dを見つけて
e と nを公開鍵 dを秘密鍵とする
暗号化は (平文^ e) % n として行い 解読は (暗号文^d) % nとする
この一連の流れが巧くゆくのはそれは p,qが素数だったからだけど
実際に必要なのは e,d,n
逆にこの一連の流れさえ巧くゆくなら p,qが素数であったかかどうか
を問い詰める必要はない。 逆に考えれば、素数らしいｐ、ｑさえあ
れば、試してみて復元出来るならそれでいい。
という事で、p,qがホントに素数かどうかより 素数らしい値を手軽に
得られる方法でいいという理屈。
002726
NGNG○ n=p*q; (e*d) % (p-1)*(q-1) == 1
0028デフォルトの名無しさん
NGNGMax32ビットの素数リストすら作れない、
メモリの確保が不可能だ。
37万くらいまでなら泥臭い方法で出来たけど。
10桁ですら、とてもとても。
0029デフォルトの名無しさん
NGNG0030デフォルトの名無しさん
NGNGこのあたりまでのリストなら
ならなんとか....
0031デフォルトの名無しさん
NGNGエラストネスのふるいを使えばどうでしょうか？
基本的にある数Nが素数かどうかはN^1/2までの約数が無いかを調べればよいので、
N^1/2を超える最小の1からのすべての素数の積(変な日本語ですいません)
までのパターンがあれば、後はNをその積で割った余りを検証すればよいのでは？
1*2*3*5*7*11*13*・・・・・・って計算していって0xFFFFを超えた数を見つければ
良いと思います。
003231
NGNGエラストテネス？恥ずかしい・・・
003331
NGNG0034デフォルトの名無しさん
NGNG中途半端な知識だな。
それはRSAな。
○RSA∈公開鍵暗号
×RSA＝公開鍵暗号
×RSA∋公開鍵暗号
003554
NGNG何の略か忘れたけど・・・
003634
NGNGRSA暗号は発案者三人の頭文字をとっている。
Rivest
Shamir
Adleman
で、こいつらが作った会社の名前もRSAセキュリティ
0037デフォルトの名無しさん
NGNG>>31
つかってるよ。
ただ、エラスネスするためには
配列やリストを用意しなければいけないんだけど
Booleanの配列で
32bit個の入れ物をつくっちまうと
メモリオバなのよ。
直接メモリのビットを見ればいいんだけどねえ。
軟弱プログラマだからやり方しらないよ。
0038荒井注
NGNGエラトステネスだバカヤロー
0039デフォルトの名無しさん
NGNGかかった（PenMMX266）ので途中で断念した。時間かけ
ればできるけど、それでもまだ１０桁だよ。１００万桁なんて
想像できん。
0040デフォルトの名無しさん
NGNG0041デフォルトの名無しさん
NGNG0042デフォルトの名無しさん
NGNGどういうアルゴリズムを使うわけ？
素数の定義からは逃げられないよ。
0043デフォルトの名無しさん
NGNG0044デフォルトの名無しさん
NGNG激しく同意。
0045デフォルトの名無しさん
NGNGめちゃくちゃややこしいし
一つの数が素数か否かを判断するのにも
膨大な時間がかかる。
この部分を分散計算するという考えそのものが思いつきもしない。
ｳﾂﾆｬ
0046デフォルトの名無しさん
NGNGBCDや１００進法で 500KByte 掛算時のΣai*biで桁溢れがおきないから大丈夫
でも、素数かどうか判定するのはどうかな？
１の言うRabinの定理で確実に素数らしい奴の予想をつけといて、それを皆で割って
みるしかないとすると、 参加者が何人いてもあんまり助けにならないよね
0047デフォルトの名無しさん
NGNG>素数とは、1、3、5、7、11、13など、
>1およびその数自体でしか割り切れない数
大学４年生にもなってはじめてしりましたよ。
0048デフォルトの名無しさん
NGNG....99999677/99999703/99999721/99999773/99999787/99999821/99999827/99999839/99999847/99999931/99999941/99999959/99999971/99999989/
0049デフォルトの名無しさん
NGNG重箱隅だが2も素数だよ。
0050デフォルトの名無しさん
NGNG0051デフォルトの名無しさん
NGNGで、メモリ直読みでエストラネスをやってみようか。
符号なし32bitMAX程度は出したいですね。
4294967295
でも
符号なし64bitMaxの18446744073709551616
は出来そうにない...
これでも100ﾏｿ桁まであと....
0052デフォルトの名無しさん
NGNGすごい桁数扱えるのね。
妙な所を
きっちりと作ってやがんな＞＞MS
0053デフォルトの名無しさん
NGNG0054デフォルトの名無しさん
NGNG0056デフォルトの名無しさん
NGNG>>53-54どこみていっているんだ？
005756
NGNG0058デフォルトの名無しさん
NGNG10^(10^6)〜(10^(10^6) + 2^24) の範囲には素数がありそうだと思える
これらの全部のフラグをたてといて、
Rabinのアルゴリズムで最後迄残った数個に絞り込むのはそれほど大変じゃない筈
後は、これが素数かどうかを証明する方法だなあ・・・
0059デフォルトの名無しさん
NGNG１と自分以外で割り切れるかどうか、試せばいいんじゃないの？
0060デフォルトの名無しさん
NGNGいや、あの、それじゃ１命令10GHzで多桁除算が出来たたとしてもさ
10^(10^6)/10^10 = 10^(10^6-10)秒かかるんだよ
10^(10^6) のオーダだとフルイに使うメモリにしたって
世界中のパソコンのハードディスク足したって足りはしないさ。
006158=60
NGNGそれから √n迄でいいから、上の計算よりずっと少なくて
√(10^(10^6) )/10^10 = 10^(1000-10)=10^990秒でいい
0062デフォルトの名無しさん
NGNGFFT方式だとどれくらい早くなるのかな？
0063デフォルトの名無しさん
NGNGハゲしく同意
フルイにかけるメモリ確保すら問題
0064デフォルトの名無しさん
NGNGだれか量子コンピュータを作ってよ。
0065デフォルトの名無しさん
NGNGfor (c=2;c<√(a+b);c++)
{
if(cが素数でないなら)continue;
if( (a+b) % c !=0) goto 素数でない
}
だったらいいんでしょ？
(a+b)%c = (a%c + b%c)%c だよね？
a=1000000! としとけば、100万迄の a%cはゼロだから
計算量は100万分の１になるんじゃないかな？
0066デフォルトの名無しさん
NGNG0067デフォルトの名無しさん
NGNGn! = n^n * e^(-n) * SQRT(2*pi*n) e^(1/(12*n) - 1 / (360 * n^3) + θ/(1260 * n^5))
(0<θ<1)
です。
006865
NGNG100万分の１じゃなくて100万分の節約ってだけだな
0069非65
NGNGLog[2.0,Gamma[1000000.0]] ≒ 1.85 10^7
1000000!は10進数で約557万桁、約1850万ビットの数ですか....
想像つきませんね。
0070デフォルトの名無しさん
NGNG100万桁よりも長い最初の素数 5万ドル
10億桁よりも長い最初の素数に対する最高金額の25万ドル
だからなあ
0071デフォルトの名無しさん
NGNGhttp://www.eff.org/awards/coop.html
狙うならもう１桁上だね
0072デフォルトの名無しさん
NGNG求めるならメルセンヌ素数を求めたいね。
0074デフォルトの名無しさん
NGNGなんで桁数が10進法で上昇しなければならないんだ？
意味わからん。
0076 １
NGNG僕が参考にしていた文献です。
これをもとにして、高速な素数判定ソフトを作りたかったんですが、
無学なもので、内容が理解できませんでした。
そういう悲しい事情です。
【参考文献】
「アルゴリズムの基礎」
五十嵐善英、西谷泰昭 共著
コロナ社
ISBN4-339-02348-5 \3,000
≪8.5 Rabinの素数判定アルゴリズム（P168−170）≫
0077デフォルトの名無しさん
NGNG32bitのintの最大値
2147483647
は素数なんだけど、これを泥臭く2.3.4....2147483646でわってみる判定方法で
やったとしても今のマシンなら数分(今試したよ)なのだから
高速素数判定ソフトを作っても
何の意味もないような気がする。
32bitや64bitを超えた数値を扱う方法は
アルゴリズム以外の面で実装が大変だ
>>1さんは
その辺りの事、わかってらっしゃいますか？
2147483647までの素数列を全て求めるときの
高速アルゴがあるのなら教えてもらいたい。
アルゴリズムというより
そのデータを保持するメモリ管理の方が知りたいな。
0078デフォルトの名無しさん
NGNGソースはこれから辿って下さい
http://search.net-newbie.com/php/function.gmp-prob-prime.html
0079１
NGNG実装が大変なんですか。あまりよくわかっていませんが、
桁が多いと、扱える変数（？）がなくなると言うことですか？
それにともなって、細工が必要になるのかな。素人でスマン。
＞2147483647までの素数列を全て求めるときの
＞高速アルゴがあるのなら教えてもらいたい。
僕が最初、試したアルゴリズムは、
６ｎ±１（ｎ≧１、ｎは自然数）を、素数の候補として、
手当たり次第に√（６ｎ±１）に達するまで、
（６ｎ±１）ｍｏｄ（既に得られてる素数）＝０
を判定する方法をとったんですが。これってダメ？
でも、これだと、2147483647までできないと思う。
やってみてもいないし。メモリ管理に至っては全然です…
0080デフォルトの名無しさん
NGNG0081デフォルトの名無しさん
NGNG○ 2^31
それから2^31迄の全ての素数を入れるのに必要なメモリは
およそ 4*( 2^31/ ln(2^31)) だいたい 400Mバイトだから フラグで記録した方が効率がいい
0082１
NGNG>>78
＞やったとしても今のマシンなら数分(今試したよ)なのだから
＞高速素数判定ソフトを作っても
＞何の意味もないような気がする。
＞その高速判定は、あくまでも確率的な判定方法ですよ
文献によると、
≪Rabinらはこの方法により、２＾ｐ−１なる形の整数が素数であるかどうかの判定を計算機を用いて行った。1<p<500の全ての整数についての判定をたった数分で行い、全ての判定は正しかったと報じている。この方法は素数であるかどうかを判定するどの決定性アルゴリズムより桁違いに速く、サンプル数ｍを適当な大きさにすることにより、誤りの確率をきわめて低くすることができる。≫
≪例えば、ｍ＝５０のとき、誤り確率は約千兆分の１となり、基本演算の所要時間が１μsで、３万時間連続稼働している間に１度の割にしか基本演算の誤動作が生じない計算機よりも信頼性が高い判定方法である。≫
とあります。
これ見て、作ってみたいなと思ったわけです。
0083１
NGNGアドバイスありがとう。そういう計算方法も知らなかったんです。理系だけど。ちょっとやってみます。
0084デフォルトの名無しさん
NGNGでいいからね
0085デフォルトの名無しさん
NGNG＞ n が奇数の合成数ならば，
＞ 次の条件を満たす整数 0＜a＜n が少なくとも (n-1) / 2 個存在する
＞ a^(n-1) ≡ 1 (mod n) が成り立たないか，または，
＞ 2i | (n - 1), 1＜gcd(a^((n-1)/2i -1), n)＜n を満たす整数 i が存在する．
X ≡Y mod n は 両辺の nの剰余が等しいという意味。
もしかして gcdという関数が判らないのか？
0086デフォルトの名無しさん
NGNGU32* GetSosu(void)
{
U32 *p;
int i, j;
p = malloc(sizeof(*p) * 2);
p[0] = 1;
p[1] = 2;
for(i = 3; i != 1; i += 2){
/* 素数判定 */
for(j = 1; (j <= p[0]) && ((i >> 1) > p[j]); j++){
if(0 != i % p[j]) break; /* 素数発見 */
}
if(j > p[0]){
/* 新しい素数を追加 */
p[0]++;
p = realloc(sizeof(*p) * (p[0] + 1));
p[p[0]] = i;
}
}
return(p);
}
仕事が忙しいから作ってみたよ。
テストはしてないんで、
コンパイルが通るかどうかは知らない。
何とか方とか難しいことは知らない。
っていうか、実は素数の定義自体も自信ない。
速度も別に気にしてない(毎回realloc()とかしてるし)。
あと、メモリ確保できなかった時の対処もしてない。
とりあえず、答えが出るのが何時になるかは知らないが
出せないなんて事は無と思う。
0087デフォルトの名無しさん
NGNGhttp://www.ml-labo.com/download.htm
これはどうだ
アルゴリズム(ラビン法): 予め決められた整数L と判定対象の整数p に対して以下を実行.
p = 素数の候補;
p = 1 + 2 k w を満たす k と w を計算（ただし, w は奇数）;
以下を L 回繰り返す {
1 < a < p なる a をランダムに生成;
y = a w mod p;
j = 0;
while ( ( j < k ) and ( y ≠ 1 ) and ( y ≠ p -1 ) ) {
y = y 2 mod p ;
j = j+1;
};
if ( ( j = k ) or (( j > 0 ) and ( y = 1 )) ) p は素数でない;
};
途中の whileが gcd だ。 意味が判らなければ gcd とユークリッド で検索してみろ
0088１
NGNGそうなんです。わかりません。gcdという関数の存在も今知りましたが、
それでも、文章の意味が分からないんです。
しかも、どこにツッコミ入れられているかもわからないんです。
つきはなさないで下さいね。
0089１
NGNGお忙しい中、ありがとうございます。
力不足ですいません。
時間かかりますけど、やってみます。
0090デフォルトの名無しさん
NGNGShorのアルゴリズムも、Rabinと同程度に確率的な判定法でしかないっす。
0091デフォルトの名無しさん
NGNG0092デフォルトの名無しさん
NGNGShorのアルゴリズムとは
0093デフォルトの名無しさん
NGNG和田秀夫 / 遊星社
最新の資料ではないがよくまとまっているし、プログラムコードも入っている。
専門家でなければ必読。
0094デフォルトの名無しさん
NGNG日本の東北大の教授とかが書いたやつなんてコードに無駄なとこ多くて読む気がうせる
0095デフォルトの名無しさん
NGNGこれだね
http://www.geocities.co.jp/SiliconValley-Oakland/8522/data/mformula.txt
モジュラベキとかgcdとか関数を増やしてけば使えるかも・・・でも
既に1000桁扱える１０進basicとかあるからなあ
アドバンテージはメモリがある限り桁数増やせる点だけど、今のアルゴリズム
じゃ掛算が遅くてダメだな
0096デフォルトの名無しさん
NGNGおおざっぱに書くと
１．Ｎを素因数分解したいとする。
----- 因数が見つかるまで２〜５を繰り返す ------
２．Nと互いに素で x<N な x を取ってくる
３．x^r == 1 (mod N) になる最小の自然数 r を計算
４．rが偶数になる確率は1/2程度。
５．rが偶数の時 gcd( N, x^(r/2)±1 ) のどっちかが、1やNでないNの約数を与える確率も1/2程度。
-------------------------------------------
ステップ３．を計算するのは量子コンピュータじゃないと死ねるので結局まだ実現不可だけど。
0097デフォルトの名無しさん
NGNG掛算で 10進 40万桁超で正しい答えが得られない筈
掛算途中の桁あげを処理しないといけない
でもこの修正をやっても実際は遅くて使い物にならないから意味ないけどね
function s100Mul(const a,b:s100):s100;
var i,ia,iaL,iaH:integer;
var cy,ecy:integer;
begin
SetLength(Result,Length(a)+Length(b));
cy:=0;
for i:=0 to High(Result) do
begin
iaL:=Max( 0,i-High(b));
iaH:=Min(High(a),i );
ecy:=0;
for ia:= iaL to iaH do
begin
cy:=cy+a[ia]*b[i-ia];
if cy>100000 then
begin
cy :=cy-100000;
ecy:=ecy +1000;
end;
end;
Result[i]:=(100+(cy mod 100)) mod 100;
cy:=ecy+(cy-Result[i]) div 100;
end;
s100CyNorm(Result,cy);
end;
0098デフォルトの名無しさん
NGNG2.3.5.7.9.
10を調べると2で割れるから素数ではなく
11を調べるとこれまでリストした全ての素数のうち
11/2=5以下のもので割り切れないから素数
2.3.5.7.9.11.
12を調べると2で割れるから素数ではなく
13を調べるとすべての素数で割り切れないから素数
2.3.5.7.9.11.13.
14を調べると2で割り切れるから
15を...3で割り切れる
16を...2で割り切れる
17を...これまでリストした全ての素数のうち
17/2=8以下のもので割り切れないから素数
2.3.5.7.9.11.13.17
とやっていくと
計算量が減り、まあまあ高速に素数リストが出来て
メモリのことも考慮しなくていいようになるかも。
ダメ？
0099デフォルトの名無しさん
NGNGここは！！！
http://www.ie.u-ryukyu.ac.jp/ie/mirror/www.iip.ie.u-ryukyu.ac.jp/yuu/c_html/sosuu.html
010099
NGNG0101デフォルトの名無しさん
NGNGそれより計算量が多かったり低速だったりする方法って思いつかないぞ。(w
√N を計算するより N/2 以下の素数全部で割ってみる方がマシってこと？
0102デフォルトの名無しさん
NGNGいや、よくわかってないもので。
>√N を計算するより N/2 以下の素数全部で割ってみる方がマシってこと？
詳しく教えて。
ルートエヌ？？
エラトステネスより98の方が現実的だと思ったのさ。
0103デフォルトの名無しさん
NGNGJavaなのに速いな....
0104デフォルトの名無しさん
NGNG＞基本的にある数Nが素数かどうかはN^1/2までの約数が無いかを調べればよいので
欝...
0105デフォルトの名無しさん
NGNGこれまでに確認した全ての素数で割る方法
N以下の素数の数はおよそ N/ln(N)
M迄求める計算量は ΣN/lnN 回の割り算
ただし、素数でなければ途中で割り切れるので、殆ど負荷がないとすれば
ΣN/ln(N) /ln(M) 回の割り算
いっぽう篩方式は およそ M/ln(M)*Σ(M/N))回の掛算
殆ど同じオーダだけど、片方は割り算、片方は掛算だから圧倒的
フルイを使う方が早い
0106デフォルトの名無しさん
NGNG判定したい値をNとしたら
N=a*b+c ただし c<a とする
つまり N/a=b+c
と表現した時の cが0になれば素数でないという事
aが小さい方から順に判定をした場合、実は bをN-1から小さくして
ら順にテストしてるとも言えるわけで、
a>bになったら既にテストした事を繰り返してる事になってしまう
その交差点はa=bの時の a^2=N つまりa=√N の点
0107デフォルトの名無しさん
NGNGスマSO、なぜかうちでは
フルイの方が圧倒的に遅い
これまでに確認した全ての素数で割る方法 でおよそ2分
フルイで10000000までの結果(>>30)を出したときがおよそ25分
どういうことじゃ？
そのうちソースだすよ（Delなんだけど
どこかにくだらないオーバーヘッドがあるのかもしれない。
出力時に全ての数値をなめなければ結果を出力できんし。
(わかりやすいように計算と結果は別々にしているんだ。)
メモリ効率の点からも
全ての素数で割る方法じゃないと
64(or63)bitMAXまでは求めるのが無理らしいです。
2^63ビット確保できないってさ。
>>106
ありがとう。めちゃためになる。
数学科出身？？
0108デフォルトの名無しさん
NGNGふるいに使うのは足し算だけで十分だと思われ。
今とりあえず組んでみたら1億までで30秒くらい。
0109デフォルトの名無しさん
NGNGhttp://arch.arch.kumamoto-u.ac.jp/hagane/yamanari/joho2/kougi3.html
ここのフルイだと
For I = 2 To Limit
If sosuu(I) = 1 Then
For J = 2 * I To Num
If J Mod I = 0 Then
sosuu(J) = 0
End If
Next J
End If
Next I
こんなんで、J Mod Iってのが入っているんだけど
このやり方がダメなの？
0110デフォルトの名無しさん
NGNGダメダメダメダメダメ。でもうちの新人に[i...j]のうちのnの倍数を
全て出力するプログラム、って課題を出したら、最初は確かに mod (%)
使ってた...
0111ャfフォルトの名末ｳしさん
NGNGある数がある数の倍数かどうかを判断するとこがダメなのか？
0112デフォルトの名無しさん
NGNGたとえば１０進数なら
３の倍数かどうかは全部の桁を足して３の倍数か見るだけでいい
５の倍数かどうかは最下位が０か５か
っていうのがあるよね
２進数の場合はどうだろう？
３の倍数 ２桁毎に足した結果の下２桁が11である
５の倍数 ４桁毎に足した結果の下４桁が1001
という事で、Nで割れるかどうかは
(2^n-1)/N で割り切れるnを探して n桁毎に足した結果がNなら
というのでいいのかな？
0113デフォルトの名無しさん
NGNGでも、>>112の疑問もお願い
0114１１２
NGNGn桁毎に足して、その結果がさらにNの倍数かどうか見なくちゃいけない
３の場合 3*17=$33 ２進で２桁毎に足したら６だから 0110ｂ さらに２進２桁毎に足したら 11ｂ
0115デフォルトの名無しさん
NGNG0116デフォルトの名無しさん
NGNGほんとにバカなんで教えてくださいませ
倍数判断はともかくとして
>>109って
For I = 2 To Limit
If sosuu(I) = 1 Then
For J = I+1 To Num
If sosuu(J) = 1 Then
If J Mod I = 0 Then
sosuu(J) = 0
End If
End If
Next J
End If
Next I
処理系依存だとは思いますが
こう書いた方が早いです…
Modの部分の倍数か否か比較の高速化以外には
高速化出来ないでしょうか。
For J = I+1 To Num
この部分が怪しいような気がする。
0117デフォルトの名無しさん
NGNG各桁ごとを+-反転させながら足した総和が3の倍数ってことになるね。
例 18 = 10010(2)
-1+0-0-1+0=0が3で割り切れるので3の倍数。
5の倍数だと、前述のようなやり方で、今度は2桁ずつセットにして考えればいい。
例 18 = 10010(2)
+1-00+10 = 11(2) = 3 だから、5で割ると3あまる。
多分実用性はないな。
0118117
NGNG訂正 +1-0+0-1+0=0 だ。
0119１１２
NGNG3 2
5 4
7 3
11 10
13 12
17 8
19 18
23 11
29 不明
31 5
37 36 <- ３２より長いからこれも使い難いな
41 20
43 14
47 23
0120１１２
NGNG３，５、１7 の倍数かどうかチェックすれば一挙に３つ計算出来る
から、多少は早くなるんじゃないかな
0121117
NGNGpが素数なら、2^(p-1)-1 が p で割り切れるよ。
0122デフォルトの名無しさん
NGNGどのくらいの桁数なの？
0123１１２
NGNGあ、ホントだ という事は 29は 28bitで循環か、なんで抜けたんだろ 浮動小数点使ったせいか
>>122 必要なだけ 100桁くらいあれば個人ならなんとか
千桁くらいあればとりあえず
0124デフォルトの名無しさん
NGNGIf sosuu(I) = 1 Then
J = I * 2
Do While J <= Num
sosuu(J) = 0
J = J + I
Loop
End If
Next I
なんかの倍数なんてチェックする必要すらないような気がするが...
0125１１２
NGNG2 最下位桁のみ
3 1桁単位に加算
5 最下位桁のみ
7 6桁単位に加算
11 2桁単位に加算
13 6桁単位に加算
37 3桁単位に加算
41 5桁単位に加算
73 8桁単位に加算
101 4桁単位に加算
137 8桁単位に加算
0126１１２
NGNGJ=I*2 にしてるけど、Jより小さい素数の倍数は不要でしょう？
説明すると、
まずは２刻みに４以上を全部消すよね？
次は３刻みなんだけど６は既に消えてるし２の倍数は消えてるから９から消す
その次は５刻みだけど
やっぱり５＊２ も５＊３も消えてるから５＊５から１０刻みに消せばいい
つまり素数N番目は N*Nから 2N刻みに消せばいい
0127１１２
NGNGで、FOR文がTO 以下なら繰り返すんだとして
Limit = Int(Sqr(Num))
For J = 4 To Num Step 2:sosuu(J) = 0:NextJ
For I = 3 To Limit Step 2
If sosuu(I) = 1 Then
For J = I*I To Num Step 2*I
sosuu(J) = 0
Next J
End If
Next I
FORについてる I*I や 2*I がループ中再評価されるんなら外に出してね
0128１１２
NGNGMAXNに設定出来る最大は FFFFfffD だから注意して
#include <stdio.h>
#define MAXN 0x7FFFfffFU
static unsigned char *tbl;
void ClrTable(unsigned int n){
unsigned int a=n>>4U;
unsigned int bpos=1<<( (n>>1u) & 7u);
tbl[a]|= bpos; //0なら素数 1ならそうでない
}
bool TestTable(unsigned int n){
if (n&1){
unsigned int a=n>>4U;
unsigned int bpos=1<<( (n>>1u) & 7u);
return 0== (tbl[a] & bpos);
}else
return (n==2);
}
main()
{unsigned n;
int SqrMAXN=0xFFFF;
for(int i=0;i<10;i++){int r=MAXN/SqrMAXN; SqrMAXN +=(r-SqrMAXN)/2;};
SqrMAXN++;
tbl=new unsigned char [1+MAXN/(8u*2u)];
if (!tbl){ printf("ERROR\n");return 1;};
for(int i=0;i<MAXN/(8u*2u);i++)tbl[i]=0;
for(n=3;n<=SqrMAXN;n+=2)
if(TestTable(n) )
{ unsigned nn=n*n , n2=2*n;
printf("%5u\n",n);
for(unsigned i=nn;i<=MAXN;i+=n2)ClrTable(i);
}
for(;n<=MAXN;n+=2)
if(TestTable(n) )printf("%5d\n",n);
return 0;
}
0129１１２
NGNG2147483477
2147483489
2147483497
2147483543
2147483549
2147483563
2147483579
2147483587
2147483629
2147483647
途中まで遅いけど、１０００を越えた付近からスクロール速度の方が先に限界になります。
ハードディスクにリダイレクトさせると１ギガを超えそうになったので、プログラムを弄って後ろの方だけ印刷させました
0130１１２
NGNG画面に出したらスクロールに３分かかった
0131デフォルトの名無しさん
NGNGそのスレをたてた「悲惨な1」です(藁)。
その後、以下の文献をみつけたので問題は解決しました。(17ページにアルゴリズムの工夫の仕方が載ってます)。
ftp://ftp.rsasecurity.com/pub/rsalabs/rsa_algorithm/rsa-oaep_spec.pdf
ところで、Webで以下のような記述をみつけました(これは素数判定ではなくて素数生成の話題)。
ttp://www.mcc.pref.miyagi.jp/people/ikuro/koramu/prime.htm
＞ １９７１年、旧ソ連のマチアセビッチは、素数全体をあるひとつのディオファントス
＞ 方程式の解として表すことにも成功しています。すなわち、すべての素数をつくる式
＞ を生み出したことになるのですが、その式は３７次で２４個の変数をもつ多項式と
＞ ２１次で２１個の変数をもつ多項式でした。
興味があり調べているのですが、なかなかこれに関連する文献がでてきません。なにかご存知の方はいませんか。
0132デフォルトの名無しさん
NGNG???
0133１１２
NGNG＞それは１９変数の多項式であり、１９個の自然数をその
多項式ｆに代入してｆの値が正であればそれが素数であり、
逆にすべての素数についてもｐ＝ｆ（x_1, x_2, ..., x_19）を満たす自然数 x_1, x_2, ..., x_19が存在する。
という事だから、あれば便利そうだけど、２０次とかの式を含んでるのなら例えば500bit前後に収まる
ものをというような要求には難しいんじゃないでしょうか？
オイラーの式 n2+n+41 の方がコントロールし易いんじゃないかと思います
0134デフォルトの名無しさん
NGNG0135デフォルトの名無しさん
NGNG１００万桁の最初の素数を見つけるような場合には意味が無いとおもわれ
0136デフォルトの名無しさん
NGNG（N−１）！＋１ ≡ 0 mod N
if ((（N−１）！＋１ )% N==0) 素数です else 素数ではありません
0137デフォルトの名無しさん
NGNG存在する訳じゃないし、それを満たせば全部素数って事もないよねぇ。
ていうか当然か。。。
0138デフォルトの名無しさん
NGNG代入しても負の値が出ることがかなり多いそうだよ。(つまり素数じゃないやつネ)
0139デフォルトの名無しさん
NGNG話それるけど、
そのようなプログラムでコンソールに連続的に表示するとき、
Ctrlを押してるとスクロール（表示速度）が速くなるんだけど
理由知っている人いる？
0140デフォルトの名無しさん
NGNGfor 文で MAXN 迄検索してる個所で MAXNが 2^32 に非常に接近するとヤバイよ
初期値より小さくなったらという条件も入れた方がいい というか MAXN=2^32-1 として
× for(unsigned i=nn;i<=MAXN;i+=n2)ClrTable(i);
○ for(unsigned i=nn;i>nn;i+=n2)ClrTable(i);
× for(;n<=MAXN;n+=2)
○ for(;n>2;n+=2)
とオーバフローしたら検出して抜ければいい
0141デフォルトの名無しさん
NGNG＞1000000! % 1000001 はもう予測出来ないな
１３６が書いているウイルソンの定理から 1000001 が素数なら
1000000! % 1000001 =1000000 です
しかし、1000001=101*9901なので階乗は割切れますから結果はゼロです
1000003が素数なので
1000002! % 1000003 =1000002
1000001! % 1000003 =1
1000000! % 1000003 =500001
と予測出来ます
0142デフォルトの名無しさん
NGNGついていくのに必至だ。
>>48で7時間って言ってたけど
それが5分になった（ｗ
うれしい！
手元のマシンにメモリが少ないのが悲しい所。
0143デフォルトの名無しさん
NGNGhttp://www.mogya.com/develop/kanjicalc/develop.html
10億桁扱えるユニットがあるので
メモリを食わない＆計算時間が1日程度の方法で求められるのなら
適当なルーチン(わかりやすくたのんます)書いてくれたら試すよ。
0144デフォルトの名無しさん
NGNG1Gマシンで100万桁の掛算だけで既に１日で出来ないと思うよ
0145デフォルトの名無しさん
NGNGそっか。やはりWiredNEWSで言われているように
分散処理で考えないとだめか。
100万桁といわず、ここでどんどん桁数を上げていきません？
0146１４４
NGNGそういうふうに動的に取るとメモリが細切れになって大きな配列が最後は取れなくなる
C=A*Bという計算をする時に Setlength(C,length(A)+Length(B))とメモリを確保してから
余ったら開放する方式の方がいい
で、かける方の結果の桁でループを回すんじゃなくて、>>97のように
かけた結果の桁でループを回すといいよ。その方が1桁以上早くなる筈
0147１４４
NGNG使い方は文字列で与えて文字列で出て来るから殆ど同じだし
入力で "2^1000+1" みたいな書き方も出来るし
内部１００進だから１０進に比べてそれだけで掛算が４倍早くなるし
どっちにしても、掛算について万桁超えたらfft使わないと苦しいと思う
解説はここ
ttp://homepage2.nifty.com/m_kamada/docs/fftmul.htm
0148デフォルトの名無しさん
NGNGunit mformula;
も数値を文字列で持っているんだね。
でも、なんとなく使い方が難しそうなのだが...
Stringで数値を保持しておいて
計算するときはFormulaNumStrして
結果もStringに保持しておけばいいって事なのかな？
0149１４４
NGNG例えば1万桁台の最初の素数を探してみたいとか
ルカステストで100万桁台のメルセンヌ素数を調べたいとか
やっぱり特定の目標にしないと、 桁溢れとか色んな問題がクリア出来ないと思う
0150１４４
NGNGfunction N100Str(str:String):String;
begin
Result:=FormulaNumStr(TN100Num,str');
end;
を自前で作って Caption:=N100Str('100^10-1');
みたいに試したからそんなに難しいとは思わなかったけどね
無理に関数使おうとすると、クラス貰って来てどうとかやっかいだね
どうせ長桁になったらインタプリタの時間なんてしれてるから インタプリタに任せた方がいいんじゃない？
小数点の使える TN100の他に 整数だけのs100 もあるみたい
あんまり試してないけど、素数見つけるならこっちでいいんじゃない？
0151デフォルトの名無しさん
NGNG今日こそは金曜日なので
月曜日まで計算させて
32bitMax値をたたき出そうと
思っていたんだが。
63bitMaxまでいっけるっかな。
最初に、"素数は決めうちでリストでもっておいて"
それにヒットするか否かで判断しろっていったけど
これだけ莫大な量だと
その方法が正しいのかすら怪しいな。
0152１４４
NGNGa,b,cを32bitとして cは素数として
(a*2^32 + b) % c が0かどうかを検査すればいいんだけど、
64bitの数としてではなく a,bのペアのままで素数を求めたい
(a*2^32 + b) %c
= ((a*2^32)%c + b) % c
= ((a%c)*(2^32%c))%c + b%c) % c
これであってるでしょうか？
後は 2^32 %c をどう効率よく計算するかと
加算時に桁溢れがあった場合の処理
s^32 %c は >>119 で調べたnを使って n<32なら
(2^32) % c = (2^(32%n)%c
nは普通 n=c-1 だけど、それより小さい場合もあるみたい
後、
(2^32) % c =(2^16)*( (2^16) % c) % と分解出来るかな？
0153デフォルトの名無しさん
NGNG色々コメントありがとうございます。
掛算のFFT化は時間が出来れば挑戦してみます
でも、素数判定に必要なのはどっちかというと
高速なモジュラ除算なんじゃないでしょうか？
>>152
64bit迄なら 64bit/32bitの割算の余りさえ求まればいいので
直接計算する命令がCPUにあります。
Delphiの場合 符号無しの64bit整数が無いのですが
無理やりint64を使って
function U64Div32Mod(x:int64; y:longword):longword;pascal;
var xx:longword absolute x ;
asm
MOV EAX,xx
MOV EDX,xx+4
DIV y
MOV EAX,EDX
end;
こうアセンブラで書く事が出来ます。
0154学生
NGNG「掛算のFFT化」とは？なんですか？
FFT？
0155デフォルトの名無しさん
NGNG私はとりあえず、泥臭い方法(>>98)で
普通の整数を使えば32bitまでは
出せそうなのでそれでやってみます。
土日計算させてたら、答え出るんじゃないかと期待
64bitは計算時間的に不可能だろうな。
0156デフォルトの名無しさん
NGNGFFT=高速フーリエ変換
0157学生
NGNGありがとう。
じゃ、
「掛算のFFT化」とどういう意味 ？
0158デフォルトの名無しさん
NGNG少しはさかのぼって読んでみたらどうか。>>147を読め。
0159153
NGNG理由はオーバフロー例外が起きてしまうからです。
Delphiのint64の割算はCPUの命令を使ってないから変だなと思って
たら、こういう落とし穴があったんですね・・・・
>>155
私は こんなテストルーチンを書きました。
// base からcnt範囲の素数の表を作る
type primeFlag=array of ByteBool;
function PrimeList(base:int64;cnt:integer):primeFlag;
var i,n,n2:longword ;
var ends:int64;
begin
ends:=base+cnt;
SetLength(result,cnt);//素数テーブル
for i:=0 to cnt-1 do Result[i]:=true;
n:=2; //2の倍数を全部消す
i:=(n-U64Div32Mod(base,n)) mod n;
while i+base <= n do inc(i,n);
while i<cnt do begin Result[i]:=false; inc(i,n); end;
for n2:=1 to $7FFFFFFF do begin //16bit迄の全部の倍数を消す
n:=n2*2+1;
if n>= ends then break;
i:=n-U64Div32Mod(base,n);
while i+base <= n do inc(i,n);
while i<cnt do begin Result[i]:=false; inc(i,n); end;
end;
end;
0160153
NGNGfunction U64Div32Mod(x:int64;y:longword):longword;
begin
Result:= x mod y;
end;
とやればなんとか動きます・・・・滅茶苦茶遅いです・・・いつ終わるのか
0161153
NGNGif int64(n)*int64(n)>ends then break; にしたら
1234567890123481 1234567890123493 1234567890123527 1234567890123583
1234567890123647 1234567890123661 1234567890123671 1234567890123697
これくらいの桁なら少し待ってれば計算してくれました。けっこう楽しいです
でも、あと２桁が限界かも
0162デフォルトの名無しさん
NGNG分散システムのモデルを考えてみるのもよいかもしれない
今日このごろ。
0163デフォルトの名無しさん
NGNGそれって強力にむずかしいんだよ。
0164デフォルトの名無しさん
NGNG難しいっていうのは、要するに小分けにされた問題だけ持って帰って、
計算結果を返してくれないクライアントの事とかを考慮していくと、
強力に難しくなるって事なんでしょうか？
問題を小分けにするのが難しいってことではないよね。
0165デフォルトの名無しさん
NGNG誰かまとめてくれ。
あと％って何の演算子？
0166デフォルトの名無しさん
NGNGつまり剰余を求めるやつ。
VBだと、"Mod"。Rexxだと"/"。
0167デフォルトの名無しさん
NGNG2^n-1 が素数ならメルセンヌ素数
2^2-1=11b=3 2^3-1=111b=7
>>119 が書いてるように 3の倍数かどうかは２桁毎に区切って加算 7は3桁毎に加算
して判断してもいいから
2^4-1=1111b は3+3=３の倍数、 2^6-1=111111bも3と7の倍数 と即座に判るから素数ではない
つまり nがまず素数じゃないと2^n-1も素数じゃない
Lucas-Lehmer テストは
ttp://www.nara-edu.ac.jp/~asait/prime.htm
X[0]=4から始まって
X[i]=X[i-1]^2-2 という数を用意し、X[n-1] % (（2^n)-1) =0 という簡単な方法で検査可能
これを使ったフリーソフトで700人の仲間と高校生が記録を打ち立てた
ところで 2^n-1に比べてX[i] は非常に大きな数になる。そこでさらに Lucasは
途中で モジュラ計算をしてもいい事を発見した
ttp://alfin.mine.utsunomiya-u.ac.jp/~niy/algo/l/lucas.html
0168デフォルトの名無しさん
NGNGつまりモジュラーな２乗と引算だけ用意すればいい
2^5-1 = 31でテストすると b
4 -> (4*4)-2=16-2=14
14->(14*14)%31-2=10-2=8
8->(8*8)%31-2=2-2=0
めでたくゼロになって31が素数と判ります
0169デフォルトの名無しさん
NGNG(2^n-1)のモジュラーの場合は 2^n 以上の値が出たら下位に回して足せばいい
例えば33%15= 10_0001b % 1111b = 上位の10と下位の0001を足して3 となる
だから掛算途中で2^n以上の桁が出たら、全部を下位にシフトして足せばいい
0170デフォルトの名無しさん
NGNGまさか掛算？
真中の桁で分割して (a+b)^2=a^2+2ab+b^2で３人で計算して持ち寄るとか？
それとも、単に検査したい2^n-1のnを分けて皆で計算する？
運良くヒットした人が発見者？
0171デフォルトの名無しさん
NGNG0172デフォルトの名無しさん
NGNGそれだけで落ちるなあ。
0173デフォルトの名無しさん
NGNGそらあんた、2^28 byte の連続したメモリブロックが必要になったんでしょ？
主記憶足りてる？
最初に SetSize で最大値を設定したかい？
俺んとこではうまく動くよ。
0174デフォルトの名無しさん
NGNGどのくらいだと足るの？
TBits.Sizeの設定はしたよ。
でも、計算中に突然落ちるんだわ。
SetSizeってTStreamあたりと絡めなければならないの？？
0175デフォルトの名無しさん
NGNGvar dt,odt:array of Word;
var cy:int64;
var i,hDt,hPg:integer;
var mask:LongWord;
procedure mul2; //２乗する処理
var w1,w2:array of Word;
var wk:LongWord;
var i,j:integer;
begin
SetLength(w1,1+hDt );//メモリを確保して
SetLength(w2,1+1+hDt );
cy:=0;
for i:=0 to hDt do begin
for j:=0 to i do cy:=cy+dt[j]*dt[i-j];
w1[i]:=cy;//下位桁だけを代入
cy:=cy shr 16;
end;
for i:=1 to hDt do begin
for j:=i to hDt do cy:=cy+dt[j]*dt[hDt+i-j];
w2[i-1]:=cy;//
cy:=cy shr 16;
end;
w2[hDt ]:=cy; cy:=cy shr 16;
w2[hDt+1]:=cy; cy:=0;
wk:=(w2[0] shl (16- hPg)) and $ffff;
cy:=cy+ wk +w1[0];
dt[0]:=cy;//下位桁だけを代入
cy:=cy shr 16;
for i:=1 to hDt do begin //上位を下位に加算する
wk:=(((w2[i] shl 16)+w2[i-1])shr hPg) and $ffff;
cy:=cy+ wk +w1[i];
dt[i]:=cy;//下位桁だけを代入
cy:=cy shr 16;
end;
wk:=((w2[hDt+1] shl 16)+w2[hDt])shr hPg ;
cy:=wk+(( (cy shl 16) +dt[hDt]) shr hPg);
while cy >0 do begin
dt[hDt]:=dt[hDt] and mask;
for i:=0 to hDt do begin //さらに溢れた上位を加算する
cy:=cy+dt[i];
dt[i]:=cy;//下位桁だけを代入
cy:=cy shr 16;
end;
cy:=( (cy shl 16)+ dt[hDt] ) shr hPg;
end;
end;
----- つづく ----
0176デフォルトの名無しさん
NGNGvar Rep:int64;
begin
hDt:=((n-1) div 16);//最大桁位置
hPg:= ((n-1) mod 16)+1 ;//ビット位置
mask:=(1 shl hPg)-1;
SetLength(dt,1+hDt );//メモリを確保して
for i:=0 to High(dt) do dt[i]:=0; dt[0]:=4; //初期値を４に
rep:=n;
while rep>0 do begin
odt:=copy(dt,0,hdt+1);
mul2; //２乗する
fsame:= dt[0]=odt[0];
fzero:=true;
for i:=1 to High(dt) do begin
fsame:=fsame and (dt[i]=odt[i]);//全て同じかチェック
fzero:=fzero and (dt[i]=0); //全て同じかチェック
if (not fsame) and (not fzero) then break;
end;
if fsame then begin Result:=false;exit;end;//もはや見込みなし
cy:=-2;
if fzero then
begin
if dt[0] =2 then begin Result:=True;exit;end;//オメデトウ
if dt[0] <2 then begin //2以下なら負数になるから
for i:=1 to High(dt) do dt[i]:=$ffff;//全部が１は0と同じ
dt[hDt]:=mask;
dt[0]:=$ffff+dt[0]-2;cy:=0; //で引算する
end;
end;
for i:=0 to hDt do begin //2を引く処理
cy:=cy+dt[i];
dt[i]:=cy;//下位桁だけを代入
cy:=(cy-dt[i]) div $10000;
if cy=0 then break;
end;
Dec(rep);
end;
Result:=false;
end;
0177175
NGNGそりゃ考えたら n^3のオーダーで計算時間が増えるからなあ、やっっぱり
FFTで掛算をさせて せめて n^2 のオーダーに落とさないと・・・・
それでもやっぱり１００万桁なんてのは簡単じゃない
0178デフォルトの名無しさん
NGNG確保してない所に書いてたり読んでたり？
0179175
NGNG今は掛算を筆算方式でするから n^3 のオーダで計算時間が増えます
掛算をFFTにすれば n*log(n)*n というオーダになるから
10秒くらいの時間でnが３万の単位まで計算出来ると思います。
目標が１００万桁なら、１日動かせば１つの値が計算出来ると思います。
ttp://www.geocities.co.jp/Technopolis/1793/mersennu.html
ここのページを見ると
次の n は 400万台かもしれません 検査しなければいけない
素数は10万〜２０万ありますから、皆で１つの素数を分担して
検査するという方式でいいんじゃないでしょうか？
・・・・・その前に FFTで掛算するルーチンを作らないといけないけど・・・・
0180デフォルトの名無しさん
NGNG確保していないところに書き込んでも、その位置に自動的に size を
拡大してくれるから、落ちることはないんだけどなぁ。
逆に、少しずつ大きなところをアクセスするようなコードだと、
毎回巨大なメモリコピーが起きるから遅く見えるかもね。
ソースコード着いてないのか？
0181デフォルトの名無しさん
NGNGことは全くしていないのですが...
遅いとか遅くないとかではなくW2ｋがブルーバックで落ちるからな（ｗ
SosuuBitList := TBits.Create;
try
SosuuBitListing(InputValue, SosuuBitList);
出力コード
finally
SosuuBitList.Free;
end;
procedure SosuuBitListing(MaxValue: Integer; SosuuBitList: TBits);
var
SqrNum: Extended;
Limit, i, j: Integer;
begin
if not Assigned(SosuuBitList) then Exit;
if MaxValue < 0 then Exit;
if MaxValue = 0 then
begin
SosuuBitList.Size := 1;
SosuuBitList.Bits[0] := False;
Exit;
end;
SqrNum := Sqrt(MaxValue);
Limit := Trunc(Int(SqrNum));
SosuuBitList.Size := MaxValue+1;
//ここは通過する
SosuuBitList.Bits[0] := False; //0は非素数
SosuuBitList.Bits[1] := False; //1は非素数
for i := 2 to MaxValue do
begin
SosuuBitList.Bits[i] := True; //一旦Trueに初期化
end; //ここでダメみたい。
for i := 2 to Limit do
begin
if SosuuBitList.Bits[i] = True then
begin
j := i * i;
while j <= MaxValue do
begin
SosuuBitList.Bits[j] := False;
j := j + i;
end;
end;
end;
end;
0182175
NGNGn=4,000,000付近を狙うのだとすると
16bit単位にデータを持ったとして 4M/16=250K 以上のFFTが必要
（実数データだけなので半分のサイズのFFTでいいとして)
ttp://www.lohninger.com/fourier.html
２乗する為には２回必要だから262144 の2.1*2=4.2秒
これを400万回も繰り返させるには1年かかる・・・やっぱり計算そのものを分割して
する方法を考えないとダメみたいな・・・
0183175
NGNGb=(2^n)-1 で nが奇数なら
a^b=2となる aは 2^((n+1)/2) と b - (2^(n+1)/2)
(a^2)%b=((b-a)^2)%b
だから、最後の(a^2-2)=0の代わりに
msbが１なら(b-a)と置き換えて 2^((n+1)/2)と比べればいい
さらに a^2 %b= 2^((n+1)/2)+2 となるaを求めて同じようにする
２台のパソコンで片方は a'=a^2-2 を 片方はa'=√(a+2)を
それぞれ n/2回求めて結果を比較すればいい
0184175
NGNGa^2 %b=2*c という事は
a^2 =2*c+b*D となる数Dがあるという事
a^2 =2*c+(2^n-1)*D
=2*c+2^N*D-D D=2*cとすれば
a^2 =2*c*2^n 2*2^n=2^(n+1)とまとめて
a^2 =c*2^(n+1) nが奇数であれば
a =c*2^((n+1)/2)
つまり
a'=√(a+2)ｊは aが奇数の間は２乗するより高速
0185175
NGNG○a'=√(a+2)ｊは aが偶数の間は２乗するより高速
でも何か変だなあ２乗の方も高速計算方法があるのか >> 184が間違っているのか？
0186175
NGNG○ a = ±√c *2^((n+1)/2)
0187デフォルトの名無しさん
NGNGこの数字は無理筋？
FFTって複素数の掛算がだいたい 2*18*262144なんでしょ？
1回の掛算は４回掛算と２回の足し算だから
300MFLOPSあれば狙える数字じゃないかと思うんだけど？
というかプレステ２なら６ギガFLOPSあるから プレステ２でやれば余裕じゃないの？
それとも精度とかの問題があるの？
0188デフォルトの名無しさん
NGNG土日通して計算しても10億までしかやってなかった
ｳﾂだ。
メルセンヌ素数やルーカステストについて
詳しいページとかない？
数学の教養もなく頭の悪いｵﾚには理解不能な事ばかり。
ちょっくらやさしい教科書でもみて
勉強したくなった。
0189デフォルトの名無しさん
NGNGフルイ方式なら
>>128 なら３０分程で終わるぞ。（ただしHDに書いた場合)
>>159のコードでも40億〜プラス１万の素数とか計算したら数秒だ
メルセンヌ素数とルカステストについては
>>167-168
さらに
>>175-176 がルカステストのサンプルも書いてるじゃないか？
まとめると、
１） 2^n-1 つまり、 ２進で１が続く数でかつ素数ならメルセンヌ素数と呼ぶ
２） nが素数でなければ2^n^-1も素数では無い
３）ルカステストは メルセンヌ素数を見つけるのに使う。素数nに対して
a=0;を初期値にして a:=(a^2-2)mod(2^n-1) をn度繰り返して途中でゼロになれば
2^n-1はメルセンヌ素数
４）a mod 2^n-1 は aを２進で表示して 2^n 以上の桁を 下位に回して足せばいい
５） -a mod b = b-a mod b
６） a^2 mod b= (b-a)^2 mod b
7）a*b は それぞれの桁数をan,bnとしたらan+bn-1桁以上の複素フーリエ変換して
求めた各要素を２乗した後、逆フーリエ変換すればいい
0190デフォルトの名無しさん
NGNGそうかーー、HDに書けばよかったか。
なんでもメモリ上にて計算するしか頭になかった。
それでメモリ壊れてOS落ちて悲しくなっていたのだ
HPのリンク先も見てませんでした
粗忽者です。
0191デフォルトの名無しさん
NGNG次の賞金は１千万桁で１０万ドルか・・・・・１千万円だな。
これを狙うとすると、
２進で３千３百万桁 ルカステストをするにはこの回数ループしなきゃいけないから
1日で計算出来る為には１秒におよそ４００回のループ（２乗と１回の引算）が出来ない
といけない。
16bitの内部表現として 2M WORD 2^21のFFTする事になる。
結局、プレステ２のFLOPSで64bit浮動小数点が出来ないと１日で１つの検査も難しい
0192デフォルトの名無しさん
NGNG>2^n^-1
これが顔文字に見えた俺は逝ってよしですか？
0193デフォルトの名無しさん
NGNG33219281
33219283
33219313
33219317
33219323
33219341
33219359
33219367
33219379 こんな調子であるぞ。
問題は、誰かがFFTを書いて数日で１つ検出出来るようになったとしても
それを分散してやった時に、賞金の配分をどうするかって所にもありそうだな
1万人が参加して一人千円じゃ配分の手間が大変だし、
宝くじ式にタマタマあてた奴が5割、ソフト書いた奴が２割、３割はサーバを提供
した奴ってのはどうだ？
0194デフォルトの名無しさん
NGNG１） メールアドレスとニックネームを登録してアクセスすると検査すべき素数Nが提供される
１０個くらい候補が出て選べるといいな
２）サーバーは登録された素数について検査中のマークを入れる
３）アプリは 検査が終わったらその値をサーバにマークする
４）ビンゴが報告されたらその値とニックネームが公開され、追試報告が何人か集まれば
それをメルセンヌ素数候補として認める
一番の課題はイタズラビンゴ報告をどう防ぐかかもな
0195デフォルトの名無しさん
NGNG0196デフォルトの名無しさん
NGNGhttp://www.nara-edu.ac.jp/~asait/prime.htm#section2
> Lucas-Lehmer テスト
> Mp = 2^p - 1 (p > 2) に対して次のような数列を考える。
> 4, 4^2 - 2 = 14, 14^2 - 2 = 194, 194^2 - 2 = 37634
> このとき Mp が p-1 番目の項を割りきるときに限り Mp は素数である。
これってどういう意味？？
全然わからない。
やさしく教えてください。
0197デフォルトの名無しさん
NGNG7 = 2^3 -1 は (3-1)番目の項 14 を割り切るので素数。
15 = 2^4-1 は (4-1)番目の項 194 を割り切らないので素数でない。
31 = 2^5-1 は (5-1)番目の項 37634 を割り切るので素数。
という判定法。
0198デフォルトの名無しさん
NGNGわかってみると。すっごーい。メルセンヌさんも>>197さんもえっらーい！
ありがとう〜♥
0199デフォルトの名無しさん
NGNGttp://homepage2.nifty.com/m_kamada/docs/fftmul.htm
ここを参考にして考えてみる。
ただし俺は飽きっぽいし、なにより割り込みが入ったら途中で放棄するかもしれん
0200199
NGNGFFTは基数４の方が若干高速だけどメンドクサイから基数２で我慢しよう
実数->FFT->複素２乗->逆FFTという事は
周波数間引型でFFTして、複素２乗した後、時間間引型で逆FFT
すれば、途中、ビット逆順にする工程を省ける筈だ
周波数間引きと時間間引き
ttp://momonga.t.u-tokyo.ac.jp/~ooura/fftman/ftmn1_2.html
0201199
NGNGを使っていた。
しかし A^2を計算させるのだから、実数の為の N/2 FFTも使える
ttp://momonga.t.u-tokyo.ac.jp/~ooura/fftman/ftmn2_1.htm
ところで、データは16bit FFTはn=21という事は16+21=37bitの
計算途中の精度が必要だ。これはdoubleの範囲に収まる。
しかし、最後の２乗はどうだろう？
(a+bi)^2=a^2-b^2 + 2*a*bi の計算で aが最大桁、 bが１桁大
なんて事になった場合a^2-b^2は精度を失う
これは問題にならないのだろうか？
0202デフォルトの名無しさん
NGNGそこだけ倍精度つかいな。
FFT の場合、最後に整数に戻る。ということを最大限に流用して
演算精度を考慮しないと、結局遅くなるよ。
いかに計算の手を抜けるか、を理論的に抑える努力が必要。
で、そこら辺は門外不出の秘伝に近いものがある。
Super PI のコードが公開されないのもそういった理由だと思ったが。
0203199
NGNGせいぜいまともにやった時の7割にしかならないらしい。
そこで、2 組の実数列の同時FFTを使ってそれぞれの２乗を求める
事を考える。
実数x,虚数部yのFFT結果の実部をA,虚部をBとして
x->FFT=( (A[n]+A[N-n]) + i*(B[n]-B[N-n]) )/2
y->FFT=( (B[n]+B[N-n]) + i*(A[n]-A[N-n]) )/2
から
4 *x^2->FFT= ((A[n]+A[N-n])^2-(B[n]-B[N-n])^2)+2i*((A[n]+A[N-n])*(B[n]-B[N-n]))
4i*y^2->FFT=i((B[n]+B[N-n])^2-(A[n]-A[N-n])^2)-2 *((B[n]+B[N-n])*(A[n]-A[N-n]))
と求められる。
つまりこの２式を加算して逆FFTすれば 実部と虚部に ２乗x4倍値として求まる
（どうせFFT->逆FFTの後 2^n=Nで割る必要があるので4Nで割ればいい)
0204199
NGNG・・・・でもやっぱり16bitの２乗は32bitだし
N桁の２乗は中間桁ではN倍されるから32+21=53bit必要だから
やっぱりlong doubleでも ダメっぽいね
8bitにするしかないな・・・結構厳しいなあ
0205199
NGNGEbitの所で分割する。E=(K+1)/2とする
上位をx下位をyとする
(x*2^E+y)^2
=x^2*2^2E+x*y*2^(E+1)+x^2
=x^2*2^(K+1)+x*y*2^( E+1)+y^2
最初の x^2*2^(K+1) はモジュラー計算すると x^2*2
=x^2*2+y^2 + x*y*2^( E+1)
で、xとyを >> 203 の方法で ２実数のFFTを使い
周波数軸上で X^2,Y^2, X*Yを求める
から、2*X^2+Y^2を 周波数軸で加算し実数部で得
られるようにし X*Yを虚部に得られるように調整
した後逆FFTする
後は正規化して加算すればいい。
ただし、FFTする時は x,yは8bit単位に分割し
計算精度は long doubleでやらなければならない
0206199
NGNGやっぱり速度上違いが出るのだろうか？
doubleが早いなら FFTはdoubleで計算して ２乗後は
Extendedで計算ればいいんだけど
0207199
NGNGまず、ビット反転ルーチン 今回は必要無いかと思ったんですが
実数２つのFFT上で掛算をする時に必要だったので
procedure MkBitRevTbl(n:integer;var tbl:array of integer);
var i,j,k,mh,m:integer;
begin
m :=1 shl n;
mh:=m div 2; //mh=msbのみ1のデータ
i:=0; // 1単位に増えるアドレス
j:=0; // ビット反転アドレス
while true do
begin
tbl[i]:=j;
inc(i);if i>=m then break;
k:=mh; //mh=msbのみ1のデータ
while k<=j do begin //上のビットが立っていたら
j:= j and (not k); //それを落としてゆく i-kでも良い筈
k:=k div 2;
end;
j:= j or k;
end;
end;
0208199
NGNGvar i,j,k,jj,nn,mh,m,h:integer;
var w,wr,wi,wr0,wi0,ww:Extended;
var dr,di:Extended;
begin
m:=1 shl N;
// 最初の１回目
w:=2*PI/m;
mh:=m div 2;
h:=1;
wr0:= cos(w); wr:=1;
wi0:= sin(w); wi:=0;
for i:=0 to mh-1 do begin {最初のループは色々高速化出来るので外に出す}
dr := x[i] - x[i+mh]; { 上位は必ずゼロなので引算省略可能 }
di := y[i] - y[i+mh];
x[i]:= x[i] + x[i+mh]; { 上位は必ずゼロなので この２行不要 }
y[i]:= y[i] + y[i+mh];
x[i+mh]:= wr * dr - wi * di;
y[i+mh]:= wr * di + wi * dr;
ww:= wr0 * wr - wi0 * wi; { sinテーブルをひくのとどちらが早い？ }
wi:= wr0 * wi + wi0 * wr;wr:=ww;
end;
for nn:=2 to N do begin { TODO : N-1,N はもっと効率をあげられる }
h:=h*2;
m:=mh;
mh:=mh div 2;
w:=2*PI/m;
wr0:= cos(w); wr:=1;
wi0:= sin(w); wi:=0;
for i:=0 to mh-1 do begin
j:=i;
for jj:=0 to h-1 do begin
k:=j+mh;
dr := x[j] - x[k];
di := y[j] - y[k];
x[j]:= x[j] + x[k];
y[j]:= y[j] + y[k];
x[k]:= wr * dr - wi * di;
y[k]:= wi * dr + wr * di;
inc(j,m);
end;
ww:= wr0 * wr - wi0 * wi; { sin/cosを計算させる }
wi:= wr0 * wi + wi0 * wr;wr:=ww;
end;{end i}
end; {end nn}
end;
0209199
NGNGprocedure ifft(n:integer;var x,y:array of Extended);
var dtN:integer;
var i,j,k,ii,jj,nn,mh,m,h:integer;
var w,wr,wi,wr0,wi0,ww:Extended;
var dr,di:Extended;
begin
mh:=1 shl (n-1);
// 最初の１回目
for ii:=0 to mh-1 do begin {最初のループは高速化出来るので外に出す}
i:=ii*2;
dr := x[i] - x[i+1];
di := y[i] - y[i+1];
x[i]:= x[i] + x[i+1]; x[i+1]:= dr;
y[i]:= y[i] + y[i+1]; y[i+1]:= di;
end;
h:=2;
mh:=mh div 2;
for ii:=0 to mh-1 do begin {次のループも高速化出来る}
i:=ii*2*h;j:=i+h;
dr:= x[j] ; {e^0}
di:= + y[j];
x[j]:= x[i] - dr;
y[j]:= y[i] - di;
x[i]:= x[i] + dr;
y[i]:= y[i] + di;
inc(i);inc(j);
dr:= + y[j]; {e^-2πi/4}
di:= -x[j] ;
x[j]:= x[i] - dr;
y[j]:= y[i] - di;
x[i]:= x[i] + dr;
y[i]:= y[i] + di;
end;
for nn:=3 to N do begin
h:=h*2;
mh:=mh div 2;
w:=PI / h ;
wr0:= cos(w);
wi0:=-sin(w);//逆変換用
for ii:=0 to mh-1 do begin
wr:=1;
wi:=0;
i:=ii*2*h; j:=i+h;
for k:=0 to h-1 do begin
dr:= wr*x[j] - wi*y[j];
di:= wi*x[j] + wr*y[j];
x[j]:= x[i] - dr;
y[j]:= y[i] - di;
x[i]:= x[i] + dr;
y[i]:= y[i] + di;
ww:= wr0 * wr - wi0 * wi; { sin/cosを計算させる }
wi:= wr0 * wi + wi0 * wr;wr:=ww;
inc(i);inc(j);
end;
end;
end;
end;
0210199
NGNGprocedure Mul_x_y(n:integer;var x,y:array of Extended);
var i,j,k,m,mh:integer;
var ar,ai,br,bi,cr,ci,C:Extended;
begin
fft(n,x,y); //FFTして
m:=(1 shl n);
C:=1/4/m;
mh:=m div 2;
for i:=1 to mh -1 do
begin
j:=BitRevTbl[i];
k:=BitRevTbl[m-i];
ar:=x[j]+x[k]; //それぞれの係数を求めて
ai:=y[j]-y[k];
br:=y[j]+y[k];
bi:=x[k]-x[j];
cr:= (ar*br-ai*bi)*C; //a*b を計算した時の実部
ci:= (ar*bi+ai*br)*C; //a*b を計算した時の虚部
x[j]:= cr ;
y[j]:= ci ;
x[k]:= cr ;
y[k]:= -ci ;
end;
x[ 0]:= 4* x[ 0]*y[ 0]*C; //a*b を計算した時の実部
y[ 0]:= 0;
x[ 1]:= 4* x[ 1]*y[ 1]*C;
y[ 1]:= 0;
ifft(n,x,y);
end;
0211デフォルトの名無しさん
NGNG0212デフォルトの名無しさん
NGNG19 17 13 11 7 5 3 = n
----------------------------------
524287 131071 8191 2047 127 31 7 = 2^n-1
----------------------------------
4 4 4 4 4 4 4 =a1
14 14 14 14 14 14 0 =a2=(a1^2-2)%(2^n-1)
194 194 194 194 67 8
37634 37634 4870 788 42 0 =a4
218767 95799 3953 701 111
510066 119121 5970 119 0
386344 66179 1857 1877
323156 53645 36 240
218526 122218 1294 282
504140 126220 3470 1736
103469 70490 128
417706 69559 0
307417 99585
382989 78221
275842 130559
85226 0
523263
0
0213デフォルトの名無しさん
NGNG×y->FFT=( (B[n]+B[N-n]) + i*(A[n]-A[N-n]) )/2
○y->FFT=( (B[n]+B[N-n]) - i*(A[n]-A[N-n]) )/2
0214199
NGNGとりあえず、必要だったルーチンを掲載
FFTをする為に2のベキで収まる数を求める関数
function iGTexp2(k:integer):integer; //k < 2^n となる 2^n
begin
k:=k or (k shr 16); //ビットシフトしながらORすると
k:=k or (k shr 8); //下位ビットが全部1になる
k:=k or (k shr 4);
k:=k or (k shr 2);
k:=k or (k shr 1);
Result:=k+1; //１を足せば 2^nになる
end;
function iLog2(m:integer):integer; //m=2^n を与えて nを返す
begin
if m>=$10000 then Result:=iLog2(m shr 16)+16
else if m>=$100 then Result:=iLog2(m shr 8)+8
else if m>=$10 then Result:=iLog2(m shr 4)+4
else if m>=4 then Result:=iLog2(m shr 2)+2
else if m>=2 then Result:=iLog2(m shr 1)+1
else Result:=0; //1か0 0ならエラーだけど知らない
end;
0215デフォルトの名無しさん
NGNGエラスト？エストラ？
そのフルイの配列確保がオンメモリでは出来なさそうなので
ＨＤにアクセスしながら計算させたいのですが
ＨＤ直接(?)アクセスする方法がわかりません。
教えてもらえないでしょうか。
環境は...Delなんですが...
0216199
NGNGHDにアクセスと難しく考えなくても そのまま
メモリを array of byte で取って SetLengthすれば自動的に
仮想記憶が働いてハードディスクをアクセスするよ
>>205 で書いた方法で
E=(K+1)/2とせずに、 8の倍数にすれば正規化処理で少し
楽出来るけど、計算精度が 最悪9bit余計に必要になるので
Extended の限界を超えてしまう。
0217199
NGNGvar a,b:array of Extended; //aが下位0〜e-1bit bは上位n-1〜e
var e:integer; //分割位置
var k8:integer; //aのmsbのある位置
var m,n:integer; //k < m=2^n となる n,m
var be,beb,bmo,ae,aeb,amo:integer;
var cy:int64; //溢れ分
//その初期化
procedure init;
begin
e:=(k+1) div 2;
m:=iGTexp2(k);
n:=iLog2(m);
k8:=k div 8;
SetLength(a,m);
SetLength(b,m);
be := 1+((k-e) div 8);//bのmsb+1を含むバイト位置の次
beb:=1 shl (e+be*8-k); //aから見てbeの位置は何倍されているか
bmo:=1 shl (k-e +8-be*8); //msb+1 のビット位置
ae :=1+((e) div 8); //aの bのlsb位置を含むバイト位置の次
aeb:=1 shl ( ae*8-e); //bから見てaeの位置は何倍されているか
amo:=1 shl (e +8-ae*8); //a[ae]における bのlsb のビット位置
end;
0218199
NGNGFFTによる掛算結果をそのまま再度２乗するとデータ精度が保て
ないので、各桁を0〜255程度に、さらに上位がゼロになるように
正規化する
ここでは a + b *2^e を保存するように変換する
ここで x*2^k=x となるので
1. a[i] + b[i+be]* 2^(e+be*8-k) ここで 2^(e+be*8-k)=beb
2. a[i+ae]*2^(ae*8-e) + b[i ] 2^( ae*8-e)=aeb
1を使って 0〜(e-1)bit をaに
2を使って e〜(k-1)bitを bに収める
procedure Normal1;
var i,w:integer;
var cym:int64;
begin
cy:=round(b[be-1]); //bのmsb+1 を含むバイト
b[be-1]:=cy mod bmo; cy:=cy div bmo;
for i:=be to k8 do begin
cy:=cy+round(a[i-be]+b[i]*beb) ; b[i]:=0;
a[i-be]:=cy and 255; cy:= cy shr 8;
end;
w:=( k8-be+1 - (ae-1) ) ;//cy のオフセット
for i:=0 to w-1 do cy:=cy shl 8;
cym:=1;
for i:=ae-1 to High(a) do begin
if cym<$1000000 then begin;// そんな上には溢れていない筈
cy:=cy +cym*round( a[i] ); //bのmsb+1 を含むバイト
cym:=cym*256;
end;
a[i]:=0;//上位に計算誤差を貯めないように完全にゼロにする
end;
a[ae-1]:=cy mod amo; cy:=cy div amo; //aのmsb+1 を含むバイト
for i:=ae to k8 do begin
cy:=cy+round(b[i-ae]+a[i]*aeb) ; a[i]:=0;
b[i-ae]:=cy and 255; cy:= cy shr 8;
end;
w:=( k8-ae+1 - (be-1) ) ;//cy のオフセット 何か怪しい？
for i:=0 to w-1 do cy:=cy shl 8;
cym:=1;
for i:=be-1 to High(a)-ae+1 do begin
if cym<$1000000 then begin;// そんな上には溢れていない筈
cy:=cy +cym*round(b[i]); //bのmsb+1 を含むバイト
cym:=cym*256;
end;
b[i]:=0;////上位に計算誤差を貯めないように完全にゼロにする
end;
b[be-1]:=cy mod bmo; cy:=cy div bmo;
end;
このcyから-2を引いてa[0]に足し込めばいい
（負数のままでも2乗は正常に働く筈）
0219199
NGNGE bitの所で分割する。E=(K+1)/2とする
上位をx下位をyとする
(b*2^E+a)^2
=b^2*2^2E+2*a*b*2^(E)+a^2 //2E=k+1なので
=b^2*2^(K+1)+2*a*b*2^( E)+a^2
var BitRevTbl:array of integer;
procedure SquareFFTMod(n:integer;var x,y:array of Extended);
var m:integer;
var i,j,k:integer;
var ar,ai:Extended;
var br,bi:Extended;
var Xr,Xi:Extended;
var Yr,Yi:Extended;
var C,C2:Extended;
begin
m:=(1 shl n);
if Length(BitRevTbl)<>m+1 then
begin
SetLength(BitRevTbl,m+1);
MkBitRevTbl(n,BitRevTbl);
BitRevTbl[m]:=0;
end;
FFT(n,x,y);
C :=1/4/m;
C2:=2/4/m;
for i:=0 to (m div 2) do
begin
j:=BitRevTbl[i];
k:=BitRevTbl[m-i];
ar:= x[j]+x[k]; //x->FFT=( (A[n]+A[N-n]) + i*(B[n]-B[N-n]) )/2
ai:= y[j]-y[k];
br:= y[j]+y[k]; //y->FFT=( (B[n]+B[N-n]) - i*(A[n]-A[N-n]) )/2
bi:=-x[j]+x[k];
//X= a^2+2*b^2 Y= a*b*2
Xr:=C *((ar*ar-ai*ai)+ 2*(br*br-bi*bi));
Xi:=C2*( ( ar*ai ) + 2*( br*bi ) );
Yr:=C2*(ar*br-ai*bi); //2*a*b
Yi:=C2*(ar*bi+ai*br);
x[j]:= Xr -Yi ; //Xを実部 Yが虚部になるよう配置する
y[j]:= Xi +Yr ;
x[k]:= Xr +Yi ; //j==kなら重複するが同値になる筈
y[k]:= -Xi +Yr ;
end;
iFFT(n,x,y);
end;
0220デフォルトの名無しさん
NGNG検査すべき素数を提供する時にパスワードを生成して渡しておくんじゃダメ？
そのまま検査された値全部にフラグとパスワードを保存するとディスクの無駄になるけど。
0221デフォルトの名無しさん
NGNG>>216さん
その確保したarray of byteの
Xビット目が0なのか1なのか判断する方法が
わからないでございます。
どうやるのでしょう。
質問ばかりですいません。
0222デフォルトの名無しさん
NGNGX ビット位置が１かどうかは
tbl[X div 8] and (1 shl (X mod 8) )
ビットを立てるなら
tbl[X div 8] |:= tbl[X div 8] or (1 shl (X mod 8) )
少しだけ高速な方法は
tbl[X shr 3] |:= tbl[X shr 3] or (1 shl (X and 7) )
ただ、Xが偶数なら素数じゃないから
予め X mod 2 =0 かどうかで判定して除けばテーブルは半分で済む
0223199
NGNG速度的に面の見積もりは出来た
k=2281付近で２５秒
k=4423付近で120秒
世界記録を狙うとなると１個の判定に１年ですまない計算だ
２桁以上計算速度を上げるアイデアが必要
0224199
NGNG○ n:=iLog2((m+7)div 8);
という事で 上の数字は1/8になった
k=2281付近で２秒
k=4423付近で8.5秒
k=9689付近で42秒
他にも細かいミスを修正した結果、とりあえずこの付近では
ttp://www.geocities.co.jp/Technopolis/1793/mersennu.html
の表の通りの判定が出来るようになった
汚いソースになったけど、誰かいるなら公開するけど
0226199
NGNG2ヶ月程で１つの数の判定が出来るだろう。
ただ、調べなければいけない数も膨大なので
出来れば1/50 せめて1/10に縮めないと・・・
0227デフォルトの名無しさん
NGNG100堕慰くらいで実行スレばいいよ。
そんなに2chプログラム人口いないか？
0228199
NGNG世界記録は22bit付近だから ワード分割サイズをAbitとしたら
2A+22bitでいいならDoubleで足りる
2(A+22)bit 必要ならExtendedが必要になる
0229199
NGNG1bit単位に分割してFFTで掛算すれば何ビットの精度が必要だろう？
もしそれが22bitでいいのなら、PS2で実行という方法もある
でもメモリが 2*4*4Mbyte必要になるけど
0230デフォルトの名無しさん
NGNGありがとー。
>ビットを立てるなら
>tbl[X div 8] := tbl[X div 8] or (1 shl (X mod 8) )
ビットを落とすのはどうしたらいいのかなあ。
論理演算とか苦手だわ
....notしてからandするのかな？
0231デフォルトの名無しさん
NGNGビットを落とすのは not して andでもいいし
すでに立っているのが判ってれば xorでもいいよ
0232デフォルトの名無しさん
NGNG>tbl[X div 8] := tbl[X div 8] or (1 shl (X mod 8) )
の(1 shl (X mod 8) )
が00010000としたら
左から4bit目を立てる為にorしてるんですよね
するってとー
>ビットを落とすのは not して andでもいいし
tbl[X div 8] := tbl[X div 8] and (not (1 shl (X mod 8)) )
つうことですか？
0233デフォルトの名無しさん
NGNGtbl[X div 8] := not ( (not tbl[X div 8]) or (1 shl (X mod 8)) )
tbl[X div 8] := tbl[X div 8] xor ((not tbl[X div 8]) and (1 shl (X mod 8)) )
さらに、setもresetも出来る方法(ビットの代入)として
flg=0 か $ffとして
tbl[X div 8] := tbl[X div 8] xor (( ftbl[X div 8] xor flg ) and (1 shl (X mod 8)) )
0234199
NGNGSingle = 32bit 仮数部23bitにして計算させてみた。
Doubleの結果と比較して k=251の時に差が出た
これから浮動小数点ならほぼ 2*bit幅+log2(k)が限界
と考えていいように思う。
つまりdouble (仮数部52bit) のままで計算可能ではないかと
であれば、>>224の値はほぼ半分になる
0235デフォルトの名無しさん
NGNGIndexビット目が0か1かを調べるのって
if (Tbl[Index div 8] and (1 shl (Index mod 8) )) = (1 shl (Index mod 8) ) then
Result := True
else
Result := False;
こうするんですね。
>>222の内容からちょっと悩みました
おかげでどうにか出来そうです。
なんでかよくわからないですが
TBits.Count := High(Integer)
するとメモリオーバーするみたいで
TByteArray: array of Byte
SetLength(TByteArray, ((High(Integer)-1) div 8) +1 );
すると大丈夫みたいっす。
メモリ確保量は同じ気がするのに、なぜだろう。
世界記録。がんばってくださいマンセー
トロイけど、このスレ呼んで勉強しーとこ。
0236199
NGNGResult := (Tbl[Index div 8 ] and (1 shl (Index mod 8) )) <>0 ;
さらに
Result := (Tbl[Index shr 3 ] and (1 shl (Index and 7) )) <>0 ;
0237デフォルトの名無しさん
NGNG0238199
NGNG>>71の賞金を貰うならさらに １００倍かかるな
100倍になると、今のままだと8年はひとつの数の判定にかかる計算
それにメモリもオンメモリで出来る限界付近だからもっと速度は落ちるかも
やっぱり1千万の賞金は簡単じゃないな
0239デフォルトの名無しさん
NGNG並のパソでそれだけいくと
もちっと速いコンピつかえば＆
数年でクロック数上がれば、どんどん賞金に知被いて....
素数の歴史上、大学のコンピを使えた高校生が
最大素数を求めたことがあったらしいけど
アルゴリズムだけOKなら大型コンピュータにて
流すことはできませんか？
0240デフォルトの名無しさん
NGNG条件は同じだから、差をつけるには やっぱりプログラム技術の差でつけないとね
それに、スパコン借りるなら成果出さないといけないだろう。
今の速度じゃ、スパコン１台では見つかるとはとても思えない
計算速度を上げるには、やっぱりFFTの固定小数点化しか
ないんじゃないかな。固定少数点なら、今のCPUなら足算は２〜３並列化
されるから巧くコーデングされれば相当高速になる筈。
データを64bit固定少数点にして、サインのテーブルは32bitに
すれば、掛算のコストは倍、たし算は３倍になるけど、
16bit単位に区切れるから計算量としては同程度
それにしても、なんで神様は 2^16+1を素数にしたのに 2^32+1とか2^64+1を素数に
しかなったんだろ
0241デフォルトの名無しさん
NGNG128bitSIMDにアセンブラレベルで対応した方が速いんじゃないかな
１クロックで４並列演算とか出来るのは、あれは単精度(32bit)の場合だっけ？
0242デフォルトの名無しさん
NGNGAMD 3DNow
http://www.amd.com/japan/K6/k6docs/j21928c.pdf
Pentium 2,3
ftp://download.intel.co.jp/jp/developer/jpdoc/iaopt_j.pdf
0243199
NGNG１回の計算 は 15/12112秒 １秒に750回計算出来ている。
１回の計算は FFT+２乗+逆FFT+２分割正規化処理
fftのサイズは 2048 = 11ステージ でFFT+逆FFTをやってるから
１回のFFTの複素積和は 750*2048*11*2=33.8M/秒
複素積和は4回の掛算と2回の加算だから
今現在、200MFLOPS 付近が出てる模様
メモリサイズは 2048*4*2=16Kbyteだからデータキャッシュはヒットしてる筈
この数字はathron800Mでの話
たしかに倍精度で2GFlopsとか出れば一挙に10倍になるんだから
P4は魅力的だな・・・でもまだ買い換え予定はないし
出来れば整数演算でどれくらいでるか試してみたいな
0244199
NGNGfor n:=1 to N do
for i:=0 to 2^(N-n) -1 do
for j:=0 to 2^(n-1) -1 do
begin
c:=2^(N-n);
k:=i+j*2^(N-n);
x :=a[k]-a[k+c];
a[k ]:=a[k]+a[k+c];
a[k+c]:=x*E(i,nn);
end;
こういう順にやっている。このままだと Nが大きくなると
キャッシュミスが大きくなると思う。そこで、
最後の j の単位毎に終わったら 次の n を先に計算させ
るように計算順番を入れ替えたい。つまり j,n,iの順にネスト
させたい。どうループを書けばいいんだろ？
再帰を使わないと書けないかな？
0245199
NGNGttp://www.pluto.ai.kyutech.ac.jp/plt/matumoto/fftb_export.c.txt
キャッシュを考慮してバタフライ演算が入りきるサイズに分解してるようです
0246199
NGNGCのfor文みたいな柔軟なループが書き難いし、 Athronはcallは
殆ど負荷にならないらしいから再帰でいいかなと
procedure ifft0(nc,st,stp:integer); // キャッシュをヒットさせる為に再帰を使って計算
var stph:integer;
begin
if nc> 0 then
begin
stph:=stp shr 1;
ifft0(nc-1,st ,stph); //２分割しながら深くネストして
ifft0(nc-1,st+stp ,stph);
ifft1(nc , st ); //一番奥から順に実行する
end;
end;
begin
ifft0(n,0,1 shl(n-1)); //呼び出しの根っこ
end;
fftの方は逆に１ブロック実行してはネストしてゆく格好で書けました
0247199
NGNG整数化する為に、
１）sin,cosの計算精度を向上させる為順計算でなく
1/2に分解して再帰する手順に要変更
しかし、問題はm=2^16を超える頃にcosはほぼ1になる場
合がある点。
これを解決するには、
２）cosではなくcc=1-cosを記録して計算時には
x*cos=x*(1-cc)=x-x*ccとする
３）基数2のFFTだとsin側も同じく1になるので
基数8に分解し sinが45°以下になるようにする
計算精度を上げる為に、sin,cosの少数点位置はMSBの左
にする。
データには負数があるので符号無し掛算命令を
利用する為、a*b のaが負数の時は(a-$100000000)*b
a*b-b<<32と処理しなければいけない。判断は遅いので
a*b-(b&(a>>32))<<32
最初のステージで少数点位置を Hだけシフト
(a*b-(b&(a>>32))<<32)>>(32-(H))
１回のステージ毎に 少数点位置は丸め後右に動かす
((a*b-(b&(a>>32))<<32)+$7fffffff)>>31
FFT終了時には 少数点はH-Nにある
２乗計算後に2(H-N)-Nだけシフトし、少数点位置をNにする
逆FFTも１回のステージ毎に 少数点位置は丸め後右に動かす
((a*b-(b&(a>>32))<<32)+$7fffffff)>>31
逆FFT終了時には整数化されている
と考えた段階でゲップ状態・・・・メンドクサイと思ったら終わりかな
0248199
NGNG○$3fffffff
0249デフォルトの名無しさん
NGNGそれとも実数でやるのがミソなんでしょうか。
ps2linux でやってみようと思って考えてたんだけど、
vu には 32bitx32bit の乗算なんてないのか(当然か...)
12x12 ならできるけど、EE はこんなのにはむいてないみたい。
今更だけど。
0250199
NGNGttp://www.geocities.co.jp/WallStreet/6351/mersenne.html
勉強してみます。
理論が良く判ってないのと、
フーリエ変換が他への応用も出来るのと浮動小数点を使う
限りは簡単だから拘っていました。
でもFFTを固定少数点でやるくらいならそっちがマシなような気もします
0251199
NGNG2^32-1のような 2^n-1 を 法にするのはどうでしょう？
2^n-1を法にすれば n乗根は必ず2=2^0になるので計算し易くなります。
というのは全ての係数が2のベキ乗になるからです
例えば2^8-1=255を法にすれば 2^i のiを表にして縦軸をビット逆順に入れ替えると
a0: 1 2 3 4 5 6 7 0 注意:0は2^n=1 の事
a4: 5 2 7 4 1 6 3 0
a2: 3 6 1 4 7 2 5 0
a6: 7 6 5 4 3 2 1 0
a1: 2 4 6 0 2 4 6 0
a5: 6 4 2 0 6 4 2 0
a3: 4 0 4 0 4 0 4 0
a7: 0 0 0 0 0 0 0 0
と非常に簡単に計算出来ますよね？
0252249
NGNGn 乗根を考えるということは
もとの数を n word に分けるということなので、
n がある程度大きくないと FFT のメリットがないような気がするのですが
どうでしょう。
あと、素数を法とした方が簡単な気がするんですけど、
素人なのであんまり自信はないです。
で、とりあえず、n が 2 の巾だと 2 分割の FFT ができて簡単なので
そういうのを探してみました。
2^N-1 の判定をする場合、
1 word = w bit, p を法として考えるとして
n*w >= 2N, p > n*2^{2w} であれば安心。
N を 2^23 より小さいとすると、
w=2, p=2^23*(2^4+4)+1(素数), n=2^23 とすると
32 bit の乗算でなんとかなる。
w=16 だと p=2^20*(2^32+2)+1, n=2^20 として、
64 bit の乗算が必要だけど繰返しは少ない。
なかなかつらいけどこんなもんでしょうか。
p=65537 で適当に作ってみたらものすごく遅かった。
0253199
NGNG60*60*24*365=31,536,000なので、およそ１回の２乗にかかる秒数が１つのルカス
テストにかかる年数になります。
たとえば 2^16384 -1 を法にして 16384乗根は 2なので
FFTにした時の各段階の係数は全て 2のベキになります
(8170bit単位に 8192個の単位に分割すれば =66Mbit迄可)
２のベキの係数なので回転操作と加算操作で計算出来ます
(掛算ではないので何ビットあっても不利ではない)
結果、16384bitのの２乗の計算が16384個必要になるけど、
これも 2^256-1の法で 64bit単位に256個に分解すれば
256bitの掛算は(32bit掛算*８＾２=)64になります。
掛算回数はわずか 64*256*16384 = 268M
ただしメモリを256Mbyte使っているので、FFTにかかるオーダは
最初の16384乗根で 14 つぎの256乗根で8 で合計14+8=22
一回のステージは２回の加算と１回のシフト、逆変換を含めて
256M/4*22*3*2 = およそ 8.5G という事で
8ギガmipsのマシンでもルカステスト１回にやはり年単位は必要
という事になってしまいます
これを短縮するには、この変換操作を数１００台のマシンで分解
して計算させるアイデアがやっぱり必要になるようですね
0254199
NGNG>n がある程度大きくないと FFT のメリットがないような気がするのですが
たしかに、2^32+1で2^24乗根なんてのがあれば最高ですよね
>あと、素数を法とした方が簡単な気がするんですけど
たしか、2^n+1が素数であれば 2^n乗根が2になるというメリットは
あるけど、それ以外にはあんまりないんじゃないかな
2^(2^n)-1 を法にすれば、 2^n乗根は 2になる筈です
(X+1)mod X＝１だから、 2^(2^n) mod(2^(2^n)-1) =1
だからFFTの全部の係数は 2のベキ乗になります。
という事は、FFT中掛算は実質必要無いという事です
たとえばFFTの初段は
a[i]+a[i+2^(n-1)] , ( a[i]+a[i+i+2^(n-1)]*2^(2^(n-1)))*2^i
です、これバイト位置をあわせて単純に加算して、最後にiビットシフトするだけです
2^m-1 のモジュラー計算は、桁溢れがあればそのまま下位に足すだけなので
2^m+1の場合よりもさらに簡単です(負数を考えなくていいから)
と考えたのだけど、ホントかどうかはまだ試していません
0255249
NGNGちょっと考えてみたんですけど、
素数でないとやっぱり苦しいんではないでしょうか。
フーリエ変換と逆変換が実際に逆になることを示すのに、
a の p を法とした位数を k として、0<i<k ならば
(a^0+a^i+...+a^{i(k-1)}) mod p = 0 となることを使いますけど、
たとえば p=2^4-1, a=2 とすると k は 4,
i=2 として a^0+a^2+a^4+a^6 mod p = 10 となってます。
何かうまくやる方法があるんでしょうか ?
0256199
NGNG15を法にしてみると FFTは
|1,1,0,0| |1,0,1,0|
|1,4,0,0| |0,1,0,1|
|0,0,1,1| |1,0,4,0|
|0,0,1,4| |0,2,0,8|
を２つ掛ければ出来るけど、どっちも逆変換しようとすると 3が出て来て
3の逆数が無いから解けないですね
やっぱりダメか・・・・
0257199
NGNG|4,-1,0,0| |8,0,-2,0|
|-1,1,0,0| |0,8,0,-1|
|0,0,4,-1| |-2,0,2,0|
|0,0,-1,1| |0,-2,0,1|
とすると３倍、６倍の結果になってトータル 18%15 = 3倍になるけど、いちおう計算出来ますね
2^8-1の場合も 15倍になる解なら同じ規則で見つかりました
ただ、この調子だと実際に計算出来るビット数はだいぶ小さくなってしまいますね
0258199
NGNG2^(n/2)-1 倍された逆FFTならできるという事になります
32bitの場合は 65535倍される訳で、これは約数になっていますから
つまり半分のビットが使えないという事みたいですね
0259199
NGNG11*101 =1111
1111 * 10001=11111111
11111111 * 100000001= て調子で 2^m+1 との掛算なんだ
0260デフォルトの名無しさん
NGNG0261デフォルトの名無しさん
NGNG常識的な速度で出す事やOnメモリでダス方法はわかったのだけど
それをファイルに落とす場合に
for i=1 to High(Integer)
IF iが素数なら
S := S + IntToStr(i)
のIntToStrが超低速で参ってます。
なんとか高速に沢山の数値を
IntToStrと同様に文字列にしてしまう方法はないでしょうか？
ファイルの保存形式が別にStringじゃなくても
いいのでIntegerのメモリイメージがそのままファイルに
落ちてくれるといいのだけど...
0262デフォルトの名無しさん
NGNGファイルマッピングを使ったらどう？
ファイルを開いて CreateFileMappingで割り当てて
MapViewOfFile でポインタを貰って
ポインタ使ってInteger配列として操作出来るとか
0263デフォルトの名無しさん
NGNGTByteArray = array of Byte;
TSosuuByteArrayRiddleEx = class(TObject)
private
FCount: Integer;
FByteArray: TByteArray;
function GetItems(Index: Integer): Boolean;
procedure SetItems(Index: Integer; const Value: Boolean); ;
function GetCount: Integer; ;
procedure SetCount(const Value: Integer);
protected
constructor Create;
public
property Count: Integer read GetCount write SetCount;
property Items[Index: Integer]: Boolean read GetItems write SetItems;
end;
FByteArrayは>>222で教えてもらったような所有の仕方(偶数除去)
してるのですが
これをファイルマッピングするのは
TSosuuByteArrayRiddleExの先頭アドレスを指定するといいのでしょうか？
…むしろFByteArrayの先頭アドレスかな。
またはこの篩クラスから
Sosuu: array of Integer
に代入してから保存した方が使いやすそう
(要はファイルサイズの問題)
かと思うのですが
その場合のファイルマッピングの方法を
もう少し教えていただけないでしょうか。
0264デフォルトの名無しさん
NGNGさっさと分散＆ブルートフォースで解決
0265デフォルトの名無しさん
NGNG１） 動的配列じゃなくて、配列のポインタにする事
２） そのポインタの実体をGetMem とか GlobalAlloc で確保してみる
３）ファイルマッピングを使ってみる それをコンポーネントにしておくと
ちょっと便利かな
というようにバージョンアップしてみましょう
32bit程度ならまだまだ ビットイメージで持っていた方がファイルサイ
ズの面でも速度面でも有利だよ。
ただ、32bitくらいの篩なら保存するよりはその都度作った方が速いかもしれないね
0266デフォルトの名無しさん
NGNGグサッ Λ＿Λ
：（ﾟД｡ ；）
――――（”∴ ） グヒャッ・・・
．：｜｜ ｜
∴ （__（＿_）
そうか・・・初心者だったのか。
> １） 動的配列じゃなくて、配列のポインタにする事
> ２） そのポインタの実体をGetMem とか GlobalAlloc で確保してみる
> ３）ファイルマッピングを使ってみる それをコンポーネントにしておくと
> ちょっと便利かな
うーーーん
全然イメージがつかめません。
具体的なコードが見られる所ってないですか？webでもマニュアルでも
動的配列って実際はポインタじゃなかったです？
>32bit程度ならまだまだ ビットイメージで持っていた方がファイルサイ
>ズの面でも速度面でも有利だよ。
2^30bit=2^27byte=2^17KB=2^7MB=128MBは必要なんだけど
速度面はともかく
2,3,5,7,11,13,....という配列とくらべて
ファイルサイズは有利じゃないでしょ。
>32bitくらいの篩なら保存するよりはその都度作った方が速いかもしれないね
Myマシンが非力なのでメモリ的につらいの。
0267デフォルトの名無しさん
NGNG素数定理を使って32bitの素数の個数を 2^32 / log(2^32) 個と見積ると、
2^32 / log(2^32) * 4 Byte ≒ 195000000 Byte ≒ 186MB 必要。
ファイルサイズでも不利だと思われる。
0268デフォルトの名無しさん
NGNGこっちが当然の書き方だったんだけどなあ
http://www.nifty.ne.jp/forum/fdelphi/faq/00160.htm
0269デフォルトの名無しさん
NGNGn* 2^n/log(2^n)=n*2^n/(log(2)*n) = 2^n/log(2)
log(2)はおよそ 0.3だから 3*2^n BITがおよそ必要なメモリ
ビットマップで記録するのは当然 2^nBIT ただし偶数を除けば その1/2
だからビットマップで記録すれば ６倍は効率が良い
2^32/log(2^32)*4はおよそ 1.8GByte
32bitの素数が欲しい場合、16bitの篩を用意しておけばいいのでは？
0270デフォルトの名無しさん
NGNGサンキュウ。
そうか、素数定理というのが抜けてたのか。
その定理よくわからないけど、検索してみる。
だけど、Integerって最大は2^31だから
もう少し減るかも。
なんか、自分の数学的素養が著しく
失われてることに気が付いたので
http://www.geocities.co.jp/Technopolis/1505/log.htm
ここから読んでみます。
プログラム初心者というより
人生初心者か？義務教育初心者....(TДT)ｳﾂﾀﾞ
>>268
D3ってもうずいぶん前なので､､今ではそのやり方のほうが
トリッキーかも。
>>269
ごめん。全然わかんないや。
素数定理というヤツかな。
>32bitの素数が欲しい場合、16bitの篩を用意しておけばいいのでは？
そうなの？そうは思えないんだけど、
0271269
NGNGだから 1/log(2)は 1.44 だった まあどっちにしてもサイズに関係なく
単純に記録するなら 常にビットマップで記録した方がメモリ容量は少ないって
事だね。 素数は案外詰まってるんだよね
0272デフォルトの名無しさん
NGNG0273デフォルトの名無しさん
NGNGhttp://www.com.mie-u.ac.jp/~kanie/tosm/humanind/yosou.htm#prime_number
数学は科学の女王で、数論は数学の女王だ。
頑張れ。
ζと対話できるようになれば一人前だ。
0274デフォルトの名無しさん
NGNG0275デフォルトの名無しさん
NGNG無理にきまってんだろ。
0276デフォルトの名無しさん
NGNG0277デフォルトの名無しさん
NGNG百人くらいを一つのグループにして、そのグループで一つの数のチェックを
担当して、 ヒットしたらそのグループ＋運営側で半々に わけるってのは？
百人なら半分でも５万円になる
0278デフォルトの名無しさん
NGNG0279デフォルトの名無しさん
NGNG0280デフォルトの名無しさん
NGNG感動的ストーリーでも捏造するか？
0281デフォルトの名無しさん
NGNG単にマシンパワー提供する程度で金が貰える
と思うバカは必要なっしんぐ。
0282デフォルトの名無しさん
NGNG0283デフォルトの名無しさん
NGNGたたん。最低でもP-500以上。
0284
NGNGhttp://kaba.2ch.net/test/read.cgi?bbs=news&key=997090075&st=112&to=112&nofirst=true
0285デフォルトの名無しさん
NGNGhttp://mersenne.org/PrimeNet/
0286デフォルトの名無しさん
NGNG>もし素数の数が有限個だとして、それを小さい順に ｐ１、ｐ２、ｐ３、・・・、ｐｍとします。
>新たに、(ｐ１・ｐ２・ｐ３・・・ｐｍ＋１)なる数を考えます。
>この数は、任意の素数ｐｎで割れませんので、素数です。
（一部抜粋）
もし、素数の数が有限個だとして、それを小さい順に ２，３とし
新たに、２Ｘ３＋１なる数を考える。この数 ７ は素数。
もし、素数の数が有限個だとして、それを小さい順に ２，３、７とし
新たに、２Ｘ３Ｘ７＋１なる数を考える。この数 ４３ は素数。
もし、素数の数が有限個だとして、それを小さい順に ２，３、７、４３とし
新たに、２Ｘ３Ｘ７Ｘ４３＋１なる数を考える。この数は素数。
これをひたすら繰り返していけばいいのでは。
0287デフォルトの名無しさん
NGNGこれって… マジで
できそうでないかい？
0288デフォルトの名無しさん
NGNG0290デフォルトの名無しさん
NGNG> 新たに、２Ｘ３＋１なる数を考える。この数 ７ は素数。
> もし、素数の数が有限個だとして、それを小さい順に ２，３、７とし
> 新たに、２Ｘ３Ｘ７＋１なる数を考える。この数 ４３ は素数。
ここがおかしい？
素数を小さい順に 2,3,7 とするのはダメ、5 も素数だから。
小さい順に 2,3,5,7 としないといけない。
だから「この7は素数」と分かったときに
3から7までの数で素数を探す必要があるのでは
0291290
NGNG> 3から7までの数で素数を探す必要があるのでは
これが大変すぎるから、このスレがあるんだよね…
ageてゴメンヨ…
0292290
NGNG210
↓
211 は素数
(* 2 3 5 7 11 ... 211)
1645783550795210387735581011435590727981167322669649249414629852197255934130751870910
↓
1645783550795210387735581011435590727981167322669649249414629852197255934130751870911は素数
0293デフォルトの名無しさん
NGNG0294デフォルトの名無しさん
NGNG> 新たに、(ｐ１・ｐ２・ｐ３・・・ｐｍ＋１)なる数を考えます。
> この数は、任意の素数ｐｎで割れませんので、素数です。
そもそも↑これ、証明できるか？
その数が「Pm より大きくて (P1・P2・P3・・・Pm ＋ １) より小さい数」で
割り切れる可能性は否定できないと思うが。そしたら素数じゃない。
0295デフォルトの名無しさん
NGNGちょっとしたことをオンラインで検証できるツールがあるぞー。
■桁数無制限電卓
http://www8.big.or.jp/~000/CyberSyndrome/rsa/calculator.html
■巨大数の素数判定のアプレット
http://studwww.rug.ac.be/~bcallens/project/priem.html
0296295
NGNG0297デフォルトの名無しさん
NGNGそれにしても速い。 これがRabinの定理？
0298デフォルトの名無しさん
NGNGうむ。
１．ふるいで素数を小さい方から1000個生成
２．とりあえずその1000個で割ってみる
３．>>25で言う"任意のa"として小さい方の素数50個を使ってRabinテスト
…というアルゴリズムなように読める。
0299デフォルトの名無しさん
NGNG背理法って知ってますか。
これは素数の数が有限個だとして矛盾を導いているんです。
0300デフォルトの名無しさん
NGNG命題の対偶という概念が分からないと
背理法で何が導き出されるか分からないのでは。
0301294
NGNGんだから、素数が無限個であることの証明としては問題ないのはわかっちょる。
だからといって、素数を無限に生成する式としては使えるとは限らないだろ？
『 「素数の数が有限個」の仮定のもとでの非合成数 』 なら
確かに作れるけど、少なくともこのスレではそんなもんイラン。
0303デフォルトの名無しさん
NGNG0304デフォルトの名無しさん
NGNG0305デフォルトの名無しさん
NGNG２ヶ月もなにもないって事は、
終ったってことさ。
それにdat落ちしても、このスレのリンク覚えておけば参照できるよ。
0306デフォルトの名無しさん
NGNG過去ログについてもいろいろ 議論出てたからなあ・・・・・・
まあ落ちたら落ちたでしょうがないとは思うけど
0307デフォルトの名無しさん
NGNG過去ログ保存（バックアップ）しておくサイト
つくってくれないかなあ。全文検索できるようにしてさ。
0308デフォルトの名無しさん
NGNGこのまま落ちると 過去ログとして見られないのね
このスレにはFFTのソースがあるので勿体無い
0309デフォルトの名無しさん
NGNG調べたいウン万桁の数をメールで送る。
数秒後か数分後か数日後、結果が返ってくるのでは？
0310デフォルトの名無しさん
NGNG0311デフォルトの名無しさん
NGNGそれはwebサービスってことですね
0312
NGNGためしに作ってみたら
1から1千万までの数を判定するのに80秒かかった。
0313デフォルトの名無しさん
NGNGCPUは？
試しに作ってみたら、P3-850で14秒だったよ＞1千万まで
0314
NGNG速ければRSAとかさくっと。
0315312
NGNGP3-866 (^^;
少し改良したら1千万まで平均25秒にまでには短縮できたＹＯ！
0316312
NGNG0317デフォルトの名無しさん
NGNG素数判定が早くても RSA は破れないぞ。
必要なのは素因子分解
0318デフォルトの名無しさん
NGNG画面表示に時間がかかってるだけじゃない?
ファイルに出力するとセレ566で2秒くらいだけど…
0319
NGNGえー、一千万まで2秒ですかあ？
僕の奴だとファイル出力すると26秒だよ。
さらに画面表示まですると、ずっと遅くなる。
言語の問題とかも有るのかな。Javaで組んだんですけど。
0320デフォルトの名無しさん
NGNG0321319
NGNGhttp://pc.2ch.net/test/read.cgi/tech/1003569074/765-768
0322デフォルトの名無しさん
NGNGの方がいいんでないかい？
あと2以外の偶数は調べないとか
0323319
NGNG>limit = ルート(target)の方がいいんでないかい？
これで9秒削れて
>あと2以外の偶数は調べないとか
これで1秒削れました(^^
15秒程度になった！ありがとです。
0324デフォルトの名無しさん
NGNGint main(){
vector<int> v_primes;
v_primes.reserve(5000000);
v_primes.push_back(2);
for(int i=3;i<10000000;i+=2){
int limit = (int)(sqrt(i)+1);
for(vector<int>::iterator it=v_primes.begin();it!=v_primes.end();++it){
if(i%*it == 0)
break;
if(*it >= limit){
v_primes.push_back(i);
break;
}
}
}
return 0;
}
0325デフォルトの名無しさん
NGNG割算を使わないのでかなり速くなるはず。
メモリたくさん使うけど常駐するわけじゃないし…
# てか1000万って篩に使うメモリを意識して選んだのかと思ってた。
0326319
NGNGおお、劇的な変化が！
一千万までの演算が15秒から一気に1秒台へ^^;
>メモリたくさん使うけど常駐するわけじゃないし…
メモリ食いますね。
ていうか1億で計算しようとするとメモリが足りなくなってプログラムが
異常終了してしまふ。改良せねば。
0327デフォルトの名無しさん
NGNG0328319
NGNGboolean[] b = new boolean[100000000];
で落ちてしまふ。
この配列で、素数ならばtrue、合成数ならばfalseを入れるように
つくってるんですが。
0329Delフサギコ
NGNGミTДT彡
(∩∩)
最近のコンパイラでは
演算の高速化のために
booleanは実質32bit変数だから
TrueとFalseを判定するのに
すごく余分(32倍)にメモリを食っているのではないでしょうか。
0330仕様書無しさん
NGNGint[] b = new int[3125000];
とかやって、フラグを立てるときは、
b[i / 32] |= 1 << (i % 32);
みたいなカンジでＯＫと。
0332デフォルトの名無しさん
NGNG0333Delフサギコ
NGNGミ,,ﾟДﾟ彡 ＜ 一旦、sagaって終了したスレだし
,;ﾞ ﾐ ＼＿＿＿＿＿＿＿＿＿＿＿＿＿
ﾐ. ミ
／゛゛゛゛
Cマガで素数問題が出題されたから、復活ageされたのかな？
0334デフォルトの名無しさん
NGNG直接無関係なアルゴリズムで盛り上がってたからな。
つか「素数判定はエラトステネスの篩」っつー結論は>>4あたりで既に出てる
わけだ。名前間違ってるけど。
0336デフォルトの名無しさん
NGNG高速な掛算が必要なのは 数千、数万とかものすごい桁の素数を検定しようとしたから
0337デフォルトの名無しさん
NGNGんなこたわかってるって。だから>>334にもそう書いてあるだろ。128bitオーバーって。
0338336
NGNGFFT掛算が有効なのは数千桁以上だから 128bit見て反応しちまったよ
0339デフォルトの名無しさん
NGNGつーか、まず組込型で演算できない部分を計算するのにFFT使う、ってことで
さらに精度を保ちつつ係数を2^nにして高速化とかそんな感じで進んでたと思
うんだが。
数千桁以降じゃないと有効じゃない、と言われても、数千桁までは何でやるんじゃい、
というのは特に議論されていなかったはず。
0340デフォルトの名無しさん
NGNGlog10 2^(8+10+10+10)*4
=152log10 2桁
あってる？
0341デフォルトの名無しさん
NGNGたしか、686以降では36bitアドレスバスだから、
Windowsにこだわりが無ければ、64GBのメモリが使える筈だけど。
実際に対応OSがあるかどうかは知らない。
物理的に4GB以上メモリをのっけたとしても、ポインタは依然32bitだ。
0343336
NGNGFFTはあくまでも高速化のためで
数千桁以下は、 確か100進で計算させてたんじゃないかい
と思って読み返したら >>175-176で16ビットで計算してるようだよ
0344336
NGNG>>340
log(2)*8*2^32 = 10343311891.9
わずか10G桁=１００億桁
0345仕様書無しさん
NGNGということは、１つのプロセスが使用可能なメモリ空間は、
ＯＳが何であろうと４ＧＢまでってこと？
アセンブリ言語は１６ビット時代までしか知らないんだけど、
セグメントレジスタは３８６で拡張されてませんか？
0346デフォルトの名無しさん
NGNG一つのセグメントあたり4Gだから、
一つのプロセスが複数のセグメントを使えば
それだけメモリが使えるよね。
OSがそれに対応してるのかとか、
そんなコードをはくコンパイラがあるかのかは知らんが。
0347345
NGNG少なくとも、WindowsとLinuxは、４ＧＢが全実装メモリの上限だったハズ。
0348デフォルトの名無しさん
NGNGどっかの企業の改造版だったかも。
0349デフォルトの名無しさん
NGNG物理メモリ上限は、Win9x系は2G、NT系は通常2G、オプションで3G。
これを越えるのはまだ無いと思われ。
>>348
改造しても(規模にもよるが)4Gだと思われ。
それ以上となると、ちょっとやそっとの改造じゃ済まない（つか、１から作った方が速い）
0350347
NGNG１プロセスでもオーバー４ＧＢだったら凄いな。
0351デフォルトの名無しさん
NGNG物理メモリは64GBまでいけるね。
一つのプロセスでは4GBまでだとおもうけど。
0352340
NGNG>わずか10G桁=１００億桁
なるほど，ミスったよ
1G=1000MByte=1000*1000KByte=1000*1000*1000Byte=1000*1000*1000*8bit
=10^9
でもって
表現可能な数字＝2^(10^9)
ってわけか
で、こいつのlog10をとるのか・・・
そうすると，大体一桁落ちて
数千万桁
あれ，どっかで計算ミスってる？
まあ高速なアルゴリズムよりも
高速な分散処理を考えないと普通に無理だよね
分散って言っても，一つの数値の判定を分散処理する
数学板から誰か連れてきてよ
でもって作る？
通信するので、言語は何でもありですね。
まあ、C++ですかね。
0353デフォルトの名無しさん
NGNG手元にあるWindows2000 RC2の説明によると、ProfessionalとServerは4GB、
Advanced Serverは8GBまでサポートしてるよ。
素数関係ないのでsage
0354仕様書無しさん
NGNGでもな〜、素数の判定にロハで協力すれ、って言って、人が集まるのかな。
SETI と比べてジミかつマニアック過ぎると思われ。
100万桁の自然数ひとつ判定するのに、どれくらいの台数のマシンが必要
なんだろか。
一つの数の判定するのにかけていい時間を1日くらいと考えると･･････
数十万台はいるんじゃなかろうか。
0355デフォルトの名無しさん
NGNGかっこいいスクリーンセーバー作って、裏で勝手に計算しちゃえよ（笑
0358デフォルトの名無しさん
NGNG>「素数判定はエラトステネスの篩」っつー結論
デマを撒くな〜
>>25
の確率性アルゴリズムを使うか、
特定の構造をもった数に対する検証のどちらかしかないぞ。
スレのテーマが「最速で素数表を作る」ってなら櫛で間違いないけど。
0359336
NGNG4Gバイトなら10G桁表現出来るって事ね
0360デフォルトの名無しさん
NGNGインターネット上のPCの共同作業で。
http://mersenne.org/status.htm
0361デフォルトの名無しさん
NGNG遠いネットワークの先にあるメモリーよりも先に地元のディスク使わないの？
0362デフォルトの名無しさん
NGNG10億までなら200秒かなとか単純に思ったら、すんげえ時間かかる。
ハードディスクがガリガリ言ってるし。
10億でもたったの10桁・・・素数計算って大変だねー。
0363デフォルトの名無しさん
NGNGうわ。見逃していました。すげーニュースだ。
まだ検証中（ダブルチェック中）なんだね。
M38 が (2^6972593) - 1 で M39 が 8000000 以上ってことだから
かなりの値だ。
0364デフォルトの名無しさん
NGNGなに、Ladinの定理?
0365デフォルトの名無しさん
NGNG0366デフォルトの名無しさん
NGNG0367デフォルトの名無しさん
NGNG１Gのデータをハードディスクから読み出すのにどれくらいかかると思ってるん
だか・・・・
0368デフォルトの名無しさん
NGNG1GBの読み込み1.3秒ですが何か？
0369デフォルトの名無しさん
NGNG(64bit,66MHzで533MB/sec)
0370デフォルトの名無しさん
NGNG0371デフォルトの名無しさん
NGNG0372
NGNGhttp://news.2ch.net/test/read.cgi/newsplus/1007792330/l50
0374デフォルトの名無しさん
NGNG0376デフォルトの名無しさん
NGNGこっちのほうが見やすいと思われ。
http://slashdot.jp/article.pl?sid=01/12/20/155216&mode=nested
0377デフォルトの名無しさん
NGNGC言語による最新アルゴリズム辞典の改造しただけだけどどうでしょう？*/
#define N 4000
int prime(unsigned int p)
{
unsigned int h, i, j, k, s,cash,cash2;
char a[N+1],x[N];
memset(a,0,p);
a[2] = 1; /* a = 4 */
for (k = 2; k < p; k++)
{
memcpy(x,a,p);
memset(a,1,p);
a[1] = 0; /* a = -2 mod M_p */
for (i = 0; i < p; i++)
{
if (x[i])
{
s = 0;
cash=p-(h=i)-1;
for (j = 0; j < cash; j++)
{
s = (s >> 1) + a[h] + x[j];
a[h++] = s & 1;
}
{
s = (s >> 1) + a[h] + x[j];
a[h++] = s & 1;
h -= p;
}
for (j++; j < p; j++)
{
s = (s >> 1) + a[h] + x[j];
a[h++] = s & 1;
}
if (s > 1)
{
cash2=h;
for(j=0;a[h]&&j<cash;j++) h++;
if(a[h])
{
h++;
memset(a+cash2,0,h-cash2);
cash2=0;
h -= p;
while(a[h]) h++;
}
memset(a+cash2,0,h-cash2);
a[h] = 1;
}
}
}
}
a[p] = 1 - a[0]; /* 番人 */
i = 1;
while (a[i] == a[0]) i++;
return (i == p);
}
0378デフォルトの名無しさん
NGNGとフト思ってしまった
0379素数が理解出来ない技術者
NGNG1億までエクセルのＶＢＡで計算しましたが、実行時間は、約２０時間、反例は
見つかりませんでした。
これって、結構早いのでしょうか？
ちなみにマシンは、ＶＡＩＯのPCG-FX55/XPで、メモリーは６４ＭＢ。
なお、ゴールドバッハの予想を知らない方への説明。
全ての偶数は、２個の素数の和で表される。
そこで、ある制限をかけて、いかなる偶数の平方根の２倍＋１よりも大きい
素数と、その他の素数の和となることは無いと予想したのですが、果たして反
例は見つかったか？結果は？
興味のある方は、掲示板に書きこんで頂けば、マクロが入ったエクセルのソ
フトの配置してあるＵＲＬを書き込みます。
なお、今のところ、私のエラトステネスの篩は、１００万の計算に１２秒かか
りますので、まだ改善の余地有り。
なお、定義に従えば、１００万の計算に８日かかります。
0380デフォルトの名無しさん
NGNGある程度ていうのも、へんな言い方だけど…
0381デフォルトの名無しさん
NGNGVBで実行してみたいです。
0382デフォルトの名無しさん
NGNGlog(2^n-1)はおよそlog(2)*Nだから 2^n-1〜2^(n+1)の間にも0.7*n個程度あるという事
0384デフォルトの名無しさん
NGNG○ だから 2^n-1〜2^(n+1)の間にも(2^n)/(0.7*n)個程度あるという事
0385デフォルトの名無しさん
NGNG0386デフォルトの名無しさん
NGNG0387デフォルトの名無しさん
NGNGここ参考にするといいでしょ
http://momonga.t.u-tokyo.ac.jp/~ooura/fftman/
特にこの絵
http://momonga.t.u-tokyo.ac.jp/~ooura/fftman/fig_bf2.gif
最後にビット反転置換をしなくても良い方法もあるけど、
すると毎回複製が必要になるのと、ループが少し複雑になるから
メモリの節約になるこの方法が好まれているのでしょう
0388デフォルトの名無しさん
NGNG0389デフォルトの名無しさん
NGNG確定的判定法で最速の物ってどれよ？
ecppだっけ？あれはどれぐらいのオーダー？
0390デフォルトの名無しさん
NGNG0391デフォルトの名無しさん
NGNG0392デフォルトの名無しさん
NGNG0393デフォルトの名無しさん
NGNG0394デフォルトの名無しさん
NGNG例えば1万個の素数を用意しといて割ってみて割り切れなかったら 素数でない確率は1万分とか
0395デフォルトの名無しさん
NGNGオオバカ
0396デフォルトの名無しさん
NGNGRabinの方法は、対数和にしかならない
0397デフォルトの名無しさん
NGNG0398デフォルトの名無しさん
NGNG0399デフォルトの名無しさん
NGNGテーブルを分解・圧縮して計算に使う部分をなるべくキャッシュに載るようにして
レジスタも使って（速度の桁が違うので）
2の倍数乗除にシフト、（プロセッサ内で勝手にシフトに切り替わるかな？^^;）
遅い割り算を逆数掛け算にしたりとかは効果ないんでしょうか？
0400デフォルトの名無しさん
NGNG大規模３Ｄレンダリングなどで使う
クラスタリングが分散にいいらしいです
0401デフォルトの名無しさん
NGNG並列計算も、完全に成功しても台数分 数が増えるだけ
それは頑張って１０台並列でももたった１桁増えるかどうかのレベル
掛算をFFTにしたような、一挙に桁数を倍に出来るようなアイデアが求められている
0402デフォルトの名無しさん
NGNG■ このスレッドは過去ログ倉庫に格納されています