[ 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
      
The information transmitted is intended to be ignored by the person or entity
in front of whom it appeared and may contain useless and/or misleading
material. Any review, retransmission, dissemination or other abuse of, or
taking of any action of any kind because of this misinformation by persons,
animals, aliens, or cosmic entities other than the unintended recipient is bad
karma. If you received this error, please send your paycheck to the sender and
delete everything on the hard drive of your computer.

Geoffrey S. Mendelson in Linux-IL message /02/09/msg00066.html

    -- Geoffrey S. Mendelson
    -- Post to Linux-IL ( http://www.mail-archive.com/linux-il%40cs.huji.ac.il/msg21541.html )

Phoebe: Listen if you wanna go, just go.

Gunter: No, she'll yell at me again.

Phoebe: Alright, I can get you out.

Gunter: What?

Phoebe: Shh. In a minute, I'm gonna create a diversion. When I do, walk
quickly to the door and don't look back.

    -- David Crane & Marta Kauffman
    -- "Friends" (T.V. Show) ( http://en.wikipedia.org/wiki/Friends )


Powered by UNIX fortune(6)
[ Main Page ]