離散フーリエ変換を用いた多倍長乗算の話 とか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
<rindolf> sleeper: why are people obsessed with one-liners?
<rindolf> It takes 3 lines - OMG what a disaster!
<Botje> rindolf: newline prices went up again
<rindolf> Botje: I buy my newlines in the black market
<dabreegster> Botje: again? drat.
* Botje reports rindolf to the newline police
<dabreegster> Botje: I know about an... (underground) operation going on
to pirate newlines.
<rindolf> Botje: I bribed a few cops in the newline police, but nice
try.
<dabreegster> Botje: Some crazy guys are trying to free newlines from
patents! They want to rid the market!
* cursor gets called up to serve in the newline jury
<rindolf> I think we need to start a campaign to lift all
restrictions off newlines.
<dabreegster> rindolf: La Resistance lives on!\n
<Botje> I already stockpiled millions of newlines
<dabreegster> Botje: We can have the one-liners destroyed by sundown
<dabreegster> Not destroyed, but... TURNED INTO TWO-LINERS! Mwuhahaha!
-- The Cost of Newlines
-- #perl, Freenode
How much wood would a woodchuck chuck if
a woodchuck would chuck wood?
-- One of Nadav Har'El's Email Signatures.