5ちゃんねる ★スマホ版★ ■掲示板に戻る■ 全部 1- 最新50  

■ このスレッドは過去ログ倉庫に格納されています

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

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

2 :デフォルトの名無しさん:2001/06/25(月) 18:27
>>1
理解できないのならあきらめるしかないでしょう。
他人をたよるのはよくありません。

3 :デフォルトの名無しさん:2001/06/25(月) 19:22
>>1
その「さっぱり」って言葉がやるきねぇなぁと他人に感じさせるので、
そういう聞き方では知識ナルシスト以外教えてくれないと思われ。

4 :デフォルトの名無しさん:2001/06/25(月) 19:26
エラストテネスの篩で我慢しろ。

5 :デフォルトの名無しさん:2001/06/25(月) 20:15
素数なんて、計算するものじゃなくて
データとして保持しておくものだよ。
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

6 :デフォルトの名無しさん:2001/06/25(月) 20:17
>>5
神父もびっくりだな

7 :デフォルトの名無しさん:2001/06/25(月) 20:19
下のほうは興味あるなぁ。
日本のチームはあるの?
京都産業大学なんか大学でチームつくってそうだな。

8 :デフォルトの名無しさん:2001/06/26(火) 01:09
100万桁の
100万って、どのくらいなんだ?

100万文字ってなんKB?

9 :デフォルトの名無しさん:2001/06/26(火) 01:17
百万桁と百万文字は違うだろ

10 :デフォルトの名無しさん:2001/06/26(火) 01:18
>>8
BCDで持つとして、100万桁=50万バイト=5000キロバイト=5メガバイト

11 :特殊早大生理論(広末涼子):2001/06/26(火) 06:50
量子コンピュータを作って計算させる

12 :デフォルトの名無しさん:2001/06/26(火) 07:00
http://piza.2ch.net/test/read.cgi?bbs=tech&key=991992977&ls=50

13 :頭プアー:2001/06/26(火) 07:01
素数ってなァに?

14 :デフォルトの名無しさん:2001/06/26(火) 07:26
>>13
整数環の素元のこと

15 :デフォルトの名無しさん:2001/06/26(火) 07:37
いや、Delで文字列を数値演算として
扱えるユニットがあるので
それで扱える桁数かな、とオモタわけ。

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

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

16 : 刺身定食:2001/06/26(火) 08:34
素数判定が速いソフトってないかな。

17 :デフォルトの名無しさん:2001/06/26(火) 10:20
だから、、、素数のリストを内部的に保持しておいて
そのリストと照らし合わせる判定方法にすれば
それこそ速い遅いの問題ではないっていってるんだが?

通じないか?

18 :デフォルトの名無しさん:2001/06/26(火) 12:27
>17
リストってお前、500bitとか1000bitとかの素数まで全部保持するつもりか?

19 :デフォルトの名無しさん:2001/06/26(火) 12:49
>>18
>32bit整数MAXまでの素数リストがあれば実用上何も問題なし。
って言ってんだろ?

20 :デフォルトのやる気なし:2001/06/26(火) 13:02
>>11
HNにちょっとだけワラタ。 (^◇^;

21 :デフォルトの名無しさん:2001/06/26(火) 13:03
>>18 の頭は、>>5 に影響されて、どっかに飛んでいきました。

22 :デフォルトの名無しさん:2001/06/26(火) 13:28
技術的に考えたいのなら
メモリ管理とか考える方が先。

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

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

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

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

23 :デフォルトの名無しさん:2001/06/26(火) 13:33
そもそも、1は
Rabinの定理というものが知りたかったのか...

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

24 :デフォルトの名無しさん:2001/06/26(火) 13:52
http://www.hotwired.co.jp/news/news/technology/story/2234.html

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

25 :デフォルトの名無しさん:2001/06/26(火) 13:58
>23

http://www.em.cache.waseda.ac.jp/research/complexity.html
確率性アルゴリズム

補題3
  n≧2 が素数である必要十分条件は 0<a<n なる任意の a に対して
     a^(n-1) ≡ 1 (mod n) が成り立つことである.
補題4
  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 が存在する.

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

26 :デフォルトの名無しさん:2001/06/26(火) 14:18
まず a ≡ b mod n は C言語風に書けば (a % n == b % n) == True の事

素数が注目されているのは暗号化技術に必要だから、
例えば、公開鍵暗号には2つの素数 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さえあ
れば、試してみて復元出来るならそれでいい。

という事で、p,qがホントに素数かどうかより 素数らしい値を手軽に
得られる方法でいいという理屈。

27 :26:2001/06/26(火) 15:38
× e*d=n; n % (p-1)*(q-1) == 1 になるように e,dを見つけて
○ n=p*q; (e*d) % (p-1)*(q-1) == 1

28 :デフォルトの名無しさん:2001/06/26(火) 20:48
だめだ。
Max32ビットの素数リストすら作れない、
メモリの確保が不可能だ。

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

29 :デフォルトの名無しさん:2001/06/26(火) 21:16
>>1 の知りたいこととは全く無関係に、やや面白いスレになりつつありますな

30 :デフォルトの名無しさん:2001/06/26(火) 21:22
9999929/9999931/9999937/9999943/9999971/9999973/9999991/

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

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

32 :31:2001/06/26(火) 21:53
>エラストネス
エラストテネス?恥ずかしい・・・

33 :31:2001/06/26(火) 21:57
て全然駄目ですねごめんなさい。あげないで

34 :デフォルトの名無しさん:2001/06/26(火) 22:19
>>26
中途半端な知識だな。
それはRSAな。
○RSA∈公開鍵暗号
×RSA=公開鍵暗号
×RSA∋公開鍵暗号

35 :54:2001/06/26(火) 22:36
RSAは会社名だよ
何の略か忘れたけど・・・

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

37 :デフォルトの名無しさん:2001/06/26(火) 23:24
30だけど
>>31
つかってるよ。

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

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

38 :荒井注:2001/06/26(火) 23:34
何だバカヤロー
エラトステネスだバカヤロー

39 :デフォルトの名無しさん:2001/06/27(水) 02:11
42億まで計算しようと思ったが1億まで計算するのに40分
かかった(PenMMX266)ので途中で断念した。時間かけ
ればできるけど、それでもまだ10桁だよ。100万桁なんて
想像できん。

40 :デフォルトの名無しさん:2001/06/27(水) 02:23
2chで共同で計算するのもいいかもね。

41 :デフォルトの名無しさん:2001/06/27(水) 03:14
任意の数を与えると、それが素数かどうかを判断する関数を用意すれば?

42 :デフォルトの名無しさん:2001/06/27(水) 03:22
>>41
どういうアルゴリズムを使うわけ?
素数の定義からは逃げられないよ。

43 :デフォルトの名無しさん:2001/06/27(水) 07:28
ルカステストとかラビンの定理とか。

44 :デフォルトの名無しさん:2001/06/27(水) 08:12
>>39
激しく同意。

45 :デフォルトの名無しさん:2001/06/27(水) 08:15
100桁の数一つ扱うのにも
めちゃくちゃややこしいし
一つの数が素数か否かを判断するのにも
膨大な時間がかかる。
この部分を分散計算するという考えそのものが思いつきもしない。
ウツニャ

46 :デフォルトの名無しさん:2001/06/27(水) 09:17
100万桁くらいの+ - * なら ややこしくないだろう。 / % は多少メンドクサイけどね。
 BCDや100進法で 500KByte 掛算時のΣai*biで桁溢れがおきないから大丈夫

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

47 :デフォルトの名無しさん:2001/06/27(水) 09:33
http://www.hotwired.co.jp/news/news/technology/story/2234.html
>素数とは、1、3、5、7、11、13など、
>1およびその数自体でしか割り切れない数

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

48 :デフォルトの名無しさん:2001/06/27(水) 09:33
7時間かかってらあ。
....99999677/99999703/99999721/99999773/99999787/99999821/99999827/99999839/99999847/99999931/99999941/99999959/99999971/99999989/

49 :デフォルトの名無しさん:2001/06/27(水) 09:35
>>47
重箱隅だが2も素数だよ。

50 :デフォルトの名無しさん:2001/06/27(水) 09:42
ギコbasicが多桁対応になってるから、これでなんとか出来ないかい?

51 :デフォルトの名無しさん:2001/06/27(水) 09:53
やるならDelで文字列数値ユニットがある。

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

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

52 :デフォルトの名無しさん:2001/06/27(水) 09:55
気が付いたけど、Windows付属電卓って
すごい桁数扱えるのね。

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

53 :デフォルトの名無しさん:2001/06/27(水) 10:15
1は素数じゃないのでは・・・。

54 :デフォルトの名無しさん:2001/06/27(水) 10:25
1は素数じゃないな。。。

55 :エラトステネス:2001/06/27(水) 10:42
>>4>>22>>32 エラストテネス
>>31 エラストネス
>>37 エラスネス
>>51 エストラネス

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

56 :デフォルトの名無しさん:2001/06/27(水) 10:48
1は素数だろ。
>>53-54どこみていっているんだ?

57 :56:2001/06/27(水) 10:49
ウソだった。氏ぬ。

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

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

59 :デフォルトの名無しさん:2001/06/27(水) 10:56
>>58
1と自分以外で割り切れるかどうか、試せばいいんじゃないの?

60 :デフォルトの名無しさん:2001/06/27(水) 11:02
>59
いや、あの、それじゃ1命令10GHzで多桁除算が出来たたとしてもさ
 10^(10^6)/10^10 = 10^(10^6-10)秒かかるんだよ

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

61 :58=60:2001/06/27(水) 11:11
正確に書くと 2^24 の範囲には素数は 平均7.3個程度ある筈

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

62 :デフォルトの名無しさん:2001/06/27(水) 11:30
100万=1M桁になると、掛算を筆算方式でしたんじゃ 1回の掛算だけで1時間かかるんじゃないの?

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

63 :デフォルトの名無しさん:2001/06/27(水) 11:46
>>60
ハゲしく同意

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

64 :デフォルトの名無しさん:2001/06/27(水) 12:45
Shorのアルゴリズムで素因数分解できるかどうかみてみるってのは?
だれか量子コンピュータを作ってよ。

65 :デフォルトの名無しさん:2001/06/27(水) 12:53
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万分の1になるんじゃないかな?

66 :デフォルトの名無しさん:2001/06/27(水) 13:08
1000000!がいくらになるか想像つきますか?

67 :デフォルトの名無しさん:2001/06/27(水) 13:11
ちなみに、スターリングの公式だと、
n! = n^n * e^(-n) * SQRT(2*pi*n) e^(1/(12*n) - 1 / (360 * n^3) + θ/(1260 * n^5))
(0<θ<1)
です。

68 :65:2001/06/27(水) 13:17
その前によく考えたら、100万迄は a%cはゼロだけど 1000000! % 1000001 はもう予測出来ないな
100万分の1じゃなくて100万分の節約ってだけだな

69 :非65:2001/06/27(水) 13:25
Log[10.0,Gamma[1000000.0]] ≒ 5.57 10^6
Log[2.0,Gamma[1000000.0]] ≒ 1.85 10^7

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

70 :デフォルトの名無しさん:2001/06/27(水) 13:45
でも課題は、
 100万桁よりも長い最初の素数 5万ドル
 10億桁よりも長い最初の素数に対する最高金額の25万ドル
だからなあ

71 :デフォルトの名無しさん:2001/06/27(水) 13:55
100万桁はもう出てると思うよ
http://www.eff.org/awards/coop.html

狙うならもう1桁上だね

72 :デフォルトの名無しさん:2001/06/27(水) 14:10
1000001桁?

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

73 :71じゃないけど:2001/06/27(水) 14:15
>>72
>>71のリンク先を見る限り
>>71の言う一桁上は
1,000,001でなく
10,000,000ってことかと。

74 :デフォルトの名無しさん:2001/06/27(水) 14:25
>>73
なんで桁数が10進法で上昇しなければならないんだ?
意味わからん。

75 :73:2001/06/27(水) 14:39
>>74
俺に言われても(^^;
>>71の言い方が悪い…。

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

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

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

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

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

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


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

78 :デフォルトの名無しさん:2001/06/27(水) 15:47
 その高速判定は、あくまでも確率的な判定方法ですよ
 ソースはこれから辿って下さい
http://search.net-newbie.com/php/function.gmp-prob-prime.html

79 ::2001/06/27(水) 16:28
>>77
 実装が大変なんですか。あまりよくわかっていませんが、
桁が多いと、扱える変数(?)がなくなると言うことですか?
それにともなって、細工が必要になるのかな。素人でスマン。

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

僕が最初、試したアルゴリズムは、
6n±1(n≧1、nは自然数)を、素数の候補として、
手当たり次第に√(6n±1)に達するまで、

(6n±1)mod(既に得られてる素数)=0

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

80 :デフォルトの名無しさん:2001/06/27(水) 16:37
2^23程度なら256M+128Mバイトのメモリを実装すればオンメモリで可能だろう

81 :デフォルトの名無しさん:2001/06/27(水) 16:44
× 2^23
○ 2^31

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

82 ::2001/06/27(水) 16:47
>>77
>>78

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

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

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

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

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

83 ::2001/06/27(水) 16:55
>>80、81
アドバイスありがとう。そういう計算方法も知らなかったんです。理系だけど。ちょっとやってみます。

84 :デフォルトの名無しさん:2001/06/27(水) 16:57
フラグで記録する時は 偶数を記録しないだけで 1/2に さらに3の倍数を記録しないのを追加すれば 2/6
でいいからね

85 :デフォルトの名無しさん:2001/06/27(水) 17:01
>>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という関数が判らないのか?

86 :デフォルトの名無しさん:2001/06/27(水) 17:09
/* 符号なし 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()とかしてるし)。
あと、メモリ確保できなかった時の対処もしてない。

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

87 :デフォルトの名無しさん:2001/06/27(水) 17:10
さらに噛み砕いた説明が欲しいなら

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

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

89 ::2001/06/27(水) 17:18
>>86、87
お忙しい中、ありがとうございます。
力不足ですいません。
時間かかりますけど、やってみます。

90 :デフォルトの名無しさん:2001/06/27(水) 17:57
>>64
Shorのアルゴリズムも、Rabinと同程度に確率的な判定法でしかないっす。

91 :デフォルトの名無しさん:2001/06/27(水) 18:08
ルカステストは?

92 :デフォルトの名無しさん:2001/06/27(水) 21:56
>>90
Shorのアルゴリズムとは

93 :デフォルトの名無しさん:2001/06/28(木) 01:24
「コンピュータと素因子分解」
和田秀夫 / 遊星社

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

94 :デフォルトの名無しさん:2001/06/28(木) 01:35
やっぱアルゴリズムとかコアな本は外人が書いたやつに限るね。
日本の東北大の教授とかが書いたやつなんてコードに無駄なとこ多くて読む気がうせる

95 :デフォルトの名無しさん:2001/06/28(木) 07:01
>>50
これだね
http://www.geocities.co.jp/SiliconValley-Oakland/8522/data/mformula.txt

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

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

96 :デフォルトの名無しさん:2001/06/28(木) 08:08
>>92
おおざっぱに書くと

1.Nを素因数分解したいとする。
----- 因数が見つかるまで2〜5を繰り返す ------
2.Nと互いに素で x<N な x を取ってくる
3.x^r == 1 (mod N) になる最小の自然数 r を計算
4.rが偶数になる確率は1/2程度。
5.rが偶数の時 gcd( N, x^(r/2)±1 ) のどっちかが、1やNでないNの約数を与える確率も1/2程度。
-------------------------------------------

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

97 :デフォルトの名無しさん:2001/06/28(木) 08:44
>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;

98 :デフォルトの名無しさん:2001/06/28(木) 09:37
複雑な事はわからないですが

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

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

ダメ?

99 :デフォルトの名無しさん:2001/06/28(木) 09:40
ガイシュツだっけ?

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

100 :99:2001/06/28(木) 09:45
>>48の方がよっぽど上だった,,,

101 :デフォルトの名無しさん:2001/06/28(木) 09:57
>>98
それより計算量が多かったり低速だったりする方法って思いつかないぞ。(w

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

102 :デフォルトの名無しさん:2001/06/28(木) 10:02
>>101
いや、よくわかってないもので。

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

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

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

103 :デフォルトの名無しさん:2001/06/28(木) 10:08
http://www.valley.ne.jp/~kanai/IT/sources/java/primeTable.html
Javaなのに速いな....

104 :デフォルトの名無しさん:2001/06/28(木) 10:11
>>31に書いてあった。
>基本的にある数Nが素数かどうかはN^1/2までの約数が無いかを調べればよいので
欝...

105 :デフォルトの名無しさん:2001/06/28(木) 12:59
計算量の比較

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

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

106 :デフォルトの名無しさん:2001/06/28(木) 13:09
√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 の点

107 :デフォルトの名無しさん:2001/06/28(木) 14:28
>>105
スマSO、なぜかうちでは
フルイの方が圧倒的に遅い

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

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

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

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


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

108 :デフォルトの名無しさん:2001/06/28(木) 14:32
>>105
ふるいに使うのは足し算だけで十分だと思われ。
今とりあえず組んでみたら1億までで30秒くらい。

109 :デフォルトの名無しさん:2001/06/28(木) 14:44
>>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ってのが入っているんだけど
このやり方がダメなの?

110 :デフォルトの名無しさん:2001/06/28(木) 15:00
>>このやり方がダメなの?
ダメダメダメダメダメ。でもうちの新人に[i...j]のうちのnの倍数を
全て出力するプログラム、って課題を出したら、最初は確かに mod (%)
使ってた...

111 :デフォルトの名無しさん:2001/06/28(木) 15:03
じゃ、教えてくださーい
ある数がある数の倍数かどうかを判断するとこがダメなのか?

112 :デフォルトの名無しさん:2001/06/28(木) 15:11
ある数の倍数かどうか

たとえば10進数なら

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

2進数の場合はどうだろう?
3の倍数 2桁毎に足した結果の下2桁が11である
5の倍数 4桁毎に足した結果の下4桁が1001

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

113 :デフォルトの名無しさん:2001/06/28(木) 15:12
あ、フルイの場合は、単に Nの倍数は Nを足すだけでいいから掛算も不要という意味ね

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

114 :112:2001/06/28(木) 15:23
あ、間違いだ
 n桁毎に足して、その結果がさらにNの倍数かどうか見なくちゃいけない

3の場合 3*17=$33 2進で2桁毎に足したら6だから 0110b さらに2進2桁毎に足したら 11b


115 :デフォルトの名無しさん:2001/06/28(木) 15:32
子ね

116 :デフォルトの名無しさん:2001/06/28(木) 15:34
すんまそ。
ほんとにバカなんで教えてくださいませ

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

117 :デフォルトの名無しさん:2001/06/28(木) 15:47
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あまる。

多分実用性はないな。

118 :117:2001/06/28(木) 15:48
> -1+0-0-1+0=0が3で割り切れるので3の倍数。
訂正 +1-0+0-1+0=0 だ。

119 :112:2001/06/28(木) 15:51
とりあえず 素数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 <- 32より長いからこれも使い難いな
 41 20
 43 14
 47 23

120 :112:2001/06/28(木) 15:57
まあ、2進法で8ビット単位に足して 全部の桁の総和が
3,5、17 の倍数かどうかチェックすれば一挙に3つ計算出来る
から、多少は早くなるんじゃないかな

121 :117:2001/06/28(木) 16:07
112がやろうとしていることが今分かった。
pが素数なら、2^(p-1)-1 が p で割り切れるよ。

122 :デフォルトの名無しさん:2001/06/28(木) 16:14
暗号化とかで計算に必要な素数ってのは、だいたい
どのくらいの桁数なの?

123 :112:2001/06/28(木) 16:22
>>121
あ、ホントだ という事は 29は 28bitで循環か、なんで抜けたんだろ 浮動小数点使ったせいか

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

124 :デフォルトの名無しさん:2001/06/28(木) 16:45
    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

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

125 :112:2001/06/28(木) 17:07
10進法で記録した場合
  2  最下位桁のみ
  3  1桁単位に加算
  5  最下位桁のみ
  7  6桁単位に加算
  11  2桁単位に加算
  13  6桁単位に加算
  37  3桁単位に加算
  41  5桁単位に加算
  73  8桁単位に加算
 101  4桁単位に加算
 137  8桁単位に加算

126 :112:2001/06/28(木) 17:42
>>124
J=I*2 にしてるけど、Jより小さい素数の倍数は不要でしょう?
説明すると、
まずは2刻みに4以上を全部消すよね?
次は3刻みなんだけど6は既に消えてるし2の倍数は消えてるから9から消す
その次は5刻みだけど
やっぱり5*2 も5*3も消えてるから5*5から10刻みに消せばいい

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

127 :112:2001/06/28(木) 17:53
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 がループ中再評価されるんなら外に出してね

128 :112:2001/06/28(木) 19:13
メモリーが256メガくらいあれば動くみたい
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;
}

129 :112:2001/06/28(木) 19:40
最後の何個か
2147483477
2147483489
2147483497
2147483543
2147483549
2147483563
2147483579
2147483587
2147483629
2147483647

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

130 :112:2001/06/28(木) 19:49
>107のいう 10000000 迄ならハードディスクに書き出しをすれば3秒で終わった
 画面に出したらスクロールに3分かかった

131 :デフォルトの名無しさん:2001/06/28(木) 19:51
>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

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

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

132 :デフォルトの名無しさん:2001/06/28(木) 20:47
http://www.fortunecity.com/emachines/e11/86/primesec.html
???

133 :112:2001/06/28(木) 21:53
http://photo-m.tp.chiba-u.ac.jp/~adeno/sci/math.htm

>それは19変数の多項式であり、19個の自然数をその
 多項式fに代入してfの値が正であればそれが素数であり、
 逆にすべての素数についてもp=f(x_1, x_2, ..., x_19)を満たす自然数 x_1, x_2, ..., x_19が存在する。

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

134 :デフォルトの名無しさん:2001/06/28(木) 22:50
素数の判定をそんなに急いでどうすんの?

135 :デフォルトの名無しさん:2001/06/28(木) 23:06
>>91 ルカステストってメルセンヌ素数かどうかのテストじゃなかったっけ?
100万桁の最初の素数を見つけるような場合には意味が無いとおもわれ

136 :デフォルトの名無しさん:2001/06/28(木) 23:12
全てが素数になる方程式なら ウィルソンの定理があるんじゃないの?
(N−1)!+1 ≡ 0 mod N

if (((N−1)!+1 )% N==0) 素数です else 素数ではありません

137 :デフォルトの名無しさん:2001/06/29(金) 00:06
>>133のいうオイラーの式って、その数列の上に全ての素数が
存在する訳じゃないし、それを満たせば全部素数って事もないよねぇ。
ていうか当然か。。。

138 :デフォルトの名無しさん:2001/06/29(金) 00:08
>>133の言ってる素数を生成する多項式っていうのは聞いたことがあるけれど、
代入しても負の値が出ることがかなり多いそうだよ。(つまり素数じゃないやつネ)

139 :デフォルトの名無しさん:2001/06/29(金) 01:16
>>130
話それるけど、
そのようなプログラムでコンソールに連続的に表示するとき、
Ctrlを押してるとスクロール(表示速度)が速くなるんだけど
理由知っている人いる?

140 :デフォルトの名無しさん:2001/06/29(金) 07:16
>>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)

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

141 :デフォルトの名無しさん:2001/06/29(金) 08:49
>>68
>1000000! % 1000001 はもう予測出来ないな
 136が書いているウイルソンの定理から 1000001 が素数なら
  1000000! % 1000001 =1000000 です
 しかし、1000001=101*9901なので階乗は割切れますから結果はゼロです
1000003が素数なので
1000002! % 1000003 =1000002
1000001! % 1000003 =1
1000000! % 1000003 =500001
と予測出来ます

142 :デフォルトの名無しさん:2001/06/29(金) 09:41
>>126-127=112さんありがとう。
ついていくのに必至だ。

>>48で7時間って言ってたけど
それが5分になった(w
うれしい!

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

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

144 :デフォルトの名無しさん:2001/06/29(金) 09:52
>>143 そのルーチンじゃ、掛算は桁数の2乗に比例するでしょ?
1Gマシンで100万桁の掛算だけで既に1日で出来ないと思うよ

145 :デフォルトの名無しさん:2001/06/29(金) 10:06
>>144
そっか。やはりWiredNEWSで言われているように
分散処理で考えないとだめか。
100万桁といわず、ここでどんどん桁数を上げていきません?

146 :144:2001/06/29(金) 10:11
いや、もっとずっと小さい桁数でダメになりそうだし
そういうふうに動的に取るとメモリが細切れになって大きな配列が最後は取れなくなる


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

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

147 :144:2001/06/29(金) 10:20
ベースにするなら >>95 のギコbasicベースの方がいいんじゃないの?
使い方は文字列で与えて文字列で出て来るから殆ど同じだし
入力で "2^1000+1" みたいな書き方も出来るし

内部100進だから10進に比べてそれだけで掛算が4倍早くなるし

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

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

148 :デフォルトの名無しさん:2001/06/29(金) 10:26
そっか。
unit mformula;
も数値を文字列で持っているんだね。

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

149 :144:2001/06/29(金) 10:30
それと、目標を作らないとチューニング出来ないと思うよ
 例えば1万桁台の最初の素数を探してみたいとか

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

150 :144:2001/06/29(金) 10:40
試したときは
function N100Str(str:String):String;
begin
  Result:=FormulaNumStr(TN100Num,str');
end;

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

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

151 :デフォルトの名無しさん:2001/06/29(金) 15:23
あれ?さがってる。

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

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

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

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

152 :144:2001/06/29(金) 18:26
64bitまで計算させたい時のことを考えたてるんだけど
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) % と分解出来るかな?

153 :デフォルトの名無しさん:2001/06/29(金) 20:23
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;

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

154 :学生:2001/06/29(金) 21:26
>>153
「掛算のFFT化」とは?なんですか?
FFT?

155 :デフォルトの名無しさん:2001/06/29(金) 21:28
がんばってください>>153

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

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

156 :デフォルトの名無しさん:2001/06/29(金) 21:38
>>154
FFT=高速フーリエ変換

157 :学生 :2001/06/29(金) 21:58

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

158 :デフォルトの名無しさん:2001/06/29(金) 22:10
>>157
少しはさかのぼって読んでみたらどうか。>>147を読め。

159 :153:2001/06/29(金) 22:12
うっ 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;

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

161 :153:2001/06/29(金) 22:36
 if n>= ends then break; を
 if int64(n)*int64(n)>ends then break; にしたら

1234567890123481 1234567890123493 1234567890123527 1234567890123583
1234567890123647 1234567890123661 1234567890123671 1234567890123697

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

162 :デフォルトの名無しさん:2001/06/29(金) 23:44
スタンドアロンで計算するのもよいけれど、
分散システムのモデルを考えてみるのもよいかもしれない
今日このごろ。

163 :デフォルトの名無しさん:2001/06/29(金) 23:46
>>162
それって強力にむずかしいんだよ。

164 :デフォルトの名無しさん:2001/06/29(金) 23:52
>>163
難しいっていうのは、要するに小分けにされた問題だけ持って帰って、
計算結果を返してくれないクライアントの事とかを考慮していくと、
強力に難しくなるって事なんでしょうか?

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

165 :デフォルトの名無しさん:2001/06/30(土) 01:15
結局、どういう理論があるんだ?
誰かまとめてくれ。

あと%って何の演算子?

166 :デフォルトの名無しさん:2001/06/30(土) 01:18
モジュロ演算子でしょ。
つまり剰余を求めるやつ。
VBだと、"Mod"。Rexxだと"/"。

167 :デフォルトの名無しさん:2001/06/30(土) 08:18
じゃ、メルセンヌ素数について
2^n-1 が素数ならメルセンヌ素数
2^2-1=11b=3 2^3-1=111b=7

>>119 が書いてるように 3の倍数かどうかは2桁毎に区切って加算 7は3桁毎に加算
して判断してもいいから

2^4-1=1111b は3+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

168 :デフォルトの名無しさん:2001/06/30(土) 08:32
Lucasテストに必要なのは  (a^2-2)%b という計算だけ

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

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

169 :デフォルトの名無しさん:2001/06/30(土) 08:37
モジュラーな掛算だけど、 2^nのモジュラーは上位桁を捨てるだけで良い

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

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

170 :デフォルトの名無しさん:2001/06/30(土) 09:19
でもこれ、どうやって分散したんだろ?
まさか掛算?
真中の桁で分割して (a+b)^2=a^2+2ab+b^2で3人で計算して持ち寄るとか?

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

171 :デフォルトの名無しさん:2001/06/30(土) 14:06
離散フーリエ変換して各級数を持ち寄る、とか?

172 :デフォルトの名無しさん:2001/06/30(土) 15:17
おっかしいなあ。TBitsで2^31個のビットを確保したら
それだけで落ちるなあ。

173 :デフォルトの名無しさん:2001/06/30(土) 16:11
>>172
そらあんた、2^28 byte の連続したメモリブロックが必要になったんでしょ?
主記憶足りてる?
最初に SetSize で最大値を設定したかい?
俺んとこではうまく動くよ。

174 :デフォルトの名無しさん:2001/06/30(土) 17:16
主記憶は
どのくらいだと足るの?

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

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

175 :デフォルトの名無しさん:2001/06/30(土) 18:16
function LucasTest(n:int64):boolean;
var dt,odt:array of Word;
var cy:int64;
var i,hDt,hPg:integer;
var mask:LongWord;

 procedure mul2; //2乗する処理
 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;

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

176 :デフォルトの名無しさん:2001/06/30(土) 18:17
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; //初期値を4に
rep:=n;
while rep>0 do begin
 odt:=copy(dt,0,hdt+1);
 mul2; //2乗する
 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;//全部が1は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;

177 :175:2001/06/30(土) 18:22
ほんの4桁の素数 n を入れたらもういつまでたっても計算が終わりません

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

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

178 :デフォルトの名無しさん:2001/06/30(土) 18:35
>174 計算中に落ちる時に何か報告ありませんか?
  確保してない所に書いてたり読んでたり?

179 :175:2001/06/30(土) 19:17
いま、 2^n-1 の n=2281付近で、1つの値の検査に1秒程度かかっています。
 今は掛算を筆算方式でするから n^3 のオーダで計算時間が増えます

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

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

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

180 :デフォルトの名無しさん:2001/06/30(土) 21:13
>>178
確保していないところに書き込んでも、その位置に自動的に size を
拡大してくれるから、落ちることはないんだけどなぁ。
逆に、少しずつ大きなところをアクセスするようなコードだと、
毎回巨大なメモリコピーが起きるから遅く見えるかもね。

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

181 :デフォルトの名無しさん:2001/06/30(土) 23:02
うーん、確保していないところに書き込んだりするような
ことは全くしていないのですが...
遅いとか遅くないとかではなくW2kがブルーバックで落ちるからな(w


  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;

182 :175:2001/06/30(土) 23:30

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

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

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

183 :175:2001/07/01(日) 08:36
Lucasテストを2分するアイデアその1

  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が1なら(b-a)と置き換えて 2^((n+1)/2)と比べればいい
 さらに a^2 %b= 2^((n+1)/2)+2 となるaを求めて同じようにする

2台のパソコンで片方は a'=a^2-2 を 片方はa'=√(a+2)を
それぞれ n/2回求めて結果を比較すればいい

184 :175:2001/07/01(日) 09:19
√(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)jは aが奇数の間は2乗するより高速

185 :175:2001/07/01(日) 09:25
×a'=√(a+2)jは aが奇数の間は2乗するより高速
○a'=√(a+2)jは aが偶数の間は2乗するより高速

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

186 :175:2001/07/02(月) 07:41
× a  = c *2^((n+1)/2)
○ a  = ±√c *2^((n+1)/2)

187 :デフォルトの名無しさん:2001/07/02(月) 10:29
結論としては、 2^18=262144 のFFTを1/6秒で出来れば1日に1つの検査が出来る=世界記録が狙えるわけね

この数字は無理筋?

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

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

188 :デフォルトの名無しさん:2001/07/02(月) 12:52
31bitMax=約20億を計算しようとしたら
土日通して計算しても10億までしかやってなかった
ウツだ。

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

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

189 :デフォルトの名無しさん:2001/07/02(月) 15:48
>>188 それはまさか1個単位に全部の素数で割る方法じゃないのか?
  フルイ方式なら
>>128 なら30分程で終わるぞ。(ただしHDに書いた場合)
>>159のコードでも40億〜プラス1万の素数とか計算したら数秒だ

メルセンヌ素数とルカステストについては
>>167-168
さらに
>>175-176 がルカステストのサンプルも書いてるじゃないか?

まとめると、

1) 2^n-1 つまり、 2進で1が続く数でかつ素数ならメルセンヌ素数と呼ぶ

2) nが素数でなければ2^n^-1も素数では無い

3)ルカステストは メルセンヌ素数を見つけるのに使う。素数nに対して
   a=0;を初期値にして a:=(a^2-2)mod(2^n-1) をn度繰り返して途中でゼロになれば
  2^n-1はメルセンヌ素数

4)a mod 2^n-1 は aを2進で表示して 2^n 以上の桁を 下位に回して足せばいい

5) -a mod b = b-a mod b

6) a^2 mod b= (b-a)^2 mod b

7)a*b は それぞれの桁数をan,bnとしたらan+bn-1桁以上の複素フーリエ変換して
 求めた各要素を2乗した後、逆フーリエ変換すればいい

190 :デフォルトの名無しさん:2001/07/02(月) 16:37
>>189
そうかーー、HDに書けばよかったか。
なんでもメモリ上にて計算するしか頭になかった。
それでメモリ壊れてOS落ちて悲しくなっていたのだ

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

粗忽者です。

191 :デフォルトの名無しさん:2001/07/03(火) 08:34
>>187 検算したら違ってるよ

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

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

192 :デフォルトの名無しさん:2001/07/03(火) 11:13
>>189
>2^n^-1
これが顔文字に見えた俺は逝ってよしですか?

193 :デフォルトの名無しさん:2001/07/03(火) 11:19
2^n-1が1千万桁になる素数nは
33219281
33219283
33219313
33219317
33219323
33219341
33219359
33219367
33219379 こんな調子であるぞ。

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

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

194 :デフォルトの名無しさん:2001/07/03(火) 11:33
サーバは、こんな感じだろうな

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

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

195 :デフォルトの名無しさん:2001/07/04(水) 10:22
age

196 :デフォルトの名無しさん:2001/07/07(土) 19:00
メルセンヌ素数について勉強中。

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

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

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

197 :デフォルトの名無しさん:2001/07/07(土) 19:43
>>196
7 = 2^3 -1 は (3-1)番目の項 14 を割り切るので素数。
15 = 2^4-1 は (4-1)番目の項 194 を割り切らないので素数でない。
31 = 2^5-1 は (5-1)番目の項 37634 を割り切るので素数。

という判定法。

198 :デフォルトの名無しさん:2001/07/07(土) 20:54
わかったーー。

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

ありがとう〜♥

199 :デフォルトの名無しさん:2001/07/08(日) 12:18
>>175-176 のコードの mul2をFFT化してみよう
 ttp://homepage2.nifty.com/m_kamada/docs/fftmul.htm
ここを参考にして考えてみる。

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

200 :199:2001/07/08(日) 12:33
>>191 の言を信じて 2^21のFFTを用意する事にする
FFTは基数4の方が若干高速だけどメンドクサイから基数2で我慢しよう

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

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

201 :199:2001/07/08(日) 13:03
>>199 のサイトではA*B を計算させる為に2組の実数列の同時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の範囲に収まる。
しかし、最後の2乗はどうだろう?
 (a+bi)^2=a^2-b^2 + 2*a*bi の計算で aが最大桁、 bが1桁大
なんて事になった場合a^2-b^2は精度を失う
これは問題にならないのだろうか?

202 :デフォルトの名無しさん:2001/07/08(日) 14:50
>>201
そこだけ倍精度つかいな。

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

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

203 :199:2001/07/08(日) 16:53
どうやら、半分のサイズのFFTはそれほど効率的ではないようだ
せいぜいまともにやった時の7割にしかならないらしい。
そこで、2 組の実数列の同時FFTを使ってそれぞれの2乗を求める
事を考える。
実数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]))
と求められる。
つまりこの2式を加算して逆FFTすれば 実部と虚部に 2乗x4倍値として求まる
(どうせFFT->逆FFTの後 2^n=Nで割る必要があるので4Nで割ればいい)

204 :199:2001/07/08(日) 17:06
>>202 倍精度っていってもlong doubleで doubleより16bit増えるだけだから
・・・・でもやっぱり16bitの2乗は32bitだし
N桁の2乗は中間桁ではN倍されるから32+21=53bit必要だから
やっぱりlong doubleでも ダメっぽいね

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

205 :199:2001/07/08(日) 18:20
で、2^K-1のモジュラーの2乗を どう計算するかだけど
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 の方法で 2実数のFFTを使い
周波数軸上で X^2,Y^2, X*Yを求める
から、2*X^2+Y^2を 周波数軸で加算し実数部で得
られるようにし X*Yを虚部に得られるように調整
した後逆FFTする
後は正規化して加算すればいい。
ただし、FFTする時は x,yは8bit単位に分割し
計算精度は long doubleでやらなければならない

206 :199:2001/07/08(日) 18:46
double で計算するのと Extended(long double)で計算するのって
やっぱり速度上違いが出るのだろうか?
doubleが早いなら FFTはdoubleで計算して 2乗後は
Extendedで計算ればいいんだけど

207 :199:2001/07/08(日) 21:30
とりあえず今日 出来た所まで

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

208 :199:2001/07/08(日) 21:32
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;
// 最初の1回目
  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]; { 上位は必ずゼロなので この2行不要 }
   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;

209 :199:2001/07/08(日) 21:36
//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);
 // 最初の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;

210 :199:2001/07/08(日) 22:07
試験用に書いた 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;

211 :デフォルトの名無しさん:2001/07/09(月) 08:26
で? 結局 ルカステストにはどう組み込んだらいいのだい? あげ

212 :デフォルトの名無しさん:2001/07/09(月) 09:15
ルカステスト 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

213 :デフォルトの名無しさん:2001/07/10(火) 08:19
細かい訂正 >>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

214 :199:2001/07/11(水) 09:36
コードはとりあえず全部書いたけど、バグだらけでデバック中

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

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

215 :デフォルトの名無しさん:2001/07/11(水) 09:47
いまさらなんですが

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

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

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

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


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

217 :199:2001/07/11(水) 10:19
//使ってる変数
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;

218 :199:2001/07/11(水) 10:24
正規化 ---これもどうも怪しい
 FFTによる掛算結果をそのまま再度2乗するとデータ精度が保て
ないので、各桁を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乗は正常に働く筈)

219 :199:2001/07/11(水) 11:00
{モジュラーな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;

220 :デフォルトの名無しさん:2001/07/11(水) 11:49
>>194
検査すべき素数を提供する時にパスワードを生成して渡しておくんじゃダメ?
そのまま検査された値全部にフラグとパスワードを保存するとディスクの無駄になるけど。

221 :デフォルトの名無しさん:2001/07/11(水) 12:25
>>215です。
>>216さん

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

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

222 :デフォルトの名無しさん:2001/07/11(水) 12:42
var tbl:array of byte; とすると

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

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

223 :199:2001/07/11(水) 15:32
まだデバックが終わってないから判断には早いかもしれないけど

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

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

224 :199:2001/07/11(水) 21:04
× n:=iLog2(m);
○ n:=iLog2((m+7)div 8);
という事で 上の数字は1/8になった
k=2281付近で2秒
k=4423付近で8.5秒
k=9689付近で42秒


他にも細かいミスを修正した結果、とりあえずこの付近では
ttp://www.geocities.co.jp/Technopolis/1793/mersennu.html
の表の通りの判定が出来るようになった

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

225 :199:2001/07/11(水) 22:13
http://www.geocities.co.jp/SiliconValley-Oakland/8522/data/LucasTest.txt

226 :199:2001/07/11(水) 22:40
プログラムが正しければ世界記録を狙った場合、
2ヶ月程で1つの数の判定が出来るだろう。

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

227 :デフォルトの名無しさん:2001/07/11(水) 22:48
>>226
100堕慰くらいで実行スレばいいよ。
そんなに2chプログラム人口いないか?

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

229 :199:2001/07/11(水) 23:00
>227
1bit単位に分割してFFTで掛算すれば何ビットの精度が必要だろう?

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

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

230 :デフォルトの名無しさん:2001/07/12(木) 03:03
遅くなりましたけど>>222
ありがとー。

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

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

231 :デフォルトの名無しさん:2001/07/12(木) 07:42
>>230
ビットを落とすのは not して andでもいいし
すでに立っているのが判ってれば xorでもいいよ

232 :デフォルトの名無しさん:2001/07/12(木) 09:39
めった、簡単なこと質問してスマソですわ。

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

つうことですか?

233 :デフォルトの名無しさん:2001/07/12(木) 13:10
他にもリセットの方法はあるよ

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

234 :199:2001/07/12(木) 19:05
実際に浮動小数点で何bit あれば判らないので
Single = 32bit 仮数部23bitにして計算させてみた。
Doubleの結果と比較して k=251の時に差が出た

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

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

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

235 :デフォルトの名無しさん:2001/07/12(木) 19:39
毎度ありがとうございます。

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 );
すると大丈夫みたいっす。

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

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

236 :199:2001/07/12(木) 20:14
こうやった方がいいかも

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

さらに

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

237 :デフォルトの名無しさん:2001/07/12(木) 20:20
>>236タンさっすが!

238 :199:2001/07/12(木) 21:16
これなら世界記録付近の計算は1件2週間〜1ヶ月程度で出来そうな感じだけど
>>71の賞金を貰うならさらに 100倍かかるな
100倍になると、今のままだと8年はひとつの数の判定にかかる計算
それにメモリもオンメモリで出来る限界付近だからもっと速度は落ちるかも

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

239 :デフォルトの名無しさん:2001/07/12(木) 23:56
にしても、ちかずいておりますが...
並のパソでそれだけいくと
もちっと速いコンピつかえば&
数年でクロック数上がれば、どんどん賞金に知被いて....

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

240 :デフォルトの名無しさん:2001/07/13(金) 08:54
>数年でクロック数上がれば、どんどん賞金に知被いて....

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

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

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

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

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

241 :デフォルトの名無しさん:2001/07/13(金) 09:07
へたに固定少数点化して命令数を増やすよりは
128bitSIMDにアセンブラレベルで対応した方が速いんじゃないかな
1クロックで4並列演算とか出来るのは、あれは単精度(32bit)の場合だっけ?

242 :デフォルトの名無しさん:2001/07/13(金) 09:53
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

243 :199:2001/07/13(金) 11:22
今、k=11213 付近判定させると 15秒で結果が出る
1回の計算 は 15/12112秒 1秒に750回計算出来ている。

1回の計算は FFT+2乗+逆FFT+2分割正規化処理
fftのサイズは 2048 = 11ステージ でFFT+逆FFTをやってるから
1回のFFTの複素積和は 750*2048*11*2=33.8M/秒
複素積和は4回の掛算と2回の加算だから
今現在、200MFLOPS 付近が出てる模様
メモリサイズは 2048*4*2=16Kbyteだからデータキャッシュはヒットしてる筈
この数字はathron800Mでの話

たしかに倍精度で2GFlopsとか出れば一挙に10倍になるんだから
P4は魅力的だな・・・でもまだ買い換え予定はないし

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

244 :199:2001/07/13(金) 13:20
今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の順にネスト
させたい。どうループを書けばいいんだろ?
再帰を使わないと書けないかな?

245 :199:2001/07/13(金) 14:51
このページを見つけました。
ttp://www.pluto.ai.kyutech.ac.jp/plt/matumoto/fftb_export.c.txt

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

246 :199:2001/07/15(日) 08:55
結局、こういう感じです。 再帰を使わなくてもいいけど、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); //2分割しながら深くネストして
   ifft0(nc-1,st+stp ,stph);
   ifft1(nc , st   ); //一番奥から順に実行する
   end;
 end;
begin
 ifft0(n,0,1 shl(n-1)); //呼び出しの根っこ
end;
fftの方は逆に1ブロック実行してはネストしてゆく格好で書けました

247 :199:2001/07/15(日) 10:24

整数化する為に、
1)sin,cosの計算精度を向上させる為順計算でなく
  1/2に分解して再帰する手順に要変更
しかし、問題はm=2^16を超える頃にcosはほぼ1になる場
合がある点。
これを解決するには、
2)cosではなくcc=1-cosを記録して計算時には
   x*cos=x*(1-cc)=x-x*ccとする
3)基数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))
1回のステージ毎に 少数点位置は丸め後右に動かす
((a*b-(b&(a>>32))<<32)+$7fffffff)>>31

FFT終了時には 少数点はH-Nにある
2乗計算後に2(H-N)-Nだけシフトし、少数点位置をNにする
逆FFTも1回のステージ毎に 少数点位置は丸め後右に動かす
((a*b-(b&(a>>32))<<32)+$7fffffff)>>31
逆FFT終了時には整数化されている

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

248 :199:2001/07/15(日) 10:26
×$7fffffff
○$3fffffff

249 :デフォルトの名無しさん:2001/07/15(日) 11:15
整数でやるんだったら、有限体に値を持つ FFT なんてどうでしょう。
それとも実数でやるのがミソなんでしょうか。

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

250 :199:2001/07/15(日) 13:51
ここにある方法ですね >> 249
ttp://www.geocities.co.jp/WallStreet/6351/mersenne.html

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

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

251 :199:2001/07/15(日) 16:45
>>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
と非常に簡単に計算出来ますよね?

252 :249:2001/07/16(月) 00:56
>>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 で適当に作ってみたらものすごく遅かった。

253 :199:2001/07/16(月) 09:26
賞金を狙う事を考えると 30,000,000 bit以上に目標が必要です。
60*60*24*365=31,536,000なので、およそ1回の2乗にかかる秒数が1つのルカス
テストにかかる年数になります。

たとえば 2^16384 -1 を法にして 16384乗根は 2なので
FFTにした時の各段階の係数は全て 2のベキになります
(8170bit単位に 8192個の単位に分割すれば =66Mbit迄可)
2のベキの係数なので回転操作と加算操作で計算出来ます
(掛算ではないので何ビットあっても不利ではない)

結果、16384bitのの2乗の計算が16384個必要になるけど、
これも 2^256-1の法で 64bit単位に256個に分解すれば
256bitの掛算は(32bit掛算*8^2=)64になります。
掛算回数はわずか 64*256*16384 = 268M

ただしメモリを256Mbyte使っているので、FFTにかかるオーダは
最初の16384乗根で 14 つぎの256乗根で8 で合計14+8=22
一回のステージは2回の加算と1回のシフト、逆変換を含めて
256M/4*22*3*2 = およそ 8.5G という事で
8ギガmipsのマシンでもルカステスト1回にやはり年単位は必要
という事になってしまいます

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

254 :199:2001/07/16(月) 09:46
>>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=1だから、 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の場合よりもさらに簡単です(負数を考えなくていいから)

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

255 :249:2001/07/16(月) 21:54
>>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 となってます。

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

256 :199:2001/07/16(月) 23:16
ホントですね。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|
を2つ掛ければ出来るけど、どっちも逆変換しようとすると 3が出て来て
3の逆数が無いから解けないですね

やっぱりダメか・・・・

257 :199:2001/07/17(火) 00:00
あ、でも、それぞれ
|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|

とすると3倍、6倍の結果になってトータル 18%15 = 3倍になるけど、いちおう計算出来ますね

2^8-1の場合も 15倍になる解なら同じ規則で見つかりました

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

258 :199:2001/07/17(火) 00:14
見つけた規則だと 2^(2^n)-1で 2^nのFFTをすると
2^(n/2)-1 倍された逆FFTならできるという事になります

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

259 :199:2001/07/17(火) 09:02
あ、そうか 2^(2^n)-1は
11*101 =1111
1111 * 10001=11111111
11111111  * 100000001= て調子で 2^m+1 との掛算なんだ

260 :デフォルトの名無しさん:2001/07/18(水) 10:22
じゃ、並列処理をどうやるかって方向でひとつ

261 :デフォルトの名無しさん:2001/07/19(木) 15:37
High(Integer)までの素数列を
常識的な速度で出す事やOnメモリでダス方法はわかったのだけど

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

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

262 :デフォルトの名無しさん:2001/07/19(木) 20:08
integerのメモリーイメージのまま出力したいなら
ファイルマッピングを使ったらどう?

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

263 :デフォルトの名無しさん:2001/07/19(木) 21:47
篩データをこんなのにして持っています
 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
に代入してから保存した方が使いやすそう
(要はファイルサイズの問題)
かと思うのですが

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

264 :デフォルトの名無しさん:2001/07/19(木) 23:43
つーかこんなとこでうだうだやってるより
さっさと分散&ブルートフォースで解決

265 :デフォルトの名無しさん:2001/07/20(金) 07:31
ああ、まだ初心者だったんだね。>>263

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

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

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

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

266 :デフォルトの名無しさん:2001/07/20(金) 14:16
>ああ、まだ初心者だったんだね。>>263

 グサッ  Λ_Λ
     :(゚Д。 ;)
――――(”∴  ) グヒャッ・・・
     .:|| |
    ∴ (__(__)

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

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

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

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

267 :デフォルトの名無しさん:2001/07/20(金) 17:32
> 2,3,5,7,11,13,....という配列

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

268 :デフォルトの名無しさん:2001/07/20(金) 18:07
うーん Delphi3くらい迄は動的配列なんてなかったから

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

269 :デフォルトの名無しさん:2001/07/20(金) 18:23
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
だからビットマップで記録すれば  6倍は効率が良い
2^32/log(2^32)*4はおよそ 1.8GByte

32bitの素数が欲しい場合、16bitの篩を用意しておけばいいのでは?

270 :デフォルトの名無しさん:2001/07/20(金) 23:39
>>267
サンキュウ。
そうか、素数定理というのが抜けてたのか。
その定理よくわからないけど、検索してみる。

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

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

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

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

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

271 :269:2001/07/21(土) 07:20
ごめん log(2)はおよそ 0.69 0.3は log10(2)の場合だった

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

272 :デフォルトの名無しさん:2001/07/31(火) 02:07
age?

273 :デフォルトの名無しさん:2001/07/31(火) 05:27
素数定理ならここなんてどうよ。
http://www.com.mie-u.ac.jp/~kanie/tosm/humanind/yosou.htm#prime_number
数学は科学の女王で、数論は数学の女王だ。
頑張れ。

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

274 :デフォルトの名無しさん:2001/08/02(木) 12:55
夏休み中に賞金狙うやついる?

275 :デフォルトの名無しさん:2001/08/02(木) 13:37
>>274
無理にきまってんだろ。

276 :デフォルトの名無しさん:2001/08/02(木) 18:14
100万人くらい参加したら可能かも・・でもどやって賞金分ける?

277 :デフォルトの名無しさん:2001/08/03(金) 06:09
1千万円くらいあるんでしょ?
百人くらいを一つのグループにして、そのグループで一つの数のチェックを
担当して、 ヒットしたらそのグループ+運営側で半々に わけるってのは?
 百人なら半分でも5万円になる

278 :デフォルトの名無しさん:2001/08/04(土) 13:23
100人いても98人は厨房と思われ

279 :デフォルトの名無しさん:2001/08/04(土) 19:13
プログラム作った人が総取りすればいいだろ。

280 :デフォルトの名無しさん:2001/08/05(日) 15:43
>>279 それじゃ誰が何の為に協力してくれるんだ?
感動的ストーリーでも捏造するか?

281 :デフォルトの名無しさん:2001/08/05(日) 16:01
じゃあ、白血病やSETIなどで協力している奴はなんなの?

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

282 :デフォルトの名無しさん:2001/08/05(日) 16:38
俺のP-133も役に立てるかな?

283 :デフォルトの名無しさん:2001/08/05(日) 17:44
>>282
たたん。最低でもP-500以上。

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

285 :デフォルトの名無しさん:2001/08/08(水) 17:41
Internet Mersenne Prime Search
http://mersenne.org/PrimeNet/

286 :デフォルトの名無しさん:2001/08/10(金) 01:30
素数が無限に存在することの証明として
>もし素数の数が有限個だとして、それを小さい順に p1、p2、p3、・・・、pmとします。
>新たに、(p1・p2・p3・・・pm+1)なる数を考えます。
>この数は、任意の素数pnで割れませんので、素数です。
(一部抜粋)

もし、素数の数が有限個だとして、それを小さい順に 2,3とし
新たに、2X3+1なる数を考える。この数 7 は素数。

もし、素数の数が有限個だとして、それを小さい順に 2,3、7とし
新たに、2X3X7+1なる数を考える。この数 43 は素数。

もし、素数の数が有限個だとして、それを小さい順に 2,3、7、43とし
新たに、2X3X7X43+1なる数を考える。この数は素数。

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

287 :デフォルトの名無しさん:2001/08/10(金) 01:33
>>286
これって… マジで
できそうでないかい?

288 :デフォルトの名無しさん:2001/08/10(金) 01:46
2*3*7*43+1=1807 1807/13=139 (;゚Д゚)

289 : ◆iK1ryBHs:2001/08/10(金) 02:03
>>287
カコワルイ

290 :デフォルトの名無しさん:2001/08/11(土) 20:02
> もし、素数の数が有限個だとして、それを小さい順に 2,3とし
> 新たに、2X3+1なる数を考える。この数 7 は素数。

> もし、素数の数が有限個だとして、それを小さい順に 2,3、7とし
> 新たに、2X3X7+1なる数を考える。この数 43 は素数。

ここがおかしい?

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

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

291 :290:2001/08/11(土) 20:07
> だから「この7は素数」と分かったときに
> 3から7までの数で素数を探す必要があるのでは
これが大変すぎるから、このスレがあるんだよね…
ageてゴメンヨ…

292 :290:2001/08/11(土) 20:18
(* 2 3 5 7)
210

211 は素数

(* 2 3 5 7 11 ... 211)
1645783550795210387735581011435590727981167322669649249414629852197255934130751870910

1645783550795210387735581011435590727981167322669649249414629852197255934130751870911は素数

293 :デフォルトの名無しさん:2001/08/11(土) 21:01
30分走らせたけど、>>292が素数かどうか判定出来んかった・・・

294 :デフォルトの名無しさん:2001/08/11(土) 22:17
> もし素数の数が有限個だとして、それを小さい順に p1、p2、p3、・・・、pmとします。
> 新たに、(p1・p2・p3・・・pm+1)なる数を考えます。
> この数は、任意の素数pnで割れませんので、素数です。

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

295 :デフォルトの名無しさん:2001/08/12(日) 00:59
>292
ちょっとしたことをオンラインで検証できるツールがあるぞー。

■桁数無制限電卓

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

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

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

296 :295:2001/08/12(日) 01:00
ごめんsageてもうた

297 :デフォルトの名無しさん:2001/08/12(日) 09:56
 >>295 のアプレットに 292の数字を入れたら is composite だそうだ

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

298 :デフォルトの名無しさん:2001/08/12(日) 11:03
>>297
うむ。
 1.ふるいで素数を小さい方から1000個生成
 2.とりあえずその1000個で割ってみる
 3.>>25で言う"任意のa"として小さい方の素数50個を使ってRabinテスト
…というアルゴリズムなように読める。

299 :デフォルトの名無しさん:2001/08/15(水) 23:24
>>294
背理法って知ってますか。
これは素数の数が有限個だとして矛盾を導いているんです。

300 :デフォルトの名無しさん:2001/08/16(木) 00:18
>>299
命題の対偶という概念が分からないと
背理法で何が導き出されるか分からないのでは。

301 :294:2001/08/16(木) 08:13
>>299
んだから、素数が無限個であることの証明としては問題ないのはわかっちょる。
だからといって、素数を無限に生成する式としては使えるとは限らないだろ?

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

302 :299:2001/08/16(木) 09:39
>>294

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

303 :デフォルトの名無しさん:01/09/24 09:49
保存アゲ

304 :デフォルトの名無しさん:01/10/08 20:55
あげ

305 :デフォルトの名無しさん:01/10/08 21:16
>>303-304
2ヶ月もなにもないって事は、
終ったってことさ。
それにdat落ちしても、このスレのリンク覚えておけば参照できるよ。

306 :デフォルトの名無しさん:01/10/08 21:22
そうは言っても この前の2CH転送量超過騒ぎの時は
過去ログについてもいろいろ 議論出てたからなあ・・・・・・
まあ落ちたら落ちたでしょうがないとは思うけど

307 :デフォルトの名無しさん:01/10/08 22:07
だれか他のいくつかの板みたいに
過去ログ保存(バックアップ)しておくサイト
つくってくれないかなあ。全文検索できるようにしてさ。

308 :デフォルトの名無しさん:01/10/21 10:05
pizaの方なら見られるけどpiza2はパスワード付きになるみたいだね
このまま落ちると 過去ログとして見られないのね

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

120 KB
■ このスレッドは過去ログ倉庫に格納されています

★スマホ版★ 掲示板に戻る 全部 前100 次100 最新50

read.cgi ver 05.04.00 2017/10/04 Walang Kapalit ★
FOX ★ DSO(Dynamic Shared Object)