最速の素数判定アルゴリズム　　

最速のアルゴリズムはRabinの定理を使ったものと聞き、
文献を読んでみたのですが、サッパリわかりません。
日本で完璧に理解している方、たまたまこのスレ
見ていたら、ついでに教えてください。
マジです。

>>1
理解できないのならあきらめるしかないでしょう。
他人をたよるのはよくありません。

>>1
その「さっぱり」って言葉がやるきねぇなぁと他人に感じさせるので、
そういう聞き方では知識ナルシスト以外教えてくれないと思われ。

エラストテネスの篩で我慢しろ。

素数なんて、計算するものじゃなくて
データとして保持しておくものだよ。
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

>>5
神父もびっくりだな

下のほうは興味あるなぁ。
日本のチームはあるの？
京都産業大学なんか大学でチームつくってそうだな。

100万桁の
100万って、どのくらいなんだ？

100万文字ってなんKB？

百万桁と百万文字は違うだろ

>>8
BCDで持つとして、100万桁=50万バイト=5000キロバイト=5メガバイト

量子コンピュータを作って計算させる

http://piza.2ch.net/test/read.cgi?bbs=tech&key=991992977&ls=50

素数ってなァに？

>>13
整数環の素元のこと

いや、Delで文字列を数値演算として
扱えるユニットがあるので
それで扱える桁数かな、とオモタわけ。

5MB程度ならいけるな、100万桁．．．．

ところで素数って正の数かしら？

素数判定が速いソフトってないかな。

だから、、、素数のリストを内部的に保持しておいて
そのリストと照らし合わせる判定方法にすれば
それこそ速い遅いの問題ではないっていってるんだが？

通じないか？

＞17
リストってお前、500bitとか1000bitとかの素数まで全部保持するつもりか？

>>18
>32bit整数MAXまでの素数リストがあれば実用上何も問題なし。
って言ってんだろ?

>>11
HNにちょっとだけワラタ。　(^◇^;

>>18 の頭は、>>5 に影響されて、どっかに飛んでいきました。

技術的に考えたいのなら
メモリ管理とか考える方が先。

エラストテネスで32bitMAXまでの素数リストを作ろうと思っても
メモリがオーバー風呂して出来なかったYO。

あと、何万桁を計算する手法も考えるべき。
32bitとかそういうレベルじゃなく

それでいて素数判定が早いというのには意味があるが
実用的かどうかというと全然実用的じゃないと思う。

学術的には意味はありそうだな。

そもそも、1は
Rabinの定理というものが知りたかったのか...

で、その定理の詳細ページどこかにない？

http://www.hotwired.co.jp/news/news/technology/story/2234.html

これって適当に答え提出して「これは素数だ」って言い張ったらどうなるのだろう？

＞２３

http://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 が存在する．

これを使って　適当に乱数で検査して　素数らしさが十分だったら、それで
いいってのが１が欲しい乱数の作成方法でしょ？

まず a ≡ b mod n は C言語風に書けば (a % n == b % n) == True の事

素数が注目されているのは暗号化技術に必要だから、
例えば、公開鍵暗号には２つの素数 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がホントに素数かどうかより　素数らしい値を手軽に
得られる方法でいいという理屈。

× e*d=n;　n % (p-1)*(q-1) == 1　になるように e,dを見つけて
○ n=p*q; (e*d) % (p-1)*(q-1) == 1

だめだ。
Max32ビットの素数リストすら作れない、
メモリの確保が不可能だ。

37万くらいまでなら泥臭い方法で出来たけど。
10桁ですら、とてもとても。

>>1 の知りたいこととは全く無関係に、やや面白いスレになりつつありますな

9999929/9999931/9999937/9999943/9999971/9999973/9999991/

このあたりまでのリストなら
ならなんとか....

>>28,30さん
エラストネスのふるいを使えばどうでしょうか？
基本的にある数Nが素数かどうかはN^1/2までの約数が無いかを調べればよいので、
N^1/2を超える最小の1からのすべての素数の積(変な日本語ですいません)
までのパターンがあれば、後はNをその積で割った余りを検証すればよいのでは？
1*2*3*5*7*11*13*・・・・・・って計算していって0xFFFFを超えた数を見つければ
良いと思います。

>エラストネス
エラストテネス？恥ずかしい・・・

て全然駄目ですねごめんなさい。あげないで

>>26
中途半端な知識だな。
それはRSAな。
○RSA∈公開鍵暗号
×RSA＝公開鍵暗号
×RSA∋公開鍵暗号

RSAは会社名だよ
何の略か忘れたけど・・・

>>35 (54??)
RSA暗号は発案者三人の頭文字をとっている。
Rivest
Shamir
Adleman
で、こいつらが作った会社の名前もRSAセキュリティ

30だけど
>>31
つかってるよ。

ただ、エラスネスするためには
配列やリストを用意しなければいけないんだけど
Booleanの配列で
32bit個の入れ物をつくっちまうと
メモリオバなのよ。

直接メモリのビットを見ればいいんだけどねえ。
軟弱プログラマだからやり方しらないよ。

何だバカヤロー
エラトステネスだバカヤロー

４２億まで計算しようと思ったが１億まで計算するのに４０分
かかった（PenMMX266）ので途中で断念した。時間かけ
ればできるけど、それでもまだ１０桁だよ。１００万桁なんて
想像できん。

２ｃｈで共同で計算するのもいいかもね。

任意の数を与えると、それが素数かどうかを判断する関数を用意すれば？

>>41
どういうアルゴリズムを使うわけ？
素数の定義からは逃げられないよ。

ルカステストとかラビンの定理とか。

>>39
激しく同意。

100桁の数一つ扱うのにも
めちゃくちゃややこしいし
一つの数が素数か否かを判断するのにも
膨大な時間がかかる。
この部分を分散計算するという考えそのものが思いつきもしない。
ｳﾂﾆｬ

100万桁くらいの+ - * なら ややこしくないだろう。 / % は多少メンドクサイけどね。
　BCDや１００進法で 500KByte 掛算時のΣai*biで桁溢れがおきないから大丈夫

でも、素数かどうか判定するのはどうかな？
１の言うRabinの定理で確実に素数らしい奴の予想をつけといて、それを皆で割って
みるしかないとすると、　参加者が何人いてもあんまり助けにならないよね

http://www.hotwired.co.jp/news/news/technology/story/2234.html
>素数とは、1、3、5、7、11、13など、
>1およびその数自体でしか割り切れない数

大学４年生にもなってはじめてしりましたよ。

7時間かかってらあ。
....99999677/99999703/99999721/99999773/99999787/99999821/99999827/99999839/99999847/99999931/99999941/99999959/99999971/99999989/

>>47
重箱隅だが2も素数だよ。

ギコbasicが多桁対応になってるから、これでなんとか出来ないかい？

やるならDelで文字列数値ユニットがある。

で、メモリ直読みでエストラネスをやってみようか。
符号なし32bitMAX程度は出したいですね。
4294967295

でも
符号なし64bitMaxの18446744073709551616
は出来そうにない...
これでも100ﾏｿ桁まであと....

気が付いたけど、Windows付属電卓って
すごい桁数扱えるのね。

妙な所を
きっちりと作ってやがんな＞＞MS

1は素数じゃないのでは・・・。

1は素数じゃないな。。。

>>4>>22>>32　エラストテネス
>>31　エラストネス
>>37　エラスネス
>>51　エストラネス

荒井注の冥福を祈りつつsage

1は素数だろ。
>>53-54どこみていっているんだ？

ウソだった。氏ぬ。

素数定理からすると、
10^(10^6)〜(10^(10^6) + 2^24) の範囲には素数がありそうだと思える
これらの全部のフラグをたてといて、
Rabinのアルゴリズムで最後迄残った数個に絞り込むのはそれほど大変じゃない筈

後は、これが素数かどうかを証明する方法だなあ・・・

>>58
１と自分以外で割り切れるかどうか、試せばいいんじゃないの？

>59
いや、あの、それじゃ１命令10GHzで多桁除算が出来たたとしてもさ
　10^(10^6)/10^10 = 10^(10^6-10)秒かかるんだよ

10^(10^6) のオーダだとフルイに使うメモリにしたって
世界中のパソコンのハードディスク足したって足りはしないさ。

正確に書くと 2^24 の範囲には素数は 平均7.3個程度ある筈

それから　√n迄でいいから、上の計算よりずっと少なくて
　 √(10^(10^6) )/10^10 = 10^(1000-10)=10^990秒でいい

100万=1M桁になると、掛算を筆算方式でしたんじゃ 1回の掛算だけで１時間かかるんじゃないの？

FFT方式だとどれくらい早くなるのかな？

>>60
ハゲしく同意

フルイにかけるメモリ確保すら問題

Shorのアルゴリズムで素因数分解できるかどうかみてみるってのは？
だれか量子コンピュータを作ってよ。

a=を大きな数として b=log(a)程度の数で

for (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万分の１になるんじゃないかな？

1000000!がいくらになるか想像つきますか？

ちなみに、スターリングの公式だと、
n! = n^n * e^(-n) * SQRT(2*pi*n) e^(1/(12*n) - 1 / (360 * n^3) + θ/(1260 * n^5))
(0<θ<1)
です。

その前によく考えたら、100万迄は a%cはゼロだけど 1000000! % 1000001 はもう予測出来ないな
100万分の１じゃなくて100万分の節約ってだけだな

Log[10.0,Gamma[1000000.0]] ≒ 5.57 10^6
Log[2.0,Gamma[1000000.0]] ≒ 1.85 10^7

1000000!は10進数で約557万桁、約1850万ビットの数ですか....
想像つきませんね。

でも課題は、
　100万桁よりも長い最初の素数　5万ドル
　10億桁よりも長い最初の素数に対する最高金額の25万ドル
だからなあ

１００万桁はもう出てると思うよ
http://www.eff.org/awards/coop.html

狙うならもう１桁上だね

1000001桁？

求めるならメルセンヌ素数を求めたいね。

>>72
>>71のリンク先を見る限り
>>71の言う一桁上は
1,000,001でなく
10,000,000ってことかと。

>>73
なんで桁数が10進法で上昇しなければならないんだ？
意味わからん。

>>74
俺に言われても(^^;
>>71の言い方が悪い…。

レス遅れてすいません。
僕が参考にしていた文献です。
これをもとにして、高速な素数判定ソフトを作りたかったんですが、
無学なもので、内容が理解できませんでした。
そういう悲しい事情です。

【参考文献】
「アルゴリズムの基礎」
五十嵐善英、西谷泰昭　共著
コロナ社
ISBN4-339-02348-5 \3,000

≪8.5　Rabinの素数判定アルゴリズム（P168−170）≫

普通のプログラマとしての意見だけど
32bitのintの最大値
2147483647
は素数なんだけど、これを泥臭く2.3.4....2147483646でわってみる判定方法で
やったとしても今のマシンなら数分(今試したよ)なのだから
高速素数判定ソフトを作っても
何の意味もないような気がする。

32bitや64bitを超えた数値を扱う方法は
アルゴリズム以外の面で実装が大変だ

>>1さんは
その辺りの事、わかってらっしゃいますか？


2147483647までの素数列を全て求めるときの
高速アルゴがあるのなら教えてもらいたい。
アルゴリズムというより
そのデータを保持するメモリ管理の方が知りたいな。

　その高速判定は、あくまでも確率的な判定方法ですよ
　ソースはこれから辿って下さい
http://search.net-newbie.com/php/function.gmp-prob-prime.html

>>77
　実装が大変なんですか。あまりよくわかっていませんが、
桁が多いと、扱える変数（？）がなくなると言うことですか？
それにともなって、細工が必要になるのかな。素人でスマン。

＞2147483647までの素数列を全て求めるときの
＞高速アルゴがあるのなら教えてもらいたい。

僕が最初、試したアルゴリズムは、
６ｎ±１（ｎ≧１、ｎは自然数）を、素数の候補として、
手当たり次第に√（６ｎ±１）に達するまで、

（６ｎ±１）ｍｏｄ（既に得られてる素数）＝０

を判定する方法をとったんですが。これってダメ？
でも、これだと、2147483647までできないと思う。
やってみてもいないし。メモリ管理に至っては全然です…

2^23程度なら256M+128Mバイトのメモリを実装すればオンメモリで可能だろう

×　2^23
○　2^31

それから2^31迄の全ての素数を入れるのに必要なメモリは
およそ 4*( 2^31/ ln(2^31)) だいたい 400Mバイトだから　フラグで記録した方が効率がいい

>>77
>>78

＞やったとしても今のマシンなら数分(今試したよ)なのだから
＞高速素数判定ソフトを作っても
＞何の意味もないような気がする。

＞その高速判定は、あくまでも確率的な判定方法ですよ

文献によると、
≪Rabinらはこの方法により、２＾ｐ−１なる形の整数が素数であるかどうかの判定を計算機を用いて行った。1<p<500の全ての整数についての判定をたった数分で行い、全ての判定は正しかったと報じている。この方法は素数であるかどうかを判定するどの決定性アルゴリズムより桁違いに速く、サンプル数ｍを適当な大きさにすることにより、誤りの確率をきわめて低くすることができる。≫

≪例えば、ｍ＝５０のとき、誤り確率は約千兆分の１となり、基本演算の所要時間が１μsで、３万時間連続稼働している間に１度の割にしか基本演算の誤動作が生じない計算機よりも信頼性が高い判定方法である。≫
とあります。

これ見て、作ってみたいなと思ったわけです。

>>80、81
アドバイスありがとう。そういう計算方法も知らなかったんです。理系だけど。ちょっとやってみます。

フラグで記録する時は 偶数を記録しないだけで 1/2に　さらに３の倍数を記録しないのを追加すれば 2/6
でいいからね

>>82　だから >>25 が書いてるだろ？
＞　　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という関数が判らないのか？

/* 符号なし 32bit の素数を全て求める(返り値:32bitの素数列へのポインタ(先頭は素数の数)) */
U32* 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()とかしてるし)。
あと、メモリ確保できなかった時の対処もしてない。

とりあえず、答えが出るのが何時になるかは知らないが
出せないなんて事は無と思う。

さらに噛み砕いた説明が欲しいなら

http://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 とユークリッド で検索してみろ

>>85
そうなんです。わかりません。gcdという関数の存在も今知りましたが、
それでも、文章の意味が分からないんです。
しかも、どこにツッコミ入れられているかもわからないんです。
つきはなさないで下さいね。

>>86、87
お忙しい中、ありがとうございます。
力不足ですいません。
時間かかりますけど、やってみます。

>>64
Shorのアルゴリズムも、Rabinと同程度に確率的な判定法でしかないっす。

ルカステストは？

>>90
Shorのアルゴリズムとは

「コンピュータと素因子分解」
和田秀夫 / 遊星社

最新の資料ではないがよくまとまっているし、プログラムコードも入っている。
専門家でなければ必読。

やっぱアルゴリズムとかコアな本は外人が書いたやつに限るね。
日本の東北大の教授とかが書いたやつなんてコードに無駄なとこ多くて読む気がうせる

>>50
これだね
http://www.geocities.co.jp/SiliconValley-Oakland/8522/data/mformula.txt

モジュラベキとかgcdとか関数を増やしてけば使えるかも・・・でも
既に1000桁扱える１０進basicとかあるからなあ

アドバンテージはメモリがある限り桁数増やせる点だけど、今のアルゴリズム
じゃ掛算が遅くてダメだな

>>92
おおざっぱに書くと

１．Ｎを素因数分解したいとする。
----- 因数が見つかるまで２〜５を繰り返す ------
２．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程度。
-------------------------------------------

ステップ３．を計算するのは量子コンピュータじゃないと死ねるので結局まだ実現不可だけど。

>95　今のままだと
　掛算で 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;

複雑な事はわからないですが

2.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

とやっていくと
計算量が減り、まあまあ高速に素数リストが出来て
メモリのことも考慮しなくていいようになるかも。

ダメ？

ガイシュツだっけ？

ここは！！！
http://www.ie.u-ryukyu.ac.jp/ie/mirror/www.iip.ie.u-ryukyu.ac.jp/yuu/c_html/sosuu.html

>>48の方がよっぽど上だった,,,

>>98
それより計算量が多かったり低速だったりする方法って思いつかないぞ。(w

√N を計算するより N/2 以下の素数全部で割ってみる方がマシってこと？

>>101
いや、よくわかってないもので。

>√N を計算するより N/2 以下の素数全部で割ってみる方がマシってこと？

詳しく教えて。
ルートエヌ？？

エラトステネスより98の方が現実的だと思ったのさ。

http://www.valley.ne.jp/~kanai/IT/sources/java/primeTable.html
Javaなのに速いな....

>>31に書いてあった。
＞基本的にある数Nが素数かどうかはN^1/2までの約数が無いかを調べればよいので
欝...

計算量の比較

これまでに確認した全ての素数で割る方法
　N以下の素数の数はおよそ N/ln(N)
　M迄求める計算量は ΣN/lnN 回の割り算
　ただし、素数でなければ途中で割り切れるので、殆ど負荷がないとすれば
ΣN/ln(N) /ln(M)　回の割り算

いっぽう篩方式は およそ M/ln(M)*Σ(M/N))回の掛算
　殆ど同じオーダだけど、片方は割り算、片方は掛算だから圧倒的
フルイを使う方が早い

√N が出て来る理由は

判定したい値を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 の点

>>105
スマSO、なぜかうちでは
フルイの方が圧倒的に遅い

これまでに確認した全ての素数で割る方法 でおよそ2分
フルイで10000000までの結果(>>30)を出したときがおよそ25分

どういうことじゃ？
そのうちソースだすよ（Delなんだけど
どこかにくだらないオーバーヘッドがあるのかもしれない。
出力時に全ての数値をなめなければ結果を出力できんし。
(わかりやすいように計算と結果は別々にしているんだ。)

メモリ効率の点からも
全ての素数で割る方法じゃないと
64(or63)bitMAXまでは求めるのが無理らしいです。

2^63ビット確保できないってさ。


>>106
ありがとう。めちゃためになる。
数学科出身？？

>>105
ふるいに使うのは足し算だけで十分だと思われ。
今とりあえず組んでみたら1億までで30秒くらい。

>>108
http://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ってのが入っているんだけど
このやり方がダメなの？

>>このやり方がダメなの？
ダメダメダメダメダメ。でもうちの新人に[i...j]のうちのnの倍数を
全て出力するプログラム、って課題を出したら、最初は確かに mod (%)
使ってた...

じゃ、教えてくださーい
ある数がある数の倍数かどうかを判断するとこがダメなのか？

ある数の倍数かどうか

たとえば１０進数なら

３の倍数かどうかは全部の桁を足して３の倍数か見るだけでいい
５の倍数かどうかは最下位が０か５か
　っていうのがあるよね

２進数の場合はどうだろう？
３の倍数　２桁毎に足した結果の下２桁が11である
５の倍数　４桁毎に足した結果の下４桁が1001

という事で、Nで割れるかどうかは
(2^n-1)/N で割り切れるnを探して n桁毎に足した結果がNなら
というのでいいのかな？

あ、フルイの場合は、単に Nの倍数は Nを足すだけでいいから掛算も不要という意味ね

でも、>>112の疑問もお願い

あ、間違いだ
　n桁毎に足して、その結果がさらにNの倍数かどうか見なくちゃいけない

３の場合 3*17=$33 ２進で２桁毎に足したら６だから 0110ｂ さらに２進２桁毎に足したら 11ｂ

子ね

すんまそ。
ほんとにバカなんで教えてくださいませ

倍数判断はともかくとして
>>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
この部分が怪しいような気がする。

112みたいな方法で2進数で3で割ったあまりを判別する方法を考えると、
各桁ごとを+-反転させながら足した総和が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あまる。

多分実用性はないな。

> -1+0-0-1+0=0が3で割り切れるので3の倍数。
訂正 +1-0+0-1+0=0 だ。

とりあえず 素数Nが割り切れる 2^n-1 を求めてみた　29は判らん
　 3　 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

まあ、２進法で８ビット単位に足して　全部の桁の総和が
３，５、１7 の倍数かどうかチェックすれば一挙に３つ計算出来る
から、多少は早くなるんじゃないかな

112がやろうとしていることが今分かった。
pが素数なら、2^(p-1)-1 が p で割り切れるよ。

暗号化とかで計算に必要な素数ってのは、だいたい
どのくらいの桁数なの？

>>121
あ、ホントだ　という事は 29は 28bitで循環か、なんで抜けたんだろ　浮動小数点使ったせいか

>>122 必要なだけ　100桁くらいあれば個人ならなんとか
　千桁くらいあればとりあえず

　　　　For I = 2 To Limit
　　　　　　If sosuu(I) = 1 Then
　　　　　　　J = I * 2
　　　　　　　Do While J <= Num
　　　　　　　　sosuu(J) = 0
　　　　　　　　J = J + I
　　　　　　　Loop
　　　　　　End If
　　　　Next I

なんかの倍数なんてチェックする必要すらないような気がするが...

１０進法で記録した場合
　　2　　最下位桁のみ
　　3　　1桁単位に加算
　　5　　最下位桁のみ
　　7　　6桁単位に加算
　 11　　2桁単位に加算
　 13　　6桁単位に加算
　 37　　3桁単位に加算
　 41　　5桁単位に加算
　 73　　8桁単位に加算
　101　　4桁単位に加算
　137　　8桁単位に加算

>>124
J=I*2　にしてるけど、Jより小さい素数の倍数は不要でしょう？
説明すると、
まずは２刻みに４以上を全部消すよね？
次は３刻みなんだけど６は既に消えてるし２の倍数は消えてるから９から消す
その次は５刻みだけど
やっぱり５＊２　も５＊３も消えてるから５＊５から１０刻みに消せばいい

つまり素数N番目は N*Nから 2N刻みに消せばいい

Stepが使えるかどうか知らないんだけど　StepN はN刻みって事ね
で、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　がループ中再評価されるんなら外に出してね

メモリーが２５６メガくらいあれば動くみたい
MAXNに設定出来る最大は 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;
}

最後の何個か
2147483477
2147483489
2147483497
2147483543
2147483549
2147483563
2147483579
2147483587
2147483629
2147483647

途中まで遅いけど、１０００を越えた付近からスクロール速度の方が先に限界になります。
ハードディスクにリダイレクトさせると１ギガを超えそうになったので、プログラムを弄って後ろの方だけ印刷させました

>107のいう 10000000 迄ならハードディスクに書き出しをすれば３秒で終わった
　画面に出したらスクロールに３分かかった

>12

そのスレをたてた「悲惨な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

＞　１９７１年、旧ソ連のマチアセビッチは、素数全体をあるひとつのディオファントス
＞　方程式の解として表すことにも成功しています。すなわち、すべての素数をつくる式
＞　を生み出したことになるのですが、その式は３７次で２４個の変数をもつ多項式と
＞　２１次で２１個の変数をもつ多項式でした。

興味があり調べているのですが、なかなかこれに関連する文献がでてきません。なにかご存知の方はいませんか。

http://www.fortunecity.com/emachines/e11/86/primesec.html
???

http://photo-m.tp.chiba-u.ac.jp/~adeno/sci/math.htm

＞それは１９変数の多項式であり、１９個の自然数をその
　多項式ｆに代入してｆの値が正であればそれが素数であり、
　逆にすべての素数についてもｐ＝ｆ（x_1, x_2, ..., x_19）を満たす自然数 x_1, x_2, ..., x_19が存在する。

という事だから、あれば便利そうだけど、２０次とかの式を含んでるのなら例えば500bit前後に収まる
ものをというような要求には難しいんじゃないでしょうか？
オイラーの式 n2+n+41 の方がコントロールし易いんじゃないかと思います

素数の判定をそんなに急いでどうすんの？

>>91 ルカステストってメルセンヌ素数かどうかのテストじゃなかったっけ？
１００万桁の最初の素数を見つけるような場合には意味が無いとおもわれ

全てが素数になる方程式なら ウィルソンの定理があるんじゃないの？
（N−１）！＋１　≡ 0 mod N

if ((（N−１）！＋１ )% N==0) 素数です　else 素数ではありません

>>133のいうオイラーの式って、その数列の上に全ての素数が
存在する訳じゃないし、それを満たせば全部素数って事もないよねぇ。
ていうか当然か。。。

>>133の言ってる素数を生成する多項式っていうのは聞いたことがあるけれど、
代入しても負の値が出ることがかなり多いそうだよ。(つまり素数じゃないやつネ)

>>130
話それるけど、
そのようなプログラムでコンソールに連続的に表示するとき、
Ctrlを押してるとスクロール（表示速度）が速くなるんだけど
理由知っている人いる？

>>128
for 文で 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)

とオーバフローしたら検出して抜ければいい

>>68
＞1000000! % 1000001 はもう予測出来ないな
　１３６が書いているウイルソンの定理から 1000001 が素数なら
　　1000000! % 1000001 =1000000 です
　しかし、1000001=101*9901なので階乗は割切れますから結果はゼロです
1000003が素数なので
1000002! % 1000003 =1000002
1000001! % 1000003 =1
1000000! % 1000003 =500001
と予測出来ます

>>126-127=112さんありがとう。
ついていくのに必至だ。

>>48で7時間って言ってたけど
それが5分になった（ｗ
うれしい！

手元のマシンにメモリが少ないのが悲しい所。

Delだけど
http://www.mogya.com/develop/kanjicalc/develop.html
10億桁扱えるユニットがあるので
メモリを食わない＆計算時間が1日程度の方法で求められるのなら
適当なルーチン(わかりやすくたのんます)書いてくれたら試すよ。

>>143　そのルーチンじゃ、掛算は桁数の２乗に比例するでしょ？
1Gマシンで100万桁の掛算だけで既に１日で出来ないと思うよ

>>144
そっか。やはりWiredNEWSで言われているように
分散処理で考えないとだめか。
100万桁といわず、ここでどんどん桁数を上げていきません？

いや、もっとずっと小さい桁数でダメになりそうだし
そういうふうに動的に取るとメモリが細切れになって大きな配列が最後は取れなくなる


C=A*Bという計算をする時に Setlength(C,length(A)+Length(B))とメモリを確保してから
余ったら開放する方式の方がいい

で、かける方の結果の桁でループを回すんじゃなくて、>>97のように
かけた結果の桁でループを回すといいよ。その方が1桁以上早くなる筈

ベースにするなら　>>95　のギコbasicベースの方がいいんじゃないの？
使い方は文字列で与えて文字列で出て来るから殆ど同じだし
入力で "2^1000+1" みたいな書き方も出来るし

内部１００進だから１０進に比べてそれだけで掛算が４倍早くなるし

どっちにしても、掛算について万桁超えたらfft使わないと苦しいと思う

解説はここ
ttp://homepage2.nifty.com/m_kamada/docs/fftmul.htm

そっか。
unit mformula;
も数値を文字列で持っているんだね。

でも、なんとなく使い方が難しそうなのだが...
Stringで数値を保持しておいて
計算するときはFormulaNumStrして
結果もStringに保持しておけばいいって事なのかな？

それと、目標を作らないとチューニング出来ないと思うよ
　例えば1万桁台の最初の素数を探してみたいとか

ルカステストで100万桁台のメルセンヌ素数を調べたいとか
やっぱり特定の目標にしないと、　桁溢れとか色んな問題がクリア出来ないと思う

試したときは
function N100Str(str:String):String;
begin
　 Result:=FormulaNumStr(TN100Num,str');
end;

を自前で作って　 Caption:=N100Str('100^10-1');
みたいに試したからそんなに難しいとは思わなかったけどね
　無理に関数使おうとすると、クラス貰って来てどうとかやっかいだね
どうせ長桁になったらインタプリタの時間なんてしれてるから インタプリタに任せた方がいいんじゃない？

小数点の使える TN100の他に 整数だけのs100 もあるみたい
あんまり試してないけど、素数見つけるならこっちでいいんじゃない？

あれ？さがってる。

今日こそは金曜日なので
月曜日まで計算させて
32bitMax値をたたき出そうと
思っていたんだが。

63bitMaxまでいっけるっかな。

最初に、"素数は決めうちでリストでもっておいて"
それにヒットするか否かで判断しろっていったけど

これだけ莫大な量だと
その方法が正しいのかすら怪しいな。

６４bitまで計算させたい時のことを考えたてるんだけど
a,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) % と分解出来るかな？

mformulaを書いたものです。
　色々コメントありがとうございます。
　掛算の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;

こうアセンブラで書く事が出来ます。　

>>153
「掛算のFFT化」とは？なんですか？
FFT？

がんばってください>>153

私はとりあえず、泥臭い方法(>>98)で
普通の整数を使えば32bitまでは
出せそうなのでそれでやってみます。
土日計算させてたら、答え出るんじゃないかと期待

64bitは計算時間的に不可能だろうな。

>>154
FFT=高速フーリエ変換

ありがとう。
じゃ、
「掛算のFFT化」とどういう意味　？

>>157
少しはさかのぼって読んでみたらどうか。>>147を読め。

うっ DIV 命令はダメでした。
理由はオーバフロー例外が起きてしまうからです。
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;

あ、とりあえず、
function U64Div32Mod(x:int64;y:longword):longword;
begin
　Result:= x mod y;
end;
とやればなんとか動きます・・・・滅茶苦茶遅いです・・・いつ終わるのか

　if n>= ends then break; を
　if int64(n)*int64(n)>ends then break;　にしたら

1234567890123481 1234567890123493 1234567890123527 1234567890123583
1234567890123647 1234567890123661 1234567890123671 1234567890123697

これくらいの桁なら少し待ってれば計算してくれました。けっこう楽しいです
でも、あと２桁が限界かも

スタンドアロンで計算するのもよいけれど、
分散システムのモデルを考えてみるのもよいかもしれない
今日このごろ。

>>162
それって強力にむずかしいんだよ。

>>163
難しいっていうのは、要するに小分けにされた問題だけ持って帰って、
計算結果を返してくれないクライアントの事とかを考慮していくと、
強力に難しくなるって事なんでしょうか？

問題を小分けにするのが難しいってことではないよね。

結局、どういう理論があるんだ？
誰かまとめてくれ。

あと％って何の演算子？

モジュロ演算子でしょ。
つまり剰余を求めるやつ。
VBだと、"Mod"。Rexxだと"/"。

じゃ、メルセンヌ素数について
2^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

Lucasテストに必要なのは　 (a^2-2)%b という計算だけ

つまりモジュラーな２乗と引算だけ用意すればいい

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が素数と判ります

モジュラーな掛算だけど、 2^nのモジュラーは上位桁を捨てるだけで良い

(2^n-1)のモジュラーの場合は 2^n　以上の値が出たら下位に回して足せばいい
例えば33%15= 10_0001b % 1111b = 上位の10と下位の0001を足して3 となる

だから掛算途中で2^n以上の桁が出たら、全部を下位にシフトして足せばいい

でもこれ、どうやって分散したんだろ？
まさか掛算？
真中の桁で分割して (a+b)^2=a^2+2ab+b^2で３人で計算して持ち寄るとか？

それとも、単に検査したい2^n-1のnを分けて皆で計算する？
運良くヒットした人が発見者？

離散フーリエ変換して各級数を持ち寄る、とか？

おっかしいなあ。TBitsで2^31個のビットを確保したら
それだけで落ちるなあ。

>>172
そらあんた、2^28 byte の連続したメモリブロックが必要になったんでしょ？
主記憶足りてる？
最初に SetSize で最大値を設定したかい？
俺んとこではうまく動くよ。

主記憶は
どのくらいだと足るの？

TBits.Sizeの設定はしたよ。
でも、計算中に突然落ちるんだわ。

SetSizeってTStreamあたりと絡めなければならないの？？

function　LucasTest(n:int64):boolean;
var 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;

----- つづく ----

var fsame,fzero:boolean;
var 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;

ほんの４桁の素数 n を入れたらもういつまでたっても計算が終わりません

そりゃ考えたら n^3のオーダーで計算時間が増えるからなあ、やっっぱり
FFTで掛算をさせて　せめて n^2 のオーダーに落とさないと・・・・

それでもやっぱり１００万桁なんてのは簡単じゃない

>174 計算中に落ちる時に何か報告ありませんか？
　　確保してない所に書いてたり読んでたり？

いま、 2^n-1 の n=2281付近で、１つの値の検査に１秒程度かかっています。
　今は掛算を筆算方式でするから n^3 のオーダで計算時間が増えます

掛算をFFTにすれば n*log(n)*n というオーダになるから
10秒くらいの時間でnが３万の単位まで計算出来ると思います。
目標が１００万桁なら、１日動かせば１つの値が計算出来ると思います。

ttp://www.geocities.co.jp/Technopolis/1793/mersennu.html
ここのページを見ると
　次の n は 400万台かもしれません　検査しなければいけない
素数は10万〜２０万ありますから、皆で１つの素数を分担して
検査するという方式でいいんじゃないでしょうか？

・・・・・その前に FFTで掛算するルーチンを作らないといけないけど・・・・

>>178
確保していないところに書き込んでも、その位置に自動的に size を
拡大してくれるから、落ちることはないんだけどなぁ。
逆に、少しずつ大きなところをアクセスするようなコードだと、
毎回巨大なメモリコピーが起きるから遅く見えるかもね。

ソースコード着いてないのか？

うーん、確保していないところに書き込んだりするような
ことは全くしていないのですが...
遅いとか遅くないとかではなく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;

n=4,000,000付近を狙うのだとすると
16bit単位にデータを持ったとして　4M/16=250K 以上のFFTが必要
（実数データだけなので半分のサイズのFFTでいいとして)

ttp://www.lohninger.com/fourier.html
２乗する為には２回必要だから262144 の2.1*2=4.2秒

これを400万回も繰り返させるには1年かかる・・・やっぱり計算そのものを分割して
する方法を考えないとダメみたいな・・・

Lucasテストを２分するアイデアその１

　　b=(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回求めて結果を比較すればいい

√(2*c) の求め方
a^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が奇数の間は２乗するより高速

×a'=√(a+2)ｊは　aが奇数の間は２乗するより高速
○a'=√(a+2)ｊは　aが偶数の間は２乗するより高速

でも何か変だなあ２乗の方も高速計算方法があるのか　>> 184が間違っているのか？

× a　　= c *2^((n+1)/2)
○ a　　= ±√c *2^((n+1)/2)

結論としては、 2^18=262144 のFFTを1/6秒で出来れば1日に１つの検査が出来る＝世界記録が狙えるわけね

この数字は無理筋？

FFTって複素数の掛算がだいたい 2*18*262144なんでしょ？
1回の掛算は４回掛算と２回の足し算だから
300MFLOPSあれば狙える数字じゃないかと思うんだけど？

というかプレステ２なら６ギガFLOPSあるから　プレステ２でやれば余裕じゃないの？
それとも精度とかの問題があるの？

31bitMax=約20億を計算しようとしたら
土日通して計算しても10億までしかやってなかった
ｳﾂだ。

メルセンヌ素数やルーカステストについて
詳しいページとかない？

数学の教養もなく頭の悪いｵﾚには理解不能な事ばかり。
ちょっくらやさしい教科書でもみて
勉強したくなった。

>>188 それはまさか1個単位に全部の素数で割る方法じゃないのか？
　　フルイ方式なら
>>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桁以上の複素フーリエ変換して
　求めた各要素を２乗した後、逆フーリエ変換すればいい

>>189
そうかーー、HDに書けばよかったか。
なんでもメモリ上にて計算するしか頭になかった。
それでメモリ壊れてOS落ちて悲しくなっていたのだ

HPのリンク先も見てませんでした

粗忽者です。

>>187 検算したら違ってるよ

次の賞金は１千万桁で１０万ドルか・・・・・１千万円だな。
これを狙うとすると、
２進で３千３百万桁 ルカステストをするにはこの回数ループしなきゃいけないから
1日で計算出来る為には１秒におよそ４００回のループ（２乗と１回の引算）が出来ない
といけない。

　16bitの内部表現として 2M WORD 2^21のFFTする事になる。
結局、プレステ２のFLOPSで64bit浮動小数点が出来ないと１日で１つの検査も難しい

>>189
>2^n^-1
これが顔文字に見えた俺は逝ってよしですか？

2^n-1が１千万桁になる素数nは
33219281
33219283
33219313
33219317
33219323
33219341
33219359
33219367
33219379　こんな調子であるぞ。

問題は、誰かがFFTを書いて数日で１つ検出出来るようになったとしても
それを分散してやった時に、賞金の配分をどうするかって所にもありそうだな
1万人が参加して一人千円じゃ配分の手間が大変だし、

宝くじ式にタマタマあてた奴が5割、ソフト書いた奴が２割、３割はサーバを提供
した奴ってのはどうだ？

サーバは、こんな感じだろうな

１）　メールアドレスとニックネームを登録してアクセスすると検査すべき素数Nが提供される
　　１０個くらい候補が出て選べるといいな
２）サーバーは登録された素数について検査中のマークを入れる
３）アプリは 検査が終わったらその値をサーバにマークする
４）ビンゴが報告されたらその値とニックネームが公開され、追試報告が何人か集まれば
　それをメルセンヌ素数候補として認める

一番の課題はイタズラビンゴ報告をどう防ぐかかもな

age

メルセンヌ素数について勉強中。

http://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 は素数である。

これってどういう意味？？
全然わからない。

やさしく教えてください。

>>196
7 = 2^3 -1 は (3-1)番目の項 14 を割り切るので素数。
15 = 2^4-1 は (4-1)番目の項 194 を割り切らないので素数でない。
31 = 2^5-1 は (5-1)番目の項 37634 を割り切るので素数。

という判定法。

わかったーー。

わかってみると。すっごーい。メルセンヌさんも>>197さんもえっらーい！

ありがとう〜♥

>>175-176 のコードの mul2をFFT化してみよう
　ttp://homepage2.nifty.com/m_kamada/docs/fftmul.htm
ここを参考にして考えてみる。

ただし俺は飽きっぽいし、なにより割り込みが入ったら途中で放棄するかもしれん

>>191 の言を信じて　2^21のFFTを用意する事にする
FFTは基数４の方が若干高速だけどメンドクサイから基数２で我慢しよう

実数->FFT->複素２乗->逆FFTという事は
周波数間引型でFFTして、複素２乗した後、時間間引型で逆FFT
すれば、途中、ビット逆順にする工程を省ける筈だ

周波数間引きと時間間引き
ttp://momonga.t.u-tokyo.ac.jp/~ooura/fftman/ftmn1_2.html

>>199 のサイトではA*B　を計算させる為に２組の実数列の同時FFT
を使っていた。
　しかし 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は精度を失う
これは問題にならないのだろうか？

>>201
そこだけ倍精度つかいな。

FFT の場合、最後に整数に戻る。ということを最大限に流用して
演算精度を考慮しないと、結局遅くなるよ。
いかに計算の手を抜けるか、を理論的に抑える努力が必要。
で、そこら辺は門外不出の秘伝に近いものがある。

Super PI のコードが公開されないのもそういった理由だと思ったが。

どうやら、半分のサイズのFFTはそれほど効率的ではないようだ
せいぜいまともにやった時の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で割ればいい)

>>202 倍精度っていってもlong doubleで doubleより16bit増えるだけだから
・・・・でもやっぱり16bitの２乗は32bitだし
N桁の２乗は中間桁ではN倍されるから32+21=53bit必要だから
やっぱりlong doubleでも ダメっぽいね

8bitにするしかないな・・・結構厳しいなあ

で、2^K-1のモジュラーの２乗を どう計算するかだけど
Ebitの所で分割する。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でやらなければならない

double で計算するのと Extended(long double)で計算するのって
やっぱり速度上違いが出るのだろうか？
doubleが早いなら FFTはdoubleで計算して ２乗後は
Extendedで計算ればいいんだけど

とりあえず今日 出来た所まで

まず、ビット反転ルーチン　今回は必要無いかと思ったんですが
実数２つの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;

procedure fft(n:integer;var x,y:array of Extended);
var 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;

//fftの逆演算ただし fft->ifftにより 2^n 倍される
procedure 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;

試験用に書いた x*y を計算するコード

procedure 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;

で？　結局　ルカステストにはどう組み込んだらいいのだい？　あげ

ルカステスト n=3〜19迄
　　　　　19　　 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

細かい訂正 >>203
×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

コードはとりあえず全部書いたけど、バグだらけでデバック中

とりあえず、必要だったルーチンを掲載

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;

いまさらなんですが

エラスト？エストラ？
そのフルイの配列確保がオンメモリでは出来なさそうなので
ＨＤにアクセスしながら計算させたいのですが
ＨＤ直接(?)アクセスする方法がわかりません。

教えてもらえないでしょうか。

環境は...Delなんですが...

>>215
　　HDにアクセスと難しく考えなくても　そのまま
メモリを array of byte で取って SetLengthすれば自動的に
仮想記憶が働いてハードディスクをアクセスするよ


>>205 で書いた方法で
E=(K+1)/2とせずに、 8の倍数にすれば正規化処理で少し
楽出来るけど、計算精度が 最悪9bit余計に必要になるので
Extended の限界を超えてしまう。

//使ってる変数
var 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;

正規化 ---これもどうも怪しい
　FFTによる掛算結果をそのまま再度２乗するとデータ精度が保て
ないので、各桁を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乗は正常に働く筈）

{モジュラーな２乗用ルーチン
E 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;

>>194
検査すべき素数を提供する時にパスワードを生成して渡しておくんじゃダメ？
そのまま検査された値全部にフラグとパスワードを保存するとディスクの無駄になるけど。

>>215です。
>>216さん

その確保したarray of byteの
Xビット目が0なのか1なのか判断する方法が
わからないでございます。

どうやるのでしょう。
質問ばかりですいません。

var tbl:array of byte; とすると

X ビット位置が１かどうかは

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 かどうかで判定して除けばテーブルは半分で済む

まだデバックが終わってないから判断には早いかもしれないけど

速度的に面の見積もりは出来た
k=2281付近で２５秒
k=4423付近で120秒

世界記録を狙うとなると１個の判定に１年ですまない計算だ
２桁以上計算速度を上げるアイデアが必要

× n:=iLog2(m);
○ 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
の表の通りの判定が出来るようになった

汚いソースになったけど、誰かいるなら公開するけど

http://www.geocities.co.jp/SiliconValley-Oakland/8522/data/LucasTest.txt

プログラムが正しければ世界記録を狙った場合、
2ヶ月程で１つの数の判定が出来るだろう。

ただ、調べなければいけない数も膨大なので
出来れば1/50　せめて1/10に縮めないと・・・

>>226
100堕慰くらいで実行スレばいいよ。
そんなに2chプログラム人口いないか？

今 Extendedで計算してるけど Doubleで良ければ速度は倍以上になる
世界記録は22bit付近だから ワード分割サイズをAbitとしたら
2A+22bitでいいならDoubleで足りる
2(A+22)bit　必要ならExtendedが必要になる

>227
1bit単位に分割してFFTで掛算すれば何ビットの精度が必要だろう？

もしそれが22bitでいいのなら、PS2で実行という方法もある

でもメモリが 2*4*4Mbyte必要になるけど

遅くなりましたけど>>222
ありがとー。

>ビットを立てるなら
>tbl[X div 8] := tbl[X div 8] or (1 shl (X mod 8) )

ビットを落とすのはどうしたらいいのかなあ。
論理演算とか苦手だわ
....notしてからandするのかな？

>>230
ビットを落とすのは not して andでもいいし
すでに立っているのが判ってれば xorでもいいよ

めった、簡単なこと質問してスマソですわ。

>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)) )

つうことですか？

他にもリセットの方法はあるよ

tbl[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)) )

実際に浮動小数点で何bit あれば判らないので
Single = 32bit 仮数部23bitにして計算させてみた。
Doubleの結果と比較して k=251の時に差が出た

これから浮動小数点ならほぼ 2*bit幅+log2(k)が限界
と考えていいように思う。

つまりdouble (仮数部52bit) のままで計算可能ではないかと

であれば、>>224の値はほぼ半分になる

毎度ありがとうございます。

Indexビット目が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 );
すると大丈夫みたいっす。

メモリ確保量は同じ気がするのに、なぜだろう。

世界記録。がんばってくださいマンセー
トロイけど、このスレ呼んで勉強しーとこ。

こうやった方がいいかも

Result := (Tbl[Index div 8 ] and (1 shl (Index mod 8) )) <>0 ;

さらに

Result := (Tbl[Index shr 3 ] and (1 shl (Index and 7) )) <>0 ;

>>236ﾀﾝさっすが！

これなら世界記録付近の計算は１件２週間〜１ヶ月程度で出来そうな感じだけど
>>71の賞金を貰うならさらに　１００倍かかるな
100倍になると、今のままだと8年はひとつの数の判定にかかる計算
それにメモリもオンメモリで出来る限界付近だからもっと速度は落ちるかも

やっぱり1千万の賞金は簡単じゃないな

にしても、ちかずいておりますが...
並のパソでそれだけいくと
もちっと速いコンピつかえば＆
数年でクロック数上がれば、どんどん賞金に知被いて....

素数の歴史上、大学のコンピを使えた高校生が
最大素数を求めたことがあったらしいけど
アルゴリズムだけOKなら大型コンピュータにて
流すことはできませんか？

>数年でクロック数上がれば、どんどん賞金に知被いて....

条件は同じだから、差をつけるには　やっぱりプログラム技術の差でつけないとね

それに、スパコン借りるなら成果出さないといけないだろう。
今の速度じゃ、スパコン１台では見つかるとはとても思えない

計算速度を上げるには、やっぱりFFTの固定小数点化しか
ないんじゃないかな。固定少数点なら、今のCPUなら足算は２〜３並列化
されるから巧くコーデングされれば相当高速になる筈。

データを64bit固定少数点にして、サインのテーブルは32bitに
すれば、掛算のコストは倍、たし算は３倍になるけど、
16bit単位に区切れるから計算量としては同程度

それにしても、なんで神様は　2^16+1を素数にしたのに 2^32+1とか2^64+1を素数に
しかなったんだろ

へたに固定少数点化して命令数を増やすよりは
128bitSIMDにアセンブラレベルで対応した方が速いんじゃないかな
１クロックで４並列演算とか出来るのは、あれは単精度(32bit)の場合だっけ？

64bit浮動小数点を SIMDで実行出来る=SSE2 は Pen4からの筈

AMD 3DNow
http://www.amd.com/japan/K6/k6docs/j21928c.pdf
Pentium 2,3
ftp://download.intel.co.jp/jp/developer/jpdoc/iaopt_j.pdf

今、k=11213 付近判定させると 15秒で結果が出る
１回の計算 は 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は魅力的だな・・・でもまだ買い換え予定はないし

出来れば整数演算でどれくらいでるか試してみたいな

今2^NのFFTは模式的に書くと
for 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の順にネスト
させたい。どうループを書けばいいんだろ？
再帰を使わないと書けないかな？

このページを見つけました。
ttp://www.pluto.ai.kyutech.ac.jp/plt/matumoto/fftb_export.c.txt

キャッシュを考慮してバタフライ演算が入りきるサイズに分解してるようです

結局、こういう感じです。　再帰を使わなくてもいいけど、Delphiは
Cの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の方は逆に１ブロック実行してはネストしてゆく格好で書けました

整数化する為に、
１）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終了時には整数化されている

と考えた段階でゲップ状態・・・・メンドクサイと思ったら終わりかな

×$7fffffff
○$3fffffff

整数でやるんだったら、有限体に値を持つ FFT なんてどうでしょう。
それとも実数でやるのがミソなんでしょうか。

ps2linux でやってみようと思って考えてたんだけど、
vu には 32bitx32bit の乗算なんてないのか(当然か...)
12x12 ならできるけど、EE はこんなのにはむいてないみたい。
今更だけど。

ここにある方法ですね　>> 249
ttp://www.geocities.co.jp/WallStreet/6351/mersenne.html

勉強してみます。
理論が良く判ってないのと、
フーリエ変換が他への応用も出来るのと浮動小数点を使う
限りは簡単だから拘っていました。

でもFFTを固定少数点でやるくらいならそっちがマシなような気もします

>>250 のページでは 2^1024+1や 2^31+1 を法にしていますが
2^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
と非常に簡単に計算出来ますよね？

>>251
n 乗根を考えるということは
もとの数を 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 で適当に作ってみたらものすごく遅かった。

賞金を狙う事を考えると 30,000,000 bit以上に目標が必要です。
60*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のマシンでもルカステスト１回にやはり年単位は必要
という事になってしまいます

これを短縮するには、この変換操作を数１００台のマシンで分解
して計算させるアイデアがやっぱり必要になるようですね

>>252
>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の場合よりもさらに簡単です(負数を考えなくていいから)

と考えたのだけど、ホントかどうかはまだ試していません

>>254
ちょっと考えてみたんですけど、
素数でないとやっぱり苦しいんではないでしょうか。
フーリエ変換と逆変換が実際に逆になることを示すのに、
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 となってます。

何かうまくやる方法があるんでしょうか ?

ホントですね。FFTは出来るけど、逆変換が出来ないですね
15を法にしてみると　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の逆数が無いから解けないですね

やっぱりダメか・・・・

あ、でも、それぞれ
|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倍になる解なら同じ規則で見つかりました

ただ、この調子だと実際に計算出来るビット数はだいぶ小さくなってしまいますね

見つけた規則だと 2^(2^n)-1で 2^nのFFTをすると
2^(n/2)-1 倍された逆FFTならできるという事になります

32bitの場合は 65535倍される訳で、これは約数になっていますから
つまり半分のビットが使えないという事みたいですね

あ、そうか 2^(2^n)-1は
11*101 =1111
1111　* 10001=11111111
11111111　　* 100000001= て調子で 2^m+1 との掛算なんだ

じゃ、並列処理をどうやるかって方向でひとつ

High(Integer)までの素数列を
常識的な速度で出す事やOnメモリでダス方法はわかったのだけど

それをファイルに落とす場合に
for i=1 to High(Integer)
　IF iが素数なら
　　S := S + IntToStr(i)
のIntToStrが超低速で参ってます。
なんとか高速に沢山の数値を
IntToStrと同様に文字列にしてしまう方法はないでしょうか？

ファイルの保存形式が別にStringじゃなくても
いいのでIntegerのメモリイメージがそのままファイルに
落ちてくれるといいのだけど...

integerのメモリーイメージのまま出力したいなら
ファイルマッピングを使ったらどう？

ファイルを開いて CreateFileMappingで割り当てて
　MapViewOfFile　でポインタを貰って
ポインタ使ってInteger配列として操作出来るとか

篩データをこんなのにして持っています
　TByteArray = 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
に代入してから保存した方が使いやすそう
(要はファイルサイズの問題)
かと思うのですが

その場合のファイルマッピングの方法を
もう少し教えていただけないでしょうか。

つーかこんなとこでうだうだやってるより
さっさと分散＆ブルートフォースで解決

ああ、まだ初心者だったんだね。>>263

　１） 動的配列じゃなくて、配列のポインタにする事
　２）　そのポインタの実体をGetMem　とか GlobalAlloc で確保してみる
　３）ファイルマッピングを使ってみる　それをコンポーネントにしておくと
　　　ちょっと便利かな

　　というようにバージョンアップしてみましょう

32bit程度ならまだまだ　ビットイメージで持っていた方がファイルサイ
ズの面でも速度面でも有利だよ。

ただ、32bitくらいの篩なら保存するよりはその都度作った方が速いかもしれないね

>ああ、まだ初心者だったんだね。>>263

　グサッ　　Λ＿Λ
　　　　　：（ﾟД｡ ；）
――――（”∴ 　）　グヒャッ・・・
　　　 　．：｜｜　｜
　　　　∴ （__（＿_）

そうか・・・初心者だったのか。

>　１） 動的配列じゃなくて、配列のポインタにする事
>　２）　そのポインタの実体をGetMem　とか GlobalAlloc で確保してみる
>　３）ファイルマッピングを使ってみる　それをコンポーネントにしておくと
>　　　ちょっと便利かな
うーーーん
全然イメージがつかめません。
具体的なコードが見られる所ってないですか？webでもマニュアルでも
動的配列って実際はポインタじゃなかったです？

>32bit程度ならまだまだ　ビットイメージで持っていた方がファイルサイ
>ズの面でも速度面でも有利だよ。
2^30bit=2^27byte=2^17KB=2^7MB=128MBは必要なんだけど
速度面はともかく
2,3,5,7,11,13,....という配列とくらべて
ファイルサイズは有利じゃないでしょ。

>32bitくらいの篩なら保存するよりはその都度作った方が速いかもしれないね
Myマシンが非力なのでメモリ的につらいの。

> 2,3,5,7,11,13,....という配列

素数定理を使って32bitの素数の個数を 2^32 / log(2^32) 個と見積ると、
2^32 / log(2^32) * 4 Byte ≒ 195000000 Byte ≒ 186MB 必要。
ファイルサイズでも不利だと思われる。

うーん　Delphi3くらい迄は動的配列なんてなかったから

こっちが当然の書き方だったんだけどなあ
http://www.nifty.ne.jp/forum/fdelphi/faq/00160.htm

nBITの素数を数値として記録するには
n* 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の篩を用意しておけばいいのでは？

>>267
サンキュウ。
そうか、素数定理というのが抜けてたのか。
その定理よくわからないけど、検索してみる。

だけど、Integerって最大は2^31だから
もう少し減るかも。

なんか、自分の数学的素養が著しく
失われてることに気が付いたので
http://www.geocities.co.jp/Technopolis/1505/log.htm
ここから読んでみます。
プログラム初心者というより
人生初心者か？義務教育初心者....(TДT)ｳﾂﾀﾞ

>>268
D3ってもうずいぶん前なので､､今ではそのやり方のほうが
トリッキーかも。

>>269
ごめん。全然わかんないや。
素数定理というヤツかな。

>32bitの素数が欲しい場合、16bitの篩を用意しておけばいいのでは？
そうなの？そうは思えないんだけど、

ごめん log(2)はおよそ 0.69 0.3は log10(2)の場合だった

だから 1/log(2)は 1.44　だった　まあどっちにしてもサイズに関係なく
単純に記録するなら　常にビットマップで記録した方がメモリ容量は少ないって
事だね。　素数は案外詰まってるんだよね

age?

素数定理ならここなんてどうよ。
http://www.com.mie-u.ac.jp/~kanie/tosm/humanind/yosou.htm#prime_number
数学は科学の女王で、数論は数学の女王だ。
頑張れ。

ζと対話できるようになれば一人前だ。

夏休み中に賞金狙うやついる？

>>274
無理にきまってんだろ。

100万人くらい参加したら可能かも・・でもどやって賞金分ける？

1千万円くらいあるんでしょ？
百人くらいを一つのグループにして、そのグループで一つの数のチェックを
担当して、　ヒットしたらそのグループ＋運営側で半々に　わけるってのは？
　百人なら半分でも５万円になる

100人いても９８人は厨房と思われ

プログラム作った人が総取りすればいいだろ。

>>279 それじゃ誰が何の為に協力してくれるんだ？
感動的ストーリーでも捏造するか？

じゃあ、白血病やSETIなどで協力している奴はなんなの？

単にマシンパワー提供する程度で金が貰える
と思うバカは必要なっしんぐ。

俺のP-133も役に立てるかな？

>>282
たたん。最低でもP-500以上。

ニュース速報板に紹介されてるよ。
http://kaba.2ch.net/test/read.cgi?bbs=news&key=997090075&st=112&to=112&nofirst=true

Internet Mersenne Prime Search
http://mersenne.org/PrimeNet/

素数が無限に存在することの証明として
>もし素数の数が有限個だとして、それを小さい順に ｐ１、ｐ２、ｐ３、・・・、ｐｍとします。
>新たに、(ｐ１・ｐ２・ｐ３・・・ｐｍ＋１)なる数を考えます。
>この数は、任意の素数ｐｎで割れませんので、素数です。
（一部抜粋）

もし、素数の数が有限個だとして、それを小さい順に ２，３とし
新たに、２Ｘ３＋１なる数を考える。この数 ７ は素数。

もし、素数の数が有限個だとして、それを小さい順に ２，３、７とし
新たに、２Ｘ３Ｘ７＋１なる数を考える。この数 ４３ は素数。

もし、素数の数が有限個だとして、それを小さい順に ２，３、７、４３とし
新たに、２Ｘ３Ｘ７Ｘ４３＋１なる数を考える。この数は素数。

これをひたすら繰り返していけばいいのでは。

>>286
これって…　マジで
できそうでないかい？

2*3*7*43+1=1807 1807/13=139 (；ﾟДﾟ)

>>287
ｶｺﾜﾙｲ

> もし、素数の数が有限個だとして、それを小さい順に ２，３とし
> 新たに、２Ｘ３＋１なる数を考える。この数 ７ は素数。

> もし、素数の数が有限個だとして、それを小さい順に ２，３、７とし
> 新たに、２Ｘ３Ｘ７＋１なる数を考える。この数 ４３ は素数。

ここがおかしい？

素数を小さい順に 2,3,7 とするのはダメ、5 も素数だから。
小さい順に 2,3,5,7 としないといけない。

だから「この7は素数」と分かったときに
3から7までの数で素数を探す必要があるのでは

> だから「この7は素数」と分かったときに
> 3から7までの数で素数を探す必要があるのでは
これが大変すぎるから、このスレがあるんだよね…
ageてゴメンヨ…

(* 2 3 5 7)
210
↓
211 は素数

(* 2 3 5 7 11 ... 211)
1645783550795210387735581011435590727981167322669649249414629852197255934130751870910
↓
1645783550795210387735581011435590727981167322669649249414629852197255934130751870911は素数

３０分走らせたけど、>>292が素数かどうか判定出来んかった・・・

> もし素数の数が有限個だとして、それを小さい順に ｐ１、ｐ２、ｐ３、・・・、ｐｍとします。
> 新たに、(ｐ１・ｐ２・ｐ３・・・ｐｍ＋１)なる数を考えます。
> この数は、任意の素数ｐｎで割れませんので、素数です。

そもそも↑これ、証明できるか？
その数が「Pm より大きくて (P1・P2・P3・・・Pm ＋ １) より小さい数」で
割り切れる可能性は否定できないと思うが。そしたら素数じゃない。

>292
ちょっとしたことをオンラインで検証できるツールがあるぞー。

■桁数無制限電卓

　http://www8.big.or.jp/~000/CyberSyndrome/rsa/calculator.html

■巨大数の素数判定のアプレット

　http://studwww.rug.ac.be/~bcallens/project/priem.html

ごめんsageてもうた

　>>295 のアプレットに 292の数字を入れたら is composite だそうだ

それにしても速い。　これがRabinの定理？

>>297
うむ。
　１．ふるいで素数を小さい方から1000個生成
　２．とりあえずその1000個で割ってみる
　３．>>25で言う"任意のa"として小さい方の素数50個を使ってRabinテスト
…というアルゴリズムなように読める。

>>294
背理法って知ってますか。
これは素数の数が有限個だとして矛盾を導いているんです。

>>299
命題の対偶という概念が分からないと
背理法で何が導き出されるか分からないのでは。

>>299
んだから、素数が無限個であることの証明としては問題ないのはわかっちょる。
だからといって、素数を無限に生成する式としては使えるとは限らないだろ？

『 「素数の数が有限個」の仮定のもとでの非合成数 』 なら
確かに作れるけど、少なくともこのスレではそんなもんイラン。

>>294

なんか誤解してたみたいだ。スマソ。

保存アゲ

あげ

>>303-304
２ヶ月もなにもないって事は、
終ったってことさ。
それにdat落ちしても、このスレのリンク覚えておけば参照できるよ。

そうは言っても　この前の２ＣＨ転送量超過騒ぎの時は
過去ログについてもいろいろ　議論出てたからなあ・・・・・・
まあ落ちたら落ちたでしょうがないとは思うけど

だれか他のいくつかの板みたいに
過去ログ保存（バックアップ）しておくサイト
つくってくれないかなあ。全文検索できるようにしてさ。

pizaの方なら見られるけどpiza2はパスワード付きになるみたいだね
このまま落ちると　過去ログとして見られないのね

このスレにはFFTのソースがあるので勿体無い

24の方法が一番速そうでない？
調べたいウン万桁の数をメールで送る。
数秒後か数分後か数日後、結果が返ってくるのでは？

保守アゲ

>>309
それはwebサービスってことですね

面白いね、素数判定。
ためしに作ってみたら
1から1千万までの数を判定するのに80秒かかった。

>>312
CPUは？
試しに作ってみたら、P3-850で14秒だったよ＞1千万まで

最速の素数判定アルゴリズム

速ければRSAとかさくっと。

>>313
P3-866 (^^;
少し改良したら1千万まで平均25秒にまでには短縮できたＹＯ！

ムズイ。25秒の壁をやぶれん。ビット演算とかしなきゃだめ？

>>314
素数判定が早くても RSA は破れないぞ。
必要なのは素因子分解

>>316
画面表示に時間がかかってるだけじゃない?
ファイルに出力するとセレ566で2秒くらいだけど…

>>318
えー、一千万まで2秒ですかあ？
僕の奴だとファイル出力すると26秒だよ。
さらに画面表示まですると、ずっと遅くなる。
言語の問題とかも有るのかな。Javaで組んだんですけど。

全部計算してからファイル出力。これで高速化だ！

↓ソースです。
http://pc.2ch.net/test/read.cgi/tech/1003569074/765-768

limit = ルート(target)
の方がいいんでないかい？
あと2以外の偶数は調べないとか

>>322
>limit = ルート(target)の方がいいんでないかい？
これで9秒削れて

>あと2以外の偶数は調べないとか
これで1秒削れました(^^

15秒程度になった！ありがとです。

作ってみた。これで約２６秒（スーパーπが１３分のマシン）。
int 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;
}

エラトステネスの篩っていう有名なアルゴリズムがあるよ。
割算を使わないのでかなり速くなるはず。
メモリたくさん使うけど常駐するわけじゃないし…

# てか1000万って篩に使うメモリを意識して選んだのかと思ってた。

>>325
おお、劇的な変化が！
一千万までの演算が15秒から一気に1秒台へ^^;

>メモリたくさん使うけど常駐するわけじゃないし…
メモリ食いますね。
ていうか1億で計算しようとするとメモリが足りなくなってプログラムが
異常終了してしまふ。改良せねば。

64Mのメモリがあれば 1億=100M は大丈夫でしょう

多分、やり方がアホなんですよね。
boolean[] b = new boolean[100000000];
で落ちてしまふ。

この配列で、素数ならばtrue、合成数ならばfalseを入れるように
つくってるんですが。

　 ∧,,∧　　二日よいで氏ぬ
　ミTДT彡
　 (∩∩)

最近のコンパイラでは
演算の高速化のために
booleanは実質32bit変数だから
TrueとFalseを判定するのに
すごく余分(32倍)にメモリを食っているのではないでしょうか。

すると、
int[] b = new int[3125000];
とかやって、フラグを立てるときは、
b[i / 32] |= 1 << (i % 32);
みたいなカンジでＯＫと。

>>329
なるほど。

>>330
おお、それを使って修正してみます。

つか、いつの間にかすごく志が低くなっている…

　　∧,,∧　 　　／￣￣￣￣￣￣￣￣￣￣￣￣￣
　 ミ,,ﾟДﾟ彡　＜　一旦、sagaって終了したスレだし
　 ,;ﾞ　　ﾐ　 　 　＼＿＿＿＿＿＿＿＿＿＿＿＿＿
　ﾐ.　　ミ
／゛゛゛゛
Cマガで素数問題が出題されたから、復活ageされたのかな？

まあ「FFTで128bitオーバーの数値を表現」とか素数判定アルゴリズムとは
直接無関係なアルゴリズムで盛り上がってたからな。
つか「素数判定はエラトステネスの篩」っつー結論は>>4あたりで既に出てる
わけだ。名前間違ってるけど。

>>329のカキコを見て反射的にコードを考えちゃったけど、
そもそも、素数を求めてどうするんだろうか。
ハッシュ表のサイズを動的に変更する……なんてことはしないよな普通。

>>334 FFT は　掛算の高速の為で
　　高速な掛算が必要なのは　数千、数万とかものすごい桁の素数を検定しようとしたから

>>336
んなこたわかってるって。だから>>334にもそう書いてあるだろ。128bitオーバーって。

そうか。　そりゃスマン
　FFT掛算が有効なのは数千桁以上だから 128bit見て反応しちまったよ

>>338
つーか、まず組込型で演算できない部分を計算するのにFFT使う、ってことで
さらに精度を保ちつつ係数を2^nにして高速化とかそんな感じで進んでたと思
うんだが。
数千桁以降じゃないと有効じゃない、と言われても、数千桁までは何でやるんじゃい、
というのは特に議論されていなかったはず。

というか、仮想空間（４G)全部使って何桁の数字が表せるの？
log10 2^(8+10+10+10)*4
=152log10 2桁
あってる？

>340
たしか、686以降では36bitアドレスバスだから、
Windowsにこだわりが無ければ、64GBのメモリが使える筈だけど。
実際に対応OSがあるかどうかは知らない。

>>341
物理的に4GB以上メモリをのっけたとしても、ポインタは依然32bitだ。

>>339
FFTはあくまでも高速化のためで
数千桁以下は、 確か100進で計算させてたんじゃないかい
と思って読み返したら >>175-176で16ビットで計算してるようだよ

　読み返してみたら結構面白かったよ

>>340
log(2)*8*2^32 = 10343311891.9
わずか10G桁=１００億桁

>>342
ということは、１つのプロセスが使用可能なメモリ空間は、
ＯＳが何であろうと４ＧＢまでってこと？
アセンブリ言語は１６ビット時代までしか知らないんだけど、
セグメントレジスタは３８６で拡張されてませんか？

>>345
一つのセグメントあたり4Gだから、
一つのプロセスが複数のセグメントを使えば
それだけメモリが使えるよね。
OSがそれに対応してるのかとか、
そんなコードをはくコンパイラがあるかのかは知らんが。

>>346
少なくとも、WindowsとLinuxは、４ＧＢが全実装メモリの上限だったハズ。

なんか、4GB以上のメモリを扱えるWindowsが有るみたいなんだけど。
どっかの企業の改造版だったかも。

>>347
物理メモリ上限は、Win9x系は2G、NT系は通常2G、オプションで3G。
これを越えるのはまだ無いと思われ。

>>348
改造しても(規模にもよるが)4Gだと思われ。
それ以上となると、ちょっとやそっとの改造じゃ済まない（つか、１から作った方が速い）

マジで！？
１プロセスでもオーバー４ＧＢだったら凄いな。

よく考えたらLinuxはPAEに対応してるよな。
物理メモリは64GBまでいけるね。
一つのプロセスでは4GBまでだとおもうけど。

>log(2)*8*2^32 = 10343311891.9
>わずか10G桁=１００億桁
なるほど，ミスったよ
1G=1000MByte=1000*1000KByte=1000*1000*1000Byte=1000*1000*1000*8bit
=10^9
でもって
表現可能な数字＝2^(10^9)
ってわけか
で、こいつのlog10をとるのか・・・
そうすると，大体一桁落ちて
数千万桁
あれ，どっかで計算ミスってる？

まあ高速なアルゴリズムよりも
高速な分散処理を考えないと普通に無理だよね
分散って言っても，一つの数値の判定を分散処理する
数学板から誰か連れてきてよ
でもって作る？
通信するので、言語は何でもありですね。
まあ、C++ですかね。

>>349
手元にあるWindows2000 RC2の説明によると、ProfessionalとServerは4GB、
Advanced Serverは8GBまでサポートしてるよ。

素数関係ないのでsage

>>352
でもな〜、素数の判定にロハで協力すれ、って言って、人が集まるのかな。
SETI と比べてジミかつマニアック過ぎると思われ。
100万桁の自然数ひとつ判定するのに、どれくらいの台数のマシンが必要
なんだろか。
一つの数の判定するのにかけていい時間を1日くらいと考えると･･････
数十万台はいるんじゃなかろうか。

>>354
かっこいいスクリーンセーバー作って、裏で勝手に計算しちゃえよ（笑

>>354
自己複製機能と感染機能を実装すれば大丈夫かと。

nimbda改造すればすぐできそうだNE!

>>334
>「素数判定はエラトステネスの篩」っつー結論

デマを撒くな〜

>>25
の確率性アルゴリズムを使うか、
特定の構造をもった数に対する検証のどちらかしかないぞ。

スレのテーマが「最速で素数表を作る」ってなら櫛で間違いないけど。

>>352　log10(256)が2.4だからバイト数のおよそ2.5倍ね

4Gバイトなら10G桁表現出来るって事ね

39番目のメルセンヌ素数が見つかったらしいな。
インターネット上のPCの共同作業で。
http://mersenne.org/status.htm

オンメモリじゃなきゃもっと行くでしょ？
遠いネットワークの先にあるメモリーよりも先に地元のディスク使わないの？

1億までの素数リストを出すのに約20秒かかった。
10億までなら200秒かなとか単純に思ったら、すんげえ時間かかる。
ハードディスクがガリガリ言ってるし。
10億でもたったの10桁・・・素数計算って大変だねー。

>>360
うわ。見逃していました。すげーニュースだ。
まだ検証中（ダブルチェック中）なんだね。
M38 が (2^6972593) - 1 で M39 が 8000000 以上ってことだから
かなりの値だ。

>>1
なに、Ladinの定理?

よくわからんが１０進数で何桁だ？

2^8000000くらいなら10進で2408240桁くらいだろう。

>>361
１Gのデータをハードディスクから読み出すのにどれくらいかかると思ってるん
だか・・・・

うちのマシンには12ディスクのRAID(FC-AL)が8並列で付いてるので
1GBの読み込み1.3秒ですが何か？

PCIバスの帯域は32bit,33MHzだと133MB/secですが何か？
(64bit,66MHzで533MB/sec)

そんな安いバス使ってませんが何か？

crc

最大の素数発見！ 13万人で検算中。
http://news.2ch.net/test/read.cgi/newsplus/1007792330/l50

あ、検算済みでした。

age

>>64
できました。
http://www.research.ibm.com/resources/news/20011219_quantum.shtml

>>64 >>375
こっちのほうが見やすいと思われ。
http://slashdot.jp/article.pl?sid=01/12/20/155216&mode=nested

/*素数のlucas判定のプログラム作ってみました。
　C言語による最新アルゴリズム辞典の改造しただけだけどどうでしょう？*/
#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);
}

素数を見つけるよりもロト６の１等の番号の組み合わせの方が知りたい
とフト思ってしまった

　エラトステネスの篩を使用して、ゴールドバッハ（ゴルトバッハ）の予想を、
1億までエクセルのＶＢＡで計算しましたが、実行時間は、約２０時間、反例は
見つかりませんでした。
　これって、結構早いのでしょうか？
　ちなみにマシンは、ＶＡＩＯのPCG-FX55/XPで、メモリーは６４ＭＢ。
　なお、ゴールドバッハの予想を知らない方への説明。
　全ての偶数は、２個の素数の和で表される。
　そこで、ある制限をかけて、いかなる偶数の平方根の２倍＋１よりも大きい
素数と、その他の素数の和となることは無いと予想したのですが、果たして反
例は見つかったか？結果は？
　興味のある方は、掲示板に書きこんで頂けば、マクロが入ったエクセルのソ
フトの配置してあるＵＲＬを書き込みます。
　なお、今のところ、私のエラトステネスの篩は、１００万の計算に１２秒かか
りますので、まだ改善の余地有り。
　なお、定義に従えば、１００万の計算に８日かかります。

ところで、ある程度以上大きくて2^n-1以外の素数ってあるの？
ある程度ていうのも、へんな言い方だけど…

＞379
VBで実行してみたいです。

>>380　 素数は Aの付近に 1/log(A)程度の頻度でありますよ
log(2^n-1)はおよそlog(2)*Nだから　2^n-1〜2^(n+1)の間にも0.7*n個程度あるという事

>>382　さんきゅ。

×　だから　2^n-1〜2^(n+1)の間にも0.7*n個程度あるという事
○　だから　2^n-1〜2^(n+1)の間にも(2^n)/(0.7*n)個程度あるという事

　ＦＦＴではなぜ、ビット反転操作をするんですか?教えてください。

刺身定食くいたくなった…

>>385
ここ参考にするといいでしょ
http://momonga.t.u-tokyo.ac.jp/~ooura/fftman/

特にこの絵
http://momonga.t.u-tokyo.ac.jp/~ooura/fftman/fig_bf2.gif

最後にビット反転置換をしなくても良い方法もあるけど、
　すると毎回複製が必要になるのと、ループが少し複雑になるから
メモリの節約になるこの方法が好まれているのでしょう

ありがとうございました＞387

確率的判定法(Rabin)の計算量オーダーはlog(n)なのは知ってるけど
確定的判定法で最速の物ってどれよ？
ecppだっけ？あれはどれぐらいのオーダー？

>>389　なんでsageで聞いたの？

俺も知りたいage

あげ

あげ

確率的判定法って　単に何度か割ってみて割り切れなかったらと変わりないように思うのだが？

例えば1万個の素数を用意しといて割ってみて割り切れなかったら　素数でない確率は1万分とか

>>394
オオバカ

>>394 確かにそうだけど、それだと1桁確度を改善しようとすると、判定時間も１桁増になってしまう

Rabinの方法は、対数和にしかならない

保存

ぷろろｒぐ

素人でよくわからんのですが
テーブルを分解・圧縮して計算に使う部分をなるべくキャッシュに載るようにして
レジスタも使って（速度の桁が違うので）
2の倍数乗除にシフト、（プロセッサ内で勝手にシフトに切り替わるかな？^^;）
遅い割り算を逆数掛け算にしたりとかは効果ないんでしょうか？

倍数→累乗　マチガえました(笑
大規模３Ｄレンダリングなどで使う
クラスタリングが分散にいいらしいです

>>399　小手先のテクニックは　同じ計算時間で　計算出来る数が少し増えるだけ
　並列計算も、完全に成功しても台数分　数が増えるだけ
　それは頑張って１０台並列でももたった１桁増えるかどうかのレベル


　掛算をFFTにしたような、一挙に桁数を倍に出来るようなアイデアが求められている

32bitMAXなら計算済みのを展開して2分探索すればあっという間だね