著作権フリーのC++高速汎用多倍長演算を作るｽﾚ

1 ：デフォルトの名無しさん：2007/11/24(土) 09:53:51

まずは高速フーリエ変換を実装する
ネット上でここが詳しい
http://homepage1.nifty.com/~umetani/ims/Algebra.html

2 ：デフォルトの名無しさん：2007/11/24(土) 09:55:52

てめぇ一人でやれ

3 ：デフォルトの名無しさん：2007/11/24(土) 10:02:03

日本じゃ著作権フリーは無理だから、外人に頼むしかないな。

4 ：1：2007/11/24(土) 10:07:13

ここに書き込んだ人が著作権を主張しないで自由に使って良いっていえばいいだけではないの?

5 ：デフォルトの名無しさん：2007/11/24(土) 10:15:38

作るスレじゃなくて作ってくださいスレの間違いだろ。
単発依頼スレを自分だけ特別に立てて良いと思った理由を述べろ。

6 ：1：2007/11/24(土) 11:03:23

ちゃんと活動するよ　まずは離散フーリエ変換のノートを書き込むからまっててくれ

7 ：1：2007/11/24(土) 11:27:30

やっぱり面倒になった

8 ：デフォルトの名無しさん：2007/11/24(土) 12:02:03

作るやつはスレ立てる前に作り始める
できんやつはスレ立てる．
若干反例があるが，これは普遍の真理．

9 ：1：2007/11/24(土) 12:02:48

なぜかというと理論的な部分は簡単なのだが、高速化がする部分が複雑で既にあるやつ使った方が良いと思うからだ

10 ：デフォルトの名無しさん：2007/11/24(土) 12:08:16

著作権フリーなものを海外から日本へ持ち込んだ場合、
その著作権は誰のものになるんだ？

11 ：デフォルトの名無しさん：2007/11/24(土) 12:08:32

まさか、FFTをやってくれるフリーなライブラリが欲しいだけってことはないよね？
ちなみに、それはもう存在してて、ググれば見つかるよ。
商用利用もOK。GPLじゃない。

12 ：1：2007/11/24(土) 12:15:50

Ｒ　環に対して
Ｒ［Ｘ］/（Ｘ＾ｎ-1）　→　R[G]　（同型のフーリエ変換）

という写像を作って
多項式環の乗法が、　右側だと簡単で計算出来て、
それを同型で戻すという寸法だ

13 ：1：2007/11/24(土) 12:28:08

>>11
では、高速離散フーリエ変換のライブラリでフリーで高速にうごくやつを組み込んで作ろう
どこかにいいやつありますか?

14 ：デフォルトの名無しさん：2007/11/24(土) 12:34:10

ググレカス

15 ：1：2007/11/24(土) 12:49:28

>>14
それより厳密な整数演算に使える高速フーリエ変換ってあるんですか?
誤差が出るやつでは駄目です　整数演算のみでやるやつです

16 ：デフォルトの名無しさん：2007/11/24(土) 16:00:29

多倍長じゃなかったのかよ

17 ：1：2007/11/24(土) 17:43:36

多倍長って2バイトや4バイトや8バイト以上の整数、小数を扱うって事だよ

18 ：デフォルトの名無しさん：2007/11/24(土) 18:30:37

つくってみるわ
とりあえず
基幹部分はC
インターフェイスはPython
でいいんだよな

19 ：デフォルトの名無しさん：2007/11/24(土) 18:34:52

ttp://homepage1.nifty.com/~umetani/ims/Tools.html

20 ：1：2007/11/24(土) 19:29:10

FFT乗算の原理


数値を多項式と見立てて、多項式の乗算、即ち、離散系の畳み込みをフーリエ空間で行い、
それを実空間に戻して、数値の乗算をしたことにするのである。
今、4096個の一次元データに4096個データの一次元フィルタを掛けるとすれば、
実空間では、4096*4096 約1678万回の乗算が必要となるが、
DFTでは、畳み込みだけでは、0.4万回で済む。前後に3回のDFTが必要となるが、
この分は、FFTで行えば8万回程度なので、約8.4万回で済むことになる。約200倍のスピードとなる。
http://www.nextftp.com/swlabo/m0_pctech/hp_ultraprecision/up_815_1.htm

21 ：1：2007/11/24(土) 19:34:58

高速フーリエ変換にかかるコストを除けば、Ｎ*Ｎの乗法がＮに出来る最強アルゴリズムだよな
3/4　Ｎ*Ｎにできるアルゴリズムがあるが断然違うよな

22 ：デフォルトの名無しさん：2007/11/24(土) 19:37:55

lispでも使ってれば？

23 ：デフォルトの名無しさん：2007/11/24(土) 19:48:59

せめてメール蘭にsageって入れようぜ。

24 ：デフォルトの名無しさん：2007/11/24(土) 22:23:19

このスレはFFTを実装するスレだったのか。
多倍長整数の加減乗除を実装するスレなのかと思った。

25 ：デフォルトの名無しさん：2007/11/24(土) 22:32:20

スレタイに汎用って入れているくらいだしね。

26 ：デフォルトの名無しさん：2007/11/24(土) 23:11:52

>>9
じゃあ低速でいいからとっとと作れよ、簡単なんだろ。

高速化するのはその後でいいよ。
最低限、自分の発言に責任を持て。

27 ：デフォルトの名無しさん：2007/11/25(日) 00:32:34

FFTが多倍長の乗算に使われる事を知らない奴ばかりで、
1がDQNあつかいされてるな。
フーリエ変換って、多倍長乗算やCTとか、
周波数とか関係のない、意外なとこで利用されてるから面白いよね。

28 ：デフォルトの名無しさん：2007/11/25(日) 02:31:31

違うわ。αでもβでもいいから物を出さないのがいけない。

29 ：デフォルトの名無しさん：2007/11/25(日) 14:53:53

>>4亀レスだが
2ちゃんねるに書き込んだものってひろゆきが権利持つんじゃないだろうか。

30 ：1：2007/11/25(日) 14:54:26

離散フーリエ変換の定義

可換環Ｒに対して、次を満たす要素ｅが存在するとする　
・　ｅは1の原始ｎ+1乗根　（すなわち、ｅのｎ+1乗は1かつ、ｎ乗以下では1ではない）
・　ｋがｎ以下ならば、　Σ( i, 0, n) E(k) ^ i =0 ただしE(k) ＝ｅ＾ｋ

このとき、Ｇ＝｛1，Ｅ（1），・・・，Ｅ（ｎ）｝は群をなす
ＧからＲへの写像全体をＣ（ｅ）と置くと、次のようにしてR上の有限次元ベクトル空間になる
内積を、(f,g)=Σ( i, 0, n) f( Ｅ（j）) g( Ｅ（-j）)で定義すると、｛λi｝が正規直交基底を与える　（ただしλi(　Ｅ（j）　)＝Ｅ（jj））
さらにＣ（ｅ）は次で与える積でR代数になる　(f*g)( Ｅ（k）) = Σ( i, 0, n) f( Ｅ（k-i）) g( Ｅ（i）)

R[X]の要素fに対して、C(e)の要素Ｆ（ｆ）を次のように定義する　Ｆ（ｆ）( Ｅ（k）) = Σ( i, 0, n) f( Ｅ（i）) Ｅ（ki）
この対応はR代数として準同型写像を与える
これから導かれる　同型　R[X]/(X^(n+1) - 1) ≡　C(e)　を離散フーリエ変換という

31 ：1：2007/11/25(日) 15:00:24

たとえば、Ｒ＝Ｚ/（5）＝｛0，1，2，3，4｝において、2は1の原始4乗根であり、もう一方も条件を満たすので
Ｒ［Ｘ］/（Ｘ^4-1）　からC(2)への離散フーリエ変換が定義できる

32 ：デフォルトの名無しさん：2007/11/25(日) 18:34:53

フーリエ変換を使って高速になるのは桁数がある程度以上大きくなってからでしょう。
だからまずフーリエ変換を使わないプログラムを作るべきでは？

33 ：1：2007/11/25(日) 18:54:30

>>32　まあそういうことですね　どこまでフリーリエ以外が速いか知らないといけないですからね
しかしフーリエで適用できるようにデータ構造は決めないと駄目なんです

34 ：1：2007/11/25(日) 19:06:03

離散フーリエ変換のまとめ

可換環Ｒに対して、次を満たす要素ｅが存在するとする　
・　ｅは1の原始ｎ乗根
・　ｋがｎ未満ならば、　Σ( i, 0, n-1) E(k) ^ i =0 ただしE(k) ＝ｅ＾ｋ

このときR[X]/(x^n - 1)の自己同型写像F

f 　→　Σ( i, 0, n-1)　f( E(i) ) X^i をフーリエ変換という
逆写像は、f 　→　(1/n) Σ( i, 0, n-1)　f( E(n-i) ) X^i で与えられる

このとき、畳み込み積　　(f*g)( Ｅ（k）) = Σ( i, 0, n-1) f( Ｅ（k-i）) g( Ｅ（i）)に対して　F(fg) = F(f) * F(g)が成り立つ
X^s * X^t = X^s (s=t) or 0 (それ以外)を満たす　この性質から畳み込み積の計算は簡単にできる

35 ：1：2007/11/25(日) 19:09:54

34は途中で30で定義していた同型R[X]/(X^n - 1) ≡　C(e)を同一視している部分があります

36 ：1：2007/11/25(日) 19:55:32

まず必要なことは、素数pとZ/(p)を乗法群と見たときの最大位数の元を決定することです
最大位数が大きいほど、計算できる桁数が上がります　しかし、pが大きすぎると計算時間がかかります
どうやって最大位数の大きい元をみつけたらいいんでしょうか?

37 ：1：2007/11/25(日) 20:00:50

p = 2^16 - 15に対して、17はZ/(p)において位数p-1の要素だと判明したのでこれを使うことにします

38 ：1：2007/11/25(日) 20:06:43

素数と位数を調べるコード

#include <iostream>
using namespace std;

main(){
int n,k,N,cnt;
for(N=(1<<16)-1;N>2;N--){
for(n=2;n*n<N;n++){k=N/n; if(k*n==N)break;}
if(n*n<=N)continue;
for(n=2;n<=N/2;n++){
k=n;
for(cnt=1;k!=1;cnt++){k*=n;if(k>N)for(;k>N;k-=N);}
if(cnt+1==N){cout<<"素数 "<<N<<"に対し要素"<<n<<"の位数は"<<cnt<<endl;break;}
}}}

39 ：1：2007/11/25(日) 20:16:07

変更してp=2^16 + 1 と位数pの要素3を使うことにします

40 ：1：2007/11/25(日) 21:55:11

FFT乗法の設計の概要を書いておきます
p=2^16 + 1 とe=3とおく

f(X) = Σa_i ・X^i　(a_i　< 2^16) を多倍長整数とします

これをフーリエ変換すると、 Σf( e^i ) ・X^i　になります

g(X)も多倍長整数とするとフーリエ変換はΣg( e^i ) ・X^i　になります

fとgの畳み込み積は、 Σf( e^i ) g( e^i ) ・X^i　となり

これをフーリエ逆変換すると、Σf( e^(p-1-i) ) g( e^((p-1-i)) ) ・X^i　になります

この係数のうち、2^16であるものは桁上がりを処理して結果が出ます

41 ：1：2007/11/25(日) 21:59:03

簡単に言うと、多倍長整数f, gに対し、
f(1)g(1) 、　f(3)g(3)、　f(3^2)g(3^2)、・・・をmod pの範囲で求めれば良いだけです

42 ：デフォルトの名無しさん：2007/11/25(日) 22:14:09

>>29
もし本当にそうなら47氏逮捕のときにひろゆきも逮捕されてるはず

43 ：1：2007/11/25(日) 22:18:37

まちがえました　フーリエ逆変換は、21846・Σf( 3^(p-1-i) ) g( 3^((p-1-i)) ) ・X^i　　です
Z/(p)での演算から　1/3　＝21846

44 ：デフォルトの名無しさん：2007/11/25(日) 22:21:16

>>42
馬鹿？
金子がWinnyの実装を2ちゃんに書き込んだかよ？

45 ：1：2007/11/26(月) 14:47:41

高速フーリエ変換をメモリも再帰も行列計算も使わない方法を解明した
それは、2進数に対応させるんだ　例えば1024*1024のサイズのフーリエ変換ならば
0〜1023の2進数ごとに一定の手続きをすることで計算が完了する

46 ：1：2007/11/26(月) 14:50:16

すこしはメモリは使うけどね　全ディレクトリの探索みたいな感じで上位ディレトリのパスだけは保存しておくような感じ

47 ：1：2007/11/26(月) 15:00:41

たとえば、
サイズは2Nの、(Xi) = (ω^ij) (xi)というフーリエ変換は、

(X2i) = (ω^2ij) (xi+xN+i)
(X2i+1) = (ω^2ij+j) (xi+xN+iω^N)
に変換できる
これを繰り返したとき、Xの添え字、ωの指数、xの足し算の仕方を
000011101などの2進数から決定できる

48 ：デフォルトの名無しさん：2007/11/26(月) 15:03:18

(^ω^i)

49 ：1：2007/11/26(月) 16:36:20

計算量を考えてみたけど、フーリエ変換にN*Nかかって、これを2回やって畳み込み積がN回だから
単純な総当たり掛け算より多いような気がする　特典は桁上がりの死よりが少ないだけでは?

50 ：1：2007/11/26(月) 16:45:06

通常のフーリエ変換だとN*Nが高速だと、 N* log Nだね
Karatsuba法だとN^1.59だね　やはり高速フーリエ変換が最速か
しかし、より詳細に求めると、 2N* log N　+　N　だから、場合によってはN^1.59のほうが速いかもしれないね

51 ：デフォルトの名無しさん：2007/11/26(月) 17:01:02

良くわからんけど、FFTWじゃ駄目なの？

52 ：1：2007/11/26(月) 17:37:49

>>51
使い方分からないし、整数演算に使えるか分からないので自作してみるよ・・・
あと、自作使用した理由が、多倍長ライブラリのコンパイルが出来なかった為だったのだ
C++言語のみで、他のライブラリに依存しないで単体で動くと導入が簡単なのだ

53 ：1：2007/11/26(月) 17:50:11

Vector stringのようにクラスの動的確保確保ってどうやるのかわかりますか?

たとえば
Lint x;
x=2;
for(int i; i<1000; i++)x=x*x;
x.printf();

としたとき　xのサイズが自動で拡張されるようにしたいのですが

54 ：デフォルトの名無しさん：2007/11/26(月) 17:54:49

std::vector使えば(ry

それで性能上困ることがあれば、
自分でmalloc/free呼んだり、OSのメモリ確保APIを直接呼んだり、
それをメモリプールしてみたり。

55 ：1：2007/11/26(月) 18:03:24

サンクス
vectorは速度的に速いと思うので、内部でvector使うことにします　自分で領域を動的に管理するとバグでやすいですし

56 ：1：2007/11/26(月) 18:25:05

すみません

a = x + y という演算子を定義するとき、x+y演算と、代入演算を分離するとしない場合より代入分
計算コストがかかると思うのですが、tasu(&a, &x, &y)という関数と同じように定義することは出来ないですか?
これなら計算結果を直接格納できるのですが

57 ：デフォルトの名無しさん：2007/11/26(月) 18:41:31

そこらへんは最適化でうまくやってくれるのでは？

58 ：デフォルトの名無しさん：2007/11/26(月) 18:41:44

つ「式テンプレート」

59 ：1：2007/11/26(月) 19:00:39

>>58
サンクス　でも使い方分からなかったよ　
なんか代入にかかる計算程度をケチっても
意味がある計算には変化なさそうなので知ってるやり方でやってみる
わざと、代入が多くなるような計算すれば違いが出るとはおもうけど

60 ：1：2007/11/26(月) 20:21:12

なんか足し算作るだけでも面倒
符号と数値を保存してあるとして、足し算が負になるか正になるかで計算を分類しなくていけない

61 ：1：2007/11/26(月) 20:35:53

計算コストがかかるけど、演算結果の符号は考慮しないで最後に調べれば統一的に処理できるけど・・・

62 ：デフォルトの名無しさん：2007/11/26(月) 20:36:04

このアホはリバーシ１かな

63 ：1：2007/11/26(月) 22:01:58

足し算まで出来た
http://kansai2channeler.hp.infoseek.co.jp/cgi-bin/joyful/img/5338.txt

64 ：デフォルトの名無しさん：2007/11/26(月) 22:10:52

>>63
コメントとインデントつけてくれるとうれしいな

65 ：1：2007/11/26(月) 22:19:29

class Lint
int fugou; //正負
vector< unsigned short int > a;　//一区切りは2バイトのunsigned整数

Lint& operator=(Lint x);　// Lint の代入

Lint& operator=(int x); //整数の代入

int operator>(Lint& x);
int operator<(Lint& x); //絶対値の比較

Lint operator+(Lint& x)　//足し算　
絶対値の大きい方を先頭にして、残りの数値を引いたり足したりする

66 ：1：2007/11/26(月) 22:24:20

除算をどうするかが問題ではある
x / y = x * (1/y) だから1/yを求められればよい
すると、符号の他に小数(の位置)を記録する必要がある
ニュートン法で1/yを計算するのが良いらしいけど
有効桁数までもとめ続けた方が良いのだろうか

67 ：デフォルトの名無しさん：2007/11/26(月) 22:28:46

やっぱりそうだ

68 ：デフォルトの名無しさん：2007/11/26(月) 22:32:57

どうでもいいけどぉ〜符号はせめてぇ〜singじゃないのかぁ〜
あとFFTを理解するの面倒でスレを読んでないんだけど多倍長整数みたいなもんだよね？
Lint& operator=(int x)を見る限り4294967296までしか行けそうにないけど、俺が馬鹿なだけ？

69 ：デフォルトの名無しさん：2007/11/26(月) 22:37:30

＞どうでもいいけどぉ〜符号はせめてぇ〜singじゃないのかぁ〜
歌ってるのか？

70 ：デフォルトの名無しさん：2007/11/26(月) 22:40:55

ああー、signだった・・・
どうでもいいからやんわりと書こうとしただけで洒落ようとかこれっぽっちも思ってなかったのにorz

71 ：デフォルトの名無しさん：2007/11/26(月) 23:30:13

どうでもいいけどぉ〜♪符号はせめてぇ〜♪singじゃないのかぁ〜♪

72 ：デフォルトの名無しさん：2007/11/26(月) 23:33:55

(´；ω；`)

73 ：デフォルトの名無しさん：2007/11/26(月) 23:38:38

>>52のような事を言ってるようなレベルの>>1じゃ、
FFTWとかFFTEより速いもんが出来上がる訳もないわな。

ついでにその整数環FFT使った多倍長乗算は普通のFFTの物より遅いでFA出てた筈。
プロセッサの設計・製造からやるんなら話が違うかもしれんがね。

74 ：デフォルトの名無しさん：2007/11/26(月) 23:48:06

>>73
しーっ、リバーシスレのようになるから触っちゃダメ。

75 ：1：2007/11/26(月) 23:49:35

ふつうのFFTって何? 　Exp(iθ)とか出てくるやつ? 整数演算しようとするのに浮動小数点や複素数演算したら誤差が出る

76 ：デフォルトの名無しさん：2007/11/26(月) 23:52:07

>>74 それもそうだな。要らん事言って済まなかった。

77 ：デフォルトの名無しさん：2007/11/26(月) 23:54:49

引数はLintやLint&よりもconst参照だろ、常考。

ところでoperator +=は？
そっち作った後、それを使って残りを定義するのが普通。
More effective C++ 22項

Lint operator +(Lint const& x, Lint const& y)
{
    return Lint(x) += y;
}

あと、operator<はreturn x > *this;で十分なはず。
さらに!(*this < x)でoperator >=、!(*this > x)でoperator <=も実装できる。
ついでにoperator ==とoperator !=も作ってしまえ。

単項operator +と-も今すぐできるはずだ。

欲を言えばswapが欲しい。

78 ：デフォルトの名無しさん：2007/11/26(月) 23:55:36

>>74
やっぱりリバーシスレと同じ匂い感じた人がいたかｗ

79 ：1：2007/11/27(火) 00:29:23

これから修正します
あと逆数の計算方法の書いてあるサイトをのせます

1/n を求めるには

http://www.nextftp.com/swlabo/m0_pctech/hp_ultraprecision/up_820.htm
http://yowaken.dip.jp/tdiary/20070428.html

80 ：1：2007/11/27(火) 01:13:47

>>77
引数1つにしてくれというエラーが出ます

Lint operator +(Lint const& x, Lint const& y) { return Lint(x) += y; }

81 ：デフォルトの名無しさん：2007/11/27(火) 01:29:47

ここまで自分でなんとかしようと思う1は偉いと思う
応援するよ

82 ：1：2007/11/27(火) 01:29:58

あと、+=を定義して、他のメンバ関数で使おうとすると、未定義というエラーを出します
+は他で使えます　+を定義してから　+=を定義する方向で行こうと思います

83 ：デフォルトの名無しさん：2007/11/27(火) 01:38:29

俺もoperatorのオーバーライドの時の定義の仕方が曖昧なんだけど
うまくまとめてるサイトってないものかね？
http://www.nitoyon.com/vc/tips/multi_precision.htm
こんな感じのサンプル上がってるとこ見て、見よう見まねで実装してるんだけど・・・

84 ：1：2007/11/27(火) 01:52:28

なんかconst付けたり付けなかったりで、未定義になったりします　とりあえず動く方法でやっていきます

85 ：1：2007/11/27(火) 02:19:09

ここまでできた　最新版　掛け算や割り算はまだ
http://kansai2channeler.hp.infoseek.co.jp/cgi-bin/joyful/img/5341.txt

86 ：1：2007/11/27(火) 02:34:46

インクリメントはこれでは駄目なんですかね? 不正な構造体の演算ってでます

opPostInc() {
int i,s=a.size();
for(i=0;i<s;i++){
if(a[i]==65535)a[i]=0;
else {a[i]++;break;}}
if(i==s){a.resize(s+1);a[s]=1;}}

87 ：1：2007/11/27(火) 02:40:19

上はD言語の宣言でした　　普通にvoid operator ++() でした

88 ：1：2007/11/27(火) 03:07:26

単純な掛け算のコード書いたけど、10進表示がまだ出来ていなくてあっているか分からない

Lint operator*(Lint & x){
Lint *p,*q;
if(*this > x){p=&x; q=this;} else {p=this; q=&x;}
int s,t,u;
s=(*p).a.size(); t=(*q).a.size(); u=s*t;
Lint y; y.a.resize(u+1,0);
if(sign + x.sign==0) y.sign=-1;
int i,n; unsigned w,k=0;
for(n=0; n < u; n++){
w=k; k=0;
for(i=0; i<=min(n,s); i++){
w += (*p).a[i] * (*q).a[n-i];
k+=(w>>16); w&=LIMIT;}
y.a[n]=w;}
if(k>0)y.a[u]=k; else y.a.resize(u,0);
return y;}

89 ：デフォルトの名無しさん：2007/11/27(火) 03:34:38

a[i]==65535
ここを
a[i]==999999999
にすれば？

90 ：1：2007/11/27(火) 03:53:38

>>89
インクリメントの部分ですか?　それはコードの中身ではなくて初めの宣言が違っていました　
opPostInc() という名前はD言語用のインクリメントの宣言でした

2^16進数を10進表示にするにはどうしたらいいんだろう? 10で割ってその余りを出力していけばいいけど10で割るのが難しい

91 ：デフォルトの名無しさん：2007/11/27(火) 03:58:05

ダメだこりゃ。

92 ：1：2007/11/27(火) 03:59:44

いい方法が分かりましたよ
小数点を導入して、値を1以下になるようにシフトさせてから10倍していきます
整数になった部分が最上位の10進表示になります　それを0にして繰り返します

93 ：1：2007/11/27(火) 05:44:19

int x=425980;
cout<< (x&65535)<<" "<<(x%65536);

この正解は10のはずなんですけど値が違います　なぜでしょうか

94 ：1：2007/11/27(火) 05:47:35

すみません　間違えました
int xx=4259850;
でした

95 ：1：2007/11/27(火) 06:16:16

>>92は間違えました　1より小さくするときに10^(-n)で割ってシフトしなければいけなかったです　結局割り算がいりますかね?　

96 ：デフォルトの名無しさん：2007/11/27(火) 07:36:29

最初から１０進でどうぞ

97 ：デフォルトの名無しさん：2007/11/27(火) 13:33:20

>>95
10000進の掛け算と足し算でできるのでは？

98 ：1：2007/11/27(火) 16:00:47

>>96 ビットシフトや&が出来なくて駄目です
>>97 詳細を教えて下さい　やっぱり割り算はいると思うのですが

xを表示するなら　x%10を出力して、x=x/10として繰り返すという方法です

99 ：デフォルトの名無しさん：2007/11/27(火) 16:57:55

10進2進変換は2進多倍長の掛け算と足し算でできるのと同じように
2進10進変換が10進多倍長の掛け算と足し算でできるということです。

100 ：デフォルトの名無しさん：2007/11/27(火) 17:00:06

>>80
そのopeartor +は非メンバ関数。

swapはそれだと意味がない。代入より効率よく実装できる可能性があるからみんな定義する。
例えばstd::vectorなら、交換は内部のポインタを取り替えるだけで済む。
（代入ならメモリのコピーをしないといけない）
そしてswapという名前のメンバ関数にするのがお約束。
void swap(Lint& x) {
a.swap(x.a);
std::swap(sign, a.sign); //<algorithm>
}
そしてグローバルにもswapを定義しておく。
void swap(Lint& x, Lint& y) {
x.swap(y);
}
単項operator+と-はこういうメンバ関数でいい。
Lint opreotar +() const {
return *this;
}
Lint opreotar -() const {
Lint x = *this;
x.sign = -1;
return x;
}
っていうかコピーコンストラクタLint(Lint const& x)作れ。対称性からLint(int x)も欲しい。

>>98
こういう整数を十進法として出力するにはそれくらいしかないと思う。

101 ：デフォルトの名無しさん：2007/11/27(火) 17:14:17

コーディングテクニックと、アルゴリズムの話が混在してて、読む気にならん・・・

102 ：デフォルトの名無しさん：2007/11/27(火) 17:27:11

>>86
インクリメントはとりあえず*this + 1でいいよ。

Lint& operator ++() {
return *this += 1;
}

Lint operator ++(int) {
Lint t = *this;
return ++t;
} //あるいはreturn ++Lint(*this);とすれば1行にまとめられる。

とまあ、C++の演算子多重定義は面倒くさいので、
こんな余計なことに頭を使わないためにも、Boost.Operators推奨。
http://www.kmonos.net/alang/boost/classes/operators.html
boost::operators<Lint>とboost::operators<Lint, int>を継承し、
そこのsampleに載っているのと同じ演算子多重定義関数だけを実装すれば、残りは全て揃う。
ただし引数取るものはconst Lint&版とint版の2つを実装すること。

103 ：デフォルトの名無しさん：2007/11/27(火) 17:38:32

boost使うと、ライブラリの価値が50%は下がるな

104 ：デフォルトの名無しさん：2007/11/27(火) 17:53:21

俺はそう思わないけど、そう思う人間が多いなら、後で自分で書き直せばいい。
今はとりあえず動く、使えるようすることに注力してほしい。

105 ：デフォルトの名無しさん：2007/11/27(火) 18:20:02

>>102
<と==を定義して残りの> != <= >=をboostに任せるってやっているけど
これは問題だと思うね。
論理的には正しいが下手をすると処理時間が倍以上かかってしまいそう。

106 ：デフォルトの名無しさん：2007/11/27(火) 20:53:28

そこら辺は大抵最適化されるよバカ

107 ：デフォルトの名無しさん：2007/11/27(火) 20:59:25

>>105
>>77にあるようにオペランドの左右を入れ替えたり、結果を!で論理反転するだけ。
コンパイル時の最適化で消えたり、残ってもせいぜい数命令もあるかどうか。
だから無視できる範囲内。

それよりも>>85の現在のoperator ==（と!=）がよっぽどまずい。
!( (*this < x) || (*this > x) )と2度も比較演算を行っており、
これこそ処理時間が倍かかる。

108 ：デフォルトの名無しさん：2007/11/27(火) 21:17:23

>>107
!( (*this < x ) | ( *this > x ) )これって最適化で、( *this >= x && *this <= x )とかにしてくれるのかな？
まぁ、倍とはいかないけど大きな期待ができるほどじゃないけど・・・

boost前提のライブラリはどうかと思うけど、楽に書き直せる範囲なら
機能実装から先にやっていけばいいと思うな。
比較に処理が2倍かかろうが、どうせ後でプロファイルするんだろうし。

109 ：デフォルトの名無しさん：2007/11/27(火) 21:26:59

>>108
うわ、気付かなかった。
そういえばoperator <って符号を見ていないよな。

110 ：105：2007/11/27(火) 22:07:29

a<=bをa<bとa==bからどうやって構成するのかな？
素直に考えると(a<b)||(a==b)になると思うが
それだと最悪2倍程度時間がかかるだろう

111 ：105：2007/11/27(火) 22:11:00

あっ、そうか!(b<a)か。
勘違いした。

112 ：デフォルトの名無しさん：2007/11/27(火) 22:18:38

>>111
どんまいｗ
operator >の一つだけで一応、全部（>= < <= == !=）実装できるよ。

113 ：105：2007/11/27(火) 23:03:58

<<112
でもそれだと処理時間が...
自分が作るとしたらまずCompare関数を作り
残りはそれを呼び出すようにするね。
それなら処理時間は遅くならないはず。

114 ：デフォルトの名無しさん：2007/11/27(火) 23:36:47

>>113
処理時間は置いといてだから一応ね。
Compare作るくらいなら==作って、そのNOTで!=でいいんじゃないか？

115 ：デフォルトの名無しさん：2007/11/27(火) 23:44:11

>>113
流石に==は実装しないと処理時間がかかる(>>107-108)けど、
比較比較の> >= < <=はどれか1つあれば、残りは速度低下なく実装できる(>>77>>107)。

例えば、a < bが定義されていたとして
>>110で挙げられたa <= bは、!(b < a)と定義すればよい。
(a<b)||(a==b)は悪いけど率直というより愚直。

Boost.Operatorsが比較演算に関してoperator <と==だけを要求しているのは
伊達にやっているからじゃないんだよ。

116 ：105：2007/11/28(水) 00:17:36

a<b <=> Compare(a,b)<0
a==b <=> Compare(a,b)==0
だからCompareを作れば基本的に１つですむということ。
a<bもa==bも実際は内部で同じような処理をするわけだから
同じようなものを2つ作る必要はないでしょう。

117 ：デフォルトの名無しさん：2007/11/28(水) 00:39:11

そっちはわかる。
つい遅くなるっていうほうにだけ反応していた。

118 ：デフォルトの名無しさん：2007/11/28(水) 00:42:22

そろそろoperatorの話は置いといて>>1を待とうぜ。
>>95から出てきてないことだし。

119 ：1：2007/11/28(水) 23:28:41

フーリエ変換は次の値を出来るだけ高速に求めることと同じです　簡易ルーチン作ってみましたがより速くできますか?

N=2^16 、p=N+1、f(X)= a0 +　・・・　+ a(N-1) ・ X^(N-1)　に対して、f(1)、f(3)、・・・、f( 3^(N-1) ) mod p　を求める

#include <iostream>
#include <vector>
#include <time.h>
#define p 65537
unsigned int mod_p(unsigned int& x){
int y=(x&65535)-(x>>16);
if(y<0)y+=p; return y;}

main(){
int N=1<<13, i,j;
cout<<"初期データ生成中・・・\n";
vector< unsigned int > a, x; a.resize(N); x.resize(N);
for(i=0;i<N;i++)a[i]=((rand()&255)<<8)+(rand()&255);

vector< vector< unsigned int > > e;
for(i=0;i<N;i++)e.resize(N); for(i=0;i<N;i++)e[i].resize(N);
j=1; for(i=0;i<N;i++){ e[0][i]=1; e[i][0]=1; e[1][i]=j; j*=3;}
for(i=1;i<N;i++) for(j=1;j<N;j++) e[i][j] = e[1][(i*j)&(N-1)];

cout<<"計算開始・・・ ";
int cl=clock();
for(i=0;i<N;i++) for(j=0;j<N;j++){
x[i] += a[j]*e[i][j]; x[i]=mod_p(x[i]);}
cl=clock()-cl; cout<<cl<<" m秒\n";}

120 ：1：2007/11/28(水) 23:30:26

間違えました　Nは任意で、pは2^16+1でした

121 ：1：2007/11/28(水) 23:48:13

いま思ったんですが、どうやら高速フーリエ変換のアルゴリズムは適用出来ないようです
p=2^16+1に対して、p-2桁の多項式を扱う場合は良いんですが、
実用的なサイズは高々1000桁くらい　少なくとも2^15以上の桁で無いと効果はなさそうです

122 ：1：2007/11/29(木) 00:03:15

pの取り方と計算可能な桁数は対応してますが、例えば1000桁どおしの掛け算ならば
2000に近い素数を見つけて2000進数に変換してから高速フーリエ変換を適用する方法がありますが
前後の基数変換と、フーリエ変換にかかるコストを考えると適用出来る桁は限られそうです
単純に計算した方が速くなるケースが多そうなため

123 ：デフォルトの名無しさん：2007/11/29(木) 00:09:53

アルゴリズムの問題ではないけど、
for(i=0;i<N;i++)e.resize(N);ってこれをループにする必要ない。

124 ：1：2007/11/29(木) 00:39:15

>>123　ｻﾝｸｽ　そうでした
まずは、この計算を高速にするコードを作ることにします　うちのパソコンは鈍いので47秒掛かります
高速フーリエ変換が正しく実装出来れば、最大　(2^15*2^15) 　/　 (2^16 * 16 * 2 + 2^16) =496倍高速になるはずです
0.1秒以内で処理が終わるはずですが・・・
　
#include <iostream>
#include <vector>
#include <time.h>
#define rnd() (((rand()&255)<<8)+(rand()&255));

main(){
int N=1<<15, i,j,n;
vector< unsigned int > a, b, c;
a.resize(N); b.resize(N); c.resize(2*N);
for(i=0;i<N;i++){a[i]=rnd();b[i]=rnd();}

unsigned int k=0, x, cl=clock();
for(n=0; n<2*N; n++){ x=k; k=0;
for(i=0; i<N; i++){x += a[i]*b[n-i]; k+=x>>16; x&=65535;}
c[n]=x;}
cl=clock()-cl; cout<<N<<"桁どおしの掛け算の処理時間 "<<(cl+0.0)/1000<<" 秒\n";}

125 ：デフォルトの名無しさん：2007/11/29(木) 10:37:25

地の文だけじゃなく、コードにもホワイトスペース入れろよ・・・

126 ：デフォルトの名無しさん：2007/11/29(木) 11:22:23

レベル下がりすぎｗ

127 ：デフォルトの名無しさん：2007/11/29(木) 19:09:33

アルゴリズムのことはよくわかんねぇんだけどさ
x += a[i]*b[n-i];
bの大きさはN個で、n = 0 -> 2N. i = 0 -> Nなんだから
n - iのとる範囲は0 -> 2Nで範囲外アクセスじゃないの？

128 ：1：2007/11/29(木) 20:00:43

>>127 修正しました あとminGWとGMPいれて速度比較することにしました　コンパイルに時間が掛かります
#include <iostream>
#include <vector>
#include <time.h>
#define rnd() (((rand()&255)<<8)+(rand()&255));

main(){
int N=1<<15, i, n;
vector< unsigned int > a, b, c;
a.resize(N); b.resize(N); c.resize(2*N);
for(i=0;i<N;i++){a[i]=rnd();b[i]=rnd();}

unsigned int k=0, x, cl=clock();
for(n=0; n<=2*N-2; n++){ x=k; k=0;
for(i=max(0,n-N+1); i<=min(N-1,n); i++)
{x+=a[i]*b[n-i]; k+=x>>16; x&=65535;}
c[n]=x;}
cl=clock()-cl; cout<<N<<"桁どおしの掛け算の処理時間 "<<(cl+0.0)/1000<<" 秒\n";}

129 ：1：2007/11/29(木) 20:27:16

GMP速すぎかも・・・　50000！の計算が2秒しかかからない

#include<iostream>
#include<gmpxx.h>
#include<time.h>
using namespace std;

mpz_class factorial(int n){
mpz_class x= 1;
for(int i=1; i<=n; i++) x*= i;
return x;}

int main(void){
int c=clock();
mpz_class x = factorial(50000);
c=clock()-c;
// cout << x.get_str() << endl;
cout<< (c+0.0)/1000<<" sec \n";
scanf("%d",&c);}

130 ：デフォルトの名無しさん：2007/11/29(木) 20:35:30

>>128
minって関数？もしくはマクロ？
毎回計算させるよりローカルに入れたほうが早くならんかな？
最適化されてもよさそうなところだけど一応。

131 ：1：2007/11/29(木) 20:46:25

>>130
#include<iostream>に入っている関数だと思われ　
なんかGMPは巨大な整数には対応していないのかも・・・　65536進数の65535桁とかは扱えないかも・・　計算途中で落ちる

132 ：デフォルトの名無しさん：2007/11/29(木) 20:49:24

>>131
たしかmin関数はtemplateでしょ。
ローカルに落としたほうが早いよ。（最適化されないのなら）

133 ：1：2007/11/29(木) 20:55:21

65 536 * log(65 536) = 315 652.8　だから、10進で30万桁だからたいしたことはないはずだけど・・
円周率などは100億桁とか計算するはずだから
>>132 上のコードは高速化を目標にしてなくて単純なコードと、フーリエ変換のコードの速度を比較する為にうｐしました

134 ：デフォルトの名無しさん：2007/11/29(木) 21:22:20

そうそう、GMPは数がでかすぎるとやけにすっぱりと諦めるよな。
たぶんメモリを喰いつくすとかの見通しをもとにしてはいるんだろうけど。

135 ：デフォルトの名無しさん：2007/11/30(金) 12:17:31

わかっている人も多いと思うが念のため書いておこう。
VC++でもGMPつかえる。
http://fp.gladman.plus.com/computing/gmp4win.htm
>>133
その程度の桁数では大丈夫なはず
ただ、桁数が32bitで表現できる数を超えるとだめらしい。

136 ：1：2007/12/01(土) 05:09:31

#include<iostream>
#include<bitset>
#include<vector>
#include<time.h>
using namespace std;
#define p 65537
#define ee 3
typedef vector< unsigned int > bsu; bsu X, e; vector < bsu > x;

unsigned int mod(unsigned int z){ int y=(z&65535)-(z>>16); return y<0?p+y:y;}

void Fourier_init(){ int k=1; e.resize(p-1); for(int i=0; i<p-1; i++){ e[i]=mod(k); k*=ee; }
X.resize(1<<16); x.resize(16); for(int i=0; i<16; i++) x[i].resize( (1<<(15-i)) );}

bsu& Fourier(bsu& y){
int i, d=1, b, N=32768; bitset<32> a;
for(i=0; i<N; i++) x[0][i]=mod( y[i] );
#define A a.to_ulong()
for(;;){ if(d<15){
b=16-d-1; N=1<<b;
for(i=0; i<N; i++) x[d][i]=mod( x[d-1][i] + e[(A<<b)&65535] * x[d-1][i+N] );
d++;}

else { X[A]=mod( x[14][0] + e[A] * x[14][1] );
X[A+32768]=mod( x[14][0] + e[A+32768] * x[14][1] );
for(d=14; d>=0; d--)if(a[d])a[d]=0; else{a[d]=1;break;}
if(d<1){if(d==0)d=1; else return X;}}}}

main(){ Fourier_init();
bsu a; a.resize(32768);
for(int i=0;i<32768;i++)a[i]=rand();
int cl=clock(); Fourier(a);Fourier(a); cl=clock()-cl;cout<<(cl+0.0)/1000;}

137 ：1：2007/12/01(土) 05:12:27

フーリエ変換のコード書きましたよ　あっているかは不明ですが
乗法で必要な分のフーリエ変換してみると0.15秒くらいかかりました　
最大桁数32768桁の変換にかかる時間
GMPで3万桁の掛け算実行するとエラーになって分からないんですけど
だれか3万桁が出来るコードをあげてもらえませんか?

138 ：1：2007/12/01(土) 05:22:36

フーリエ変換は
f(X)= a0 +　・・・　+ a(N-1) ・ X^(N-1)　、N=2^16 、p=N+1
に対して、f(1)、f(3)、・・・、f( 3^(N-1) ) mod p　を求める事と同じです
3^nはあらかじめ計算しておけば、単純に計算するとN*N程度の乗法が必要ですが
{ 3^n }の周期性から、f(1)、f(3)、・・・、f( 3^(N-1) )を一斉に求めようとすると
N * log2_N程度の乗法で済みます

139 ：デフォルトの名無しさん：2007/12/01(土) 05:28:37

modとか作らなくても%でよくね？
最近のコンパイラやCPUならそっちのほうが早い気がする。
まぁmodも最適化されて%と同じになるかも知れんが。

GMPとかFFTの話に入れなくてこんなことしか言えない・・・

140 ：デフォルトの名無しさん：2007/12/01(土) 05:31:25

いろいろ駄目すぎて、お話にならない

141 ：1：2007/12/01(土) 05:56:26

>>139
ほんとですね　%でも同じ速度になりました　しかし最適化をうまくやらないコンパイラ用にあった方がいいと思ってます
minGWの最新版で最大の最適化でやりました

142 ：1：2007/12/01(土) 06:00:18

N*Nの計算量のルーチンでは9.8秒かかるところが、フーリエ変換では0.15秒ですから600倍以上の高速化になりました
しかし、答えが違いました　N*Nの方はGMPと答えが一致しています

143 ：1：2007/12/01(土) 06:08:12

重大なバグがある気配　なんかFourierを複数個呼び出しても速度が鈍くなってない・・・ちゃんと計算していないって事ですよね・・・

144 ：デフォルトの名無しさん：2007/12/01(土) 06:14:41

最新のソースはどれ？

145 ：1：2007/12/01(土) 06:36:04

まだ全然出来ていないので最新版とか無いんですよ　一応動くやつが出来たらうｐします　その他上の指摘の変更をしたやつを

146 ：デフォルトの名無しさん：2007/12/01(土) 07:47:01

どう見ても同一犯による犯行です。触るだけ無駄。

147 ：デフォルトの名無しさん：2007/12/01(土) 11:52:18

source forge とか使え。本気で作るんなら、どうせsubversionとか使うようになるだろ。
ここに張られただけのソースを見て、自分で整形までしようとは思わん

148 ：デフォルトの名無しさん：2007/12/01(土) 20:58:52

　　L　I　S　P　で　や　れ

149 ：デフォルトの名無しさん：2007/12/04(火) 06:34:58

起きっぱで眠いテンションに任せて色々書いてみた。

>>33
データ構造は結局どうなってるの？構造決まらないとFFTで積でないよ？
多倍長浮動小数だとしたら指数部を記録しないと。

>>45
そんなことは既に考えられている。というより実質ω^iの全配列を作るようなのは無駄。

>>75
加減乗除のうち加減乗は各要素が整数になることは間違いない。
なので浮動小数を使った一般的なFFTを使っても整数化誤差が例えば0.1未満など保証されていれば問題ない。
除算については除算をどう行うか（筆算のようなやり方か、ニュートン法を使うか）、どの程度の誤差を
認めるか（整数演算でキッチリ末尾までか、浮動小数演算として末尾が１違うくらいは認めるか）次第で変わる。
（>>79を書いてるということは浮動小数型？）

>>90,92-94
自分が引用してきたサイトでも見直して配列と多倍長数とがどう対応してるか再考したほうがいい。
ついでにいうと>>98や>>122の発言を見てるとB進数という言葉を勘違いしてる気がする。

>>120-122
N(要素数)が任意でp固定なんてあり得ない。ωは1の原始p乗根のはずなのでFFTに使う要素数は
p-1で固定しなきゃいけない。ということはNがどんなに小さくてもサイズp-1のFFTと同じ時間がかかる。

>>1
ついでなので。リバーシのときと同じく他人の言うことをちゃんと聴くこと。聞くではなく聴くこと。
「まずは動くモノを作る」というのなら再帰使おうが％使おうが、FFTWやGMPより遅かろうが問題ないんだから。
まずはデバグのしやすさ、見易さ優先。アルゴリズムと違うレベルでの速度を求めるのはそれから。

150 ：1：2007/12/04(火) 07:06:32

pが　2^16+1ならば、最大桁は2^15以下にしないと桁あふれのため計算が間違えます
あとpは固定でも桁数が2^n ごとに小さいサイズの計算で済みますよ　
例えば、桁が8桁しか無ければ8*8行列のフーリエ変換を何度かやれば済みます
pは大きくしたり、小さくしたりすると不便なんです　小さくしたら2^16進数が扱えなくなり大きくすると、32bit内に計算がおさまらなくなります
工夫すれば、別のpでも良いんですが、簡単の為に固定してます

151 ：1：2007/12/04(火) 07:15:37

>>149
効率は良くないですけど、浮動小数点は使わずに格納可能な限り全ての桁を記録しようと思います
たとえば、1/3なら0.3333333・・・なら65536桁まで記録して小数点の位置も記録します
整数と同じ扱いにしようと思います
データ構造は、2^16進数で、0個に一番小さい数を入れます　
もし、非常に大きい数を扱うならば、同じライブラリを2度使って65536^32786進数でやろうとおもいます

152 ：1：2007/12/04(火) 08:57:33

シンプルに書いたけど値が一致しない　どこが駄目なんだろう
#include <iostream>
#include <vector>
using namespace std;
typedef vector< unsigned int > Lint;
#define p 65537
#define N 65536
#define M 1000
unsigned int mod_p(unsigned int z){ int y=(z&65535)-(z>>16); return y<0?p+y:y;}
Lint e; void Fourier_init(){ e.resize(p); int k=1; for(int i=0; i<p; i++){ e[i]=k%p; k*=3;}}

Lint Fourier(Lint f){ Lint F; F.resize(N,0);
for(int n=0; n<N; n++) for(int i=0; i<f.size(); i++)
F[n] = mod_p(F[n] + mod_p( e[(n*i)%N] * f[i] )); return F;}

main(){
Fourier_init(); Lint f, g; f.resize(M,0); g.resize(M,0);
for(int i=0; i<M; i++){ f[i]=rand(); g[i]=rand(); }
// FFT乗法
Lint F=Fourier(f), G=Fourier(g); Lint a; a.resize(N);
for(int i=0; i<N; i++)a[i]=mod_p(21846 * mod_p(F[N-i] * G[N-i]));
for(int i=0; i<N; i++)if(a[i]==N){a[i]=0; a[i+1]++;}
//普通の乗法
Lint b; b.resize(N,0); unsigned int x, k=0;
for(int n=0; n<=2*M-2; n++){ x=k; k=0;
for(int i=max(0,n-M+1); i<=min(M-1,n); i++)
{x+=f[i]*g[n-i]; k+=x>>16; x&=65535;}
b[n]=x;}
for(int i=0;i<100;i++)cout<<a[i]<<" "<<b[i]<<"\n";}

153 ：1：2007/12/04(火) 09:16:45

よりシンプルにしてみたけど駄目だ
#include <iostream>
#include <vector>
using namespace std;
typedef vector< unsigned int > Lint;
#define p 65537
#define N 65536
#define M 2
Lint e; void Fourier_init(){ e.resize(p); int k=1; for(int i=0; i<p; i++){ e[i]=k%p; k*=3;}}

Lint Fourier(Lint f){ Lint F; F.resize(N,0);
for(int n=0; n<N; n++) for(int i=0; i<f.size(); i++)
F[n]=(F[n] + ( e[(n*i)%N] * f[i] )%p)%p; return F;}

main(){
Fourier_init(); Lint f, g; f.resize(M,0); g.resize(M,0);f[0]=f[1]=g[0]=g[1]=2;

// FFT乗法
Lint F=Fourier(f), G=Fourier(g); Lint a; a.resize(N);
for(int i=0; i<N; i++)a[i]=(21846 * (F[N-i] * G[N-i])%p)%p;
for(int i=0; i<N; i++)if(a[i]==N){a[i]=0; a[i+1]++;}
//普通の乗法
Lint b; b.resize(N,0); unsigned int x, k=0;
for(int n=0; n<=2*M-2; n++){ x=k; k=0;
for(int i=max(0,n-M+1); i<=min(M-1,n); i++)
{x+=f[i]*g[n-i]; k+=x>>16; x&=65535;}
b[n]=x;}
for(int i=0;i<10;i++)cout<<a[i]<<" "<<b[i]<<"\n";}

154 ：デフォルトの名無しさん：2007/12/16(日) 00:17:52

がんばれー

155 ：デフォルトの名無しさん：2007/12/16(日) 18:09:35

とりあえず
（１）F=Fourier(f),
（２）G=Frourier(g)
としたとき
（３）A=FxG (要素毎の積)
を求めた後に
（４）a=(1/N)Fourier^(-1)(A)
とやる必要があるな。>>153では（３）と（４）が抜けてる。

156 ：デフォルトの名無しさん：2007/12/28(金) 15:00:05

> pは大きくしたり、小さくしたりすると不便なんです　小さくしたら2^16進数が扱えなくなり大きくすると、32bit内に計算がおさまらなくなります
> 工夫すれば、別のpでも良いんですが、簡単の為に固定してます

なら精度考慮してdoubleとかを使うFFTのほうが良くない？
あと、2^16進数使うなら大体(2^16)^2*要素数を表現できるだけの精度が必要になるはずなんだけど。

なんかあちこち変じゃない？

157 ：デフォルトの名無しさん：2008/01/12(土) 21:35:14

1から5000までの数字を全部掛けたときに
末尾に並ぶ0の個数を求めるのをおながいします

158 ：デフォルトの名無しさん：2008/01/12(土) 22:04:34

1000の倍数が5個、100の倍数が50個、10の倍数が500個、555個
5000の倍数が1個、500の倍数が10個、50の倍数が100個、5の倍数が1000個
全部足して1666個と適当に予想して見る。

ってか>>1なのかな？

159 ：デフォルトの名無しさん：2008/01/12(土) 22:53:50

>>157
全く同じ趣旨の問題がVBAスレかどっかでみた。
それは兎も角、鼬害。

>>158
末尾に0が並ぶということは、並んだ数をnとすると10のn乗の倍数ということ。
つまり、2のn乗*5のn乗の倍数になる。
このnを求めればいいわけだが、2の倍数は過剰に存在しているから
5の倍数がどれだけあって、それぞれ5の何乗の倍数か判ればいいことになる。
即ち、5の倍数の個数+25の倍数の個数+125の倍数の個数+……
だと思ったんだけどね。

160 ：デフォルトの名無しさん：2008/01/13(日) 00:17:05

#include<stdio.h>
int main(void){
int i,n=1,o=0;
for(i=2;i<5001;i++){
n*=i;
while(n%10==0)
n/=10,o++;
n%=10;
}
printf("%d\n",o);
return 0;
}

1240になったけど、これであってる？

161 ：デフォルトの名無しさん：2008/01/13(日) 00:44:23

>>160
それ途中でオーバーフローしないか？

162 ：デフォルトの名無しさん：2008/01/13(日) 03:37:54

下一桁だけ計算してるから
大丈夫なんじゃない？

163 ：デフォルトの名無しさん：2008/01/13(日) 03:41:54

>>159 の言う通りに計算してみると
5000 / 5 = 1000
5000 / 25 = 200
5000 / 125 = 40
5000 / 625 = 8
5000 / 3125 = 1

1249個になるね

164 ：デフォルトの名無しさん：2008/01/13(日) 09:28:47

計算してみた。
--
4228577926605543522201064200233584405390786674626646748849782402181358052708108200690899047871706387
5370847466573006854458784860666838127363372108937727876312793903630584621606439044789869822398719297
.
.全部で16300桁以上
.
3713673959956789408847059768595145050172415177460173514309909726155093783347200000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
00000000000000000000000000
--
ということで、1249だった。

165 ：デフォルトの名無しさん：2008/01/13(日) 15:32:26

>>160のコードをn%=100にしたら1249になった。

166 ：デフォルトの名無しさん：2008/01/19(土) 03:20:07

1じゃないけと。
剰余類によるFFTとdoubleによるFFTで乗算作ってみた。

C++で作ったんだけど、どうがんばっても、doubleのほうが早い。

とにかく法mがネック。
任意のmを用いることができないので設計が大変。

2^k+1の形にしておかないと、フーリエできないし、
素数じゃないといけないし、剰余が早くないといけない。
俺は
2^32(2^32-1)+1 2^32*(2^32-4)+1 2^16(2^16-16)+1 2^16(2^16-33)+1
などを使った。

で、剰余は結構早くできるようになったんだが、それでも遅い。
剰余が入るとdoubleのFFTに比べて大体2倍になる。

つぎに最大値の問題。
乗算後の各項の係数の最大値は「項数*(基数-1)^2」になるので、
法mはそれ以上であるか、もしくは法を二つ用意して2回FFTして、
中国剰余定理から求めないといけない。これでさらに2倍かかる。

結局、doubleで精度を考慮して作ったのと剰余類でがんばって作ったのだと、
4倍くらい速度差がでる。

>>73の言うとおりだった。
俺の実装が悪いだけかも知れんので、できればFAのソースが欲しい。