icontofig | 发布于 2017-04-01 10:00:46 | 阅读量 299 | FFT
发布于 2017-04-01 10:00:46 | FFT
#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <algorithm>
#include <complex>
using namespace std;
const int maxn = 3e5+5;
const double pi = acos(-1);
int n = 0;
int n1,n2;
int get_num(){
	int num = 0;
	char c;
	bool flag = false;
	while((c = getchar()) == ' ' || c == '\r' || c == '\n');
	if(c == '-')
		flag = true;
	else num = c - '0';
	while(isdigit(c = getchar()))
		num = num * 10 + c - '0';
	return (flag ? -1 : 1) * num;
}
struct FFT{
	complex<double>omega[maxn],iomega[maxn];
	void init(const int n){
		for(int i = 0;i < n;++i){
			omega[i] = complex<double>(cos(2*pi*i/n),sin(2*pi*i/n));
			iomega[i] = conj(omega[i]);
		}
	}
	void trans(complex<double> *a,const int n,const complex<double> *omega){
		int k = 0;
		while((1 << k) < n)k++;
		for(int i = 0;i < n;++i){
			int t = 0;
			for(int j = 0;j < k;j++)if(i&(1 << j))t |= (1 << (k-j-1));
			if(i < t)swap(a[i],a[t]);
		}
		for(int l = 2;l <= n;l *= 2){
			int m = l / 2;
			for(complex<double> *p = a;p != a + n;p += l){
				for(int i = 0;i < m;++i){
					complex<double> x = omega[n/l*i] * p[m+i];
					p[m+i] = p[i] - x;
					p[i] += x;
				}
			}
		}
	}
	void dft(complex<double> *a,const int n){
		trans(a,n,omega);
	}
	void idft(complex<double> *a,const int n){
		trans(a,n,iomega);
		for(int i = 0;i < n;++i)a[i] /= n;
	}
}fft;
int a1[maxn],a2[maxn],ret[maxn];
void mul(const int *a,const int k1,const int *b,const int k2,int *res){
	int n = 1;
	while(n < k1 + k2)n <<= 1;
	static complex<double>c1[maxn],c2[maxn];
	for(int i = 0;i < k1;++i)c1[i].real(a[i]);
	for(int i = 0;i < k2;++i)c2[i].real(b[i]);
	fft.init(n);
	fft.dft(c1,n);fft.dft(c2,n);
	for(int i = 0;i < n;++i)c1[i] *= c2[i];
	fft.idft(c1,n);
	for(int i = 0;i < k1 + k2 - 1;++i)res[i] = static_cast<int>(floor(c1[i].real() + 0.5));
}
int main(){
	n1 = get_num();
	n2 = get_num();
	for(int i = 0;i <= n1;++i)a1[i] = get_num();
	for(int i = 0;i <= n2;++i)a2[i] = get_num();
	mul(a1,n1+1,a2,n2+1,ret);
	for(int i = 0;i <= n1+n2;++i)printf("%d ",ret[i]);
	return 0;
}

链接博客讲解@:Menci‘s Blog FFT


内容更新于: 2017-04-05 14:31:58
链接地址: http://blog.leanote.com/post/icontofig/4dbd29de626b

上一篇: BZOJ 3196 二逼平衡树

下一篇: NOI2005 维修数列

299 人读过
立即登录, 发表评论.
没有帐号? 立即注册
0 条评论
文档导航