[ Main Page ]

FFTW3 多倍長乗算

離散フーリエ変換を用いた多倍長乗算の話 とか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 )


Powered by UNIX fortune(6)
[ Main Page ]