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

1 ：07/11/24 〜 最終レス ：11/07/10

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

2 ：

てめぇ一人でやれ

3 ：

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

4 ：

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

5 ：

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

6 ：

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

7 ：

やっぱり面倒になった

8 ：

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

9 ：

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

10 ：

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

11 ：

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

12 ：

Ｒ　環に対して
Ｒ［Ｘ］/（Ｘ＾ｎ-1）　→　R[G]　（同型のフーリエ変換）
という写像を作って
多項式環の乗法が、　右側だと簡単で計算出来て、
それを同型で戻すという寸法だ

13 ：

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

14 ：

ググレカス

15 ：

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

16 ：

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

17 ：

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

18 ：

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

19 ：

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

20 ：

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 ：

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

22 ：

lispでも使ってれば？

23 ：

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

24 ：

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

25 ：

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

26 ：

>>9
じゃあ低速でいいからとっとと作れよ、簡単なんだろ。
高速化するのはその後でいいよ。
最低限、自分の発言に責任を持て。

27 ：

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

28 ：

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

29 ：

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

30 ：

離散フーリエ変換の定義
可換環Ｒに対して、次を満たす要素ｅが存在するとする　
・　ｅは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 ：

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

32 ：

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

33 ：

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

34 ：

離散フーリエ変換のまとめ
可換環Ｒに対して、次を満たす要素ｅが存在するとする　
・　ｅは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 ：

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

36 ：

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

37 ：

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

38 ：

素数と位数を調べるコード
#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 ：

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

40 ：

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 ：

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

42 ：

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

43 ：

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

44 ：

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

45 ：

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

46 ：

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

47 ：

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

48 ：

(^ω^i)

49 ：

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

50 ：

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

51 ：

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

52 ：

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

53 ：

Vector stringのようにクラスの動的確保確保ってどうやるのかわかりますか?
たとえば
Lint x;
x=2;
for(int i; i<1000; i++)x=x*x;
x.printf();
としたとき　xのサイズが自動で拡張されるようにしたいのですが

54 ：

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

55 ：

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

56 ：

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

57 ：

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

58 ：

つ「式テンプレート」

59 ：

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

60 ：

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

61 ：

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

62 ：

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

63 ：

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

64 ：

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

65 ：

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 ：

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

67 ：

やっぱりそうだ

68 ：

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

69 ：

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

70 ：

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

71 ：

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

72 ：

(´；ω；`)

73 ：

>>52のような事を言ってるようなレベルの>>1じゃ、
FFTWとかFFTEより速いもんが出来上がる訳もないわな。
ついでにその整数環FFT使った多倍長乗算は普通のFFTの物より遅いでFA出てた筈。
プロセッサの設計・製造からやるんなら話が違うかもしれんがね。

74 ：

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

75 ：

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

76 ：

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

77 ：

引数は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 ：

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

79 ：

これから修正します
あと逆数の計算方法の書いてあるサイトをのせます
1/n を求めるには
http://www.nextftp.com/swlabo/m0_pctech/hp_ultraprecision/up_820.htm
http://yowaken.dip.jp/tdiary/20070428.html

80 ：

>>77
引数1つにしてくれというエラーが出ます
Lint operator +(Lint const& x, Lint const& y) { return Lint(x) += y; }

81 ：

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

82 ：

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

83 ：

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

84 ：

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

85 ：

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

86 ：

インクリメントはこれでは駄目なんですかね? 不正な構造体の演算ってでます
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 ：

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

88 ：

単純な掛け算のコード書いたけど、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 ：

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

90 ：

>>89
インクリメントの部分ですか?　それはコードの中身ではなくて初めの宣言が違っていました　
opPostInc() という名前はD言語用のインクリメントの宣言でした
2^16進数を10進表示にするにはどうしたらいいんだろう? 10で割ってその余りを出力していけばいいけど10で割るのが難しい

91 ：

ダメだこりゃ。

92 ：

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

93 ：

int x=425980;
cout<< (x&65535)<<" "<<(x%65536);
この正解は10のはずなんですけど値が違います　なぜでしょうか

94 ：

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

95 ：

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

96 ：

最初から１０進でどうぞ

97 ：

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

98 ：

>>96 ビットシフトや&が出来なくて駄目です
>>97 詳細を教えて下さい　やっぱり割り算はいると思うのですが
xを表示するなら　x%10を出力して、x=x/10として繰り返すという方法です

99 ：

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

100read 1read

1read 100read

TOP カテ一覧 スレ一覧 2ch元 削除依頼 ▲

Borland Developer Studio 2006 No.13 (266)
Excel VBA 質問スレ Part24 (908)
NullPointerExceptionを「ぬるぽ」と呼ぶスレ6 (330)
Mathematicaプログラミング　質問箱　その1 (332)
C言語を勉強しあおう (252)
[無料でﾗｸﾗｸJava帳票作成] JasperReports使い集合 (315)
--log9.info------------------
スロット、パRのパロディ同人について 3 (229)
Twitter愚痴スレ@同人板 その74 (211)
俺が書き込むと途端に静かなる」に定評のある奴ちょっとこい＠同人 (672)
(゜д゜≡゜д゜)今同人板には自分しかいない20 (681)
友達の新刊の本当の感想 その２ (368)
【アザゼル・サブマス】柴犬【牛歩レポ漫画】 (539)
【海鮮ｲﾝｽﾏｽ】獣コミュニティ118【蟷螂面ｶﾞﾘ作家】 (391)
【郵便事故偽装で中止】key注意喚起スレ【タイバニ空折結婚式】 (255)
◯忍者の玉子同人15井 (1001)
同人絡みで愚痴りたい好かれてるけど嫌いなキャラ・作品・カプ 13 (534)
□■□最終幻想同人　ＤＩＳＣ18　□■□ (835)
【ZXCT】最終幻想7座蔵オンリー公式パクグッズ問題：4【ZINC】 (349)
【戸田】蛹（NAZ/本城塁）に注意15【ｽｶｲｼﾞｬｯｸ】 (912)
アスペルガーだけど同人活動 (268)
おお才辰り　35イニング (976)
田中芳樹総合スレッド5 (468)
--log55.com------------------
★2ch.scは何故失敗したのか
★クロール批判要望スレ
★削ジェンヌに文句ある人集合
★迷惑行為報告担当 - 小さな親切募集中 2
★2ch.scへの要望スレ Part3
★かっこう観測所
★スレ立て人キャップ
★2ch.scニュース系板観測所