離散フーリエ変換を用いた多倍長乗算の話 とかFFTを用いて多倍長乗算をしてみたとか FFTによる多倍精度の乗算とか 色々あるけれど、FFTWを使って実装したものが見当たらないので。 x[]とy[]の長さは2の累乗だと計算が早いらしい。また、x[]とy[]は後ろ半分以上を0で埋めます。
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <fftw3.h>
int main(int argc, char* argv[])
{
// 231 x 11
float x[] = {2.0f, 3.0f, 1.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f,};
float y[] = {1.0f, 1.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f,};
int n = sizeof(x)/sizeof(float);
int m = sizeof(y)/sizeof(float);
{
fftwf_plan p, q;
p = fftwf_plan_r2r_1d(n, x, x, FFTW_R2HC, FFTW_ESTIMATE);
q = fftwf_plan_r2r_1d(n, y, y, FFTW_R2HC, FFTW_ESTIMATE);
fftwf_execute(p);
fftwf_execute(q);
fftwf_destroy_plan(p);
fftwf_destroy_plan(q);
}
x[0] *= y[0];
x[n/2] *= y[n/2];
for (int j = 1; j < n/2; j++)
{
float e = x[j];
float d = x[n - j];
float f = y[j];
float g = y[n - j];
x[n - j] = e*g+f*d;
x[j] = e*f-d*g;
}
{
fftwf_plan p;
p = fftwf_plan_r2r_1d(n, x, x, FFTW_HC2R, FFTW_ESTIMATE);
fftwf_execute(p);
fftwf_destroy_plan(p);
}
for (int j = 0; j < n; j++)
{
x[j] = x[j]/n;
}
float xx = x[0];
for(int j = 0;j < n;j ++)
{
// 桁の処理は省略
fprintf(stderr, "[%4d]%f\n", j, x[j]);
}
return 0;
}
cc -lfftw3f -lm
[ 0]2.000000
[ 1]5.000000
[ 2]4.000000
[ 3]1.000000
[ 4]-0.000000
[ 5]0.000000
[ 6]0.000000
[ 7]0.000000
bzr is slower than Subversion in combination with Sourceforge.
[Dazjorz](http://dazjorz.com/) (17-September-2009)
-- Dazjorz
-- Chat with Shlomi Fish
There is no IGLU Cabal. Some Politically-Correct ISPs block E-mail which comes
from Hell. Unfortunately, some of the brightest minds needed for the IGLU
Cabal languish in Hell.
Omer Zak in Hackers-IL message No. 2203
("Do you want to send E-mail from hell?")
-- Omer Zak
-- Hackers-IL Message No. 2203 ( http://tech.groups.yahoo.com/group/hackers-il/message/2203 )