FFT快速傅里叶模板

标签：算法   数学   

FFT快速傅里叶模板……

/*	use way：
		assign : h(x) = f(x) * g(x)  f(x):len1  g(x):len2
		1.  len = 1;  while（len < 2 * len1 || len < 2 * len2) len <<= 1;  
		2.  for i=0 to len1-1 : x1[i](f(i),0)  for i=len1 to len-1 : x1[i](0.0)   g(x) is same.....
		3.  fft(x1,len,1) fft(x2,len,1)
		4.  for i=0 to len-1 : x1[i] = x1[i] * x2[i]
		5.  fft(x1,len,-1)
		6.  ans[i] = (long long)(x[i].a + 0.5)
		ps : goback should from len1 + len2 - 1 ,but not len !  I don't know why.....
*/


#include <vector>
#include <list>
#include <map>
#include <set>
#include <stack>
#include <algorithm>
#include <iostream>
#include <cstdio>
#include <cmath>
#include <cstring>
#include <cstdlib>

using namespace std;

const double pi = atan(1.0) * 4;

struct complex {
	double a, b;
	complex(double aa = 0.0, double bb = 0.0) { a = aa; b = bb; }
	complex operator +(const complex &e) { return complex(a + e.a, b + e.b); }
	complex operator -(const complex &e) { return complex(a - e.a, b - e.b); }
	complex operator *(const complex &e) { return complex(a * e.a - b * e.b, a * e.b + b * e.a); }
};

void change(complex * y, long long len) {
	long long i, j, k;
	for (i = 1, j = len / 2; i < len - 1; i++) {
		if (i < j) swap(y[i], y[j]);
		k = len / 2;
		while (j >= k) {
			j -= k;
			k /= 2;
		}
		if (j < k) j += k;
	}
}

void fft(complex *y, long long len, long long on) {
	change(y, len);
	for (int h = 2; h <= len; h <<= 1) {
		complex wn(cos(-on * 2 * pi / h), sin(-on * 2 * pi / h));
		for (int j = 0; j < len; j += h) {
			complex w(1, 0);
			for (int k = j; k < j + h / 2; k++) {
				complex u = y[k];
				complex t = w * y[k + h / 2];
				y[k] = u + t;
				y[k + h / 2] = u - t;
				w = w * wn;
			}
		}
	}
	if (on == -1)
		for (int i = 0; i < len; i++)
			y[i].a /= len;
}

int main(){
	return 0;
}

FFT快速傅里叶模板

标签：算法   数学   

原文地址：http://blog.csdn.net/frosero/article/details/47430763