離散フーリエ変換を用いた多倍長乗算の話 とか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
Q: What is the difference between a duck?
A: One leg is both the same.
For example, I was taught in college that one ought to figure out a program
completely on paper before even going near a computer. I found that I did not
program this way. I found that I liked to program sitting in front of a
computer, not a piece of paper. Worse still, instead of patiently writing out
a complete program and assuring myself it was correct, I tended to just spew
out code that was hopelessly broken, and gradually beat it into shape.
Debugging, I was taught, was a kind of final pass where you caught typos and
oversights. The way I worked, it seemed like programming consisted of
debugging.
For a long time I felt bad about this, just as I once felt bad that I didn't
hold my pencil the way they taught me to in elementary school. If I had only
looked over at the other makers, the painters or the architects, I would have
realized that there was a name for what I was doing: sketching. As far as I
can tell, the way they taught me to program in college was all wrong. You
should figure out programs as you're writing them, just as writers and
painters and architects do.
-- Paul Graham
-- "Hackers and Painters" ( http://www.paulgraham.com/hp.html )