FFT 入门

推荐博客 : https://oi.men.ci/fft-notes/

卷积的理解 : https://www.zhihu.com/question/22298352?rf=21686447

 

题目链接 :http://uoj.ac/problem/34

 

这是一道模板题。

给你两个多项式,请输出乘起来后的多项式。
输入格式

第一行两个整数 nn 和 mm,分别表示两个多项式的次数。

第二行 n+1n+1 个整数,表示第一个多项式的 00 到 nn 次项系数。

第三行 m+1m+1 个整数,表示第二个多项式的 00 到 mm 次项系数。
输出格式

一行 n+m+1n+m+1 个整数,表示乘起来后的多项式的 00 到 n+mn+m 次项系数。
样例一
input

1 2
1 2
1 2 1

output

1 4 5 2

explanation

(1+2x)?(1+2x+x2)=1+4x+5x2+2x3(1+2x)?(1+2x+x2)=1+4x+5x2+2x3。
限制与约定

0≤n,m≤105,保证输入中的系数大于等于 0 且小于等于 9。

时间限制:1s1s

空间限制:256MB

 

题意  :  给你两个多项式的系数,从 0 到 n 给出,求这两个多项式相乘后的系数,从小到大输出

思路分析 : 裸的 FFT ,参考kuangbin 的板子

    就是要注意以下数组的大小,main中的 len 是 2^k , 因此当m+n = 2e5 左右时,此时 2^k = 260000+ , 因此要注意数组的大小

代码示例:

#include<bits/stdc++.h>
using namespace std;
#define ll long long
const int maxn = 2e5+63000;
const double pi = acos(-1.0);
int n, m;
struct Complex{
	double x, y;
	Complex (double _x=0, double _y=0):x(_x), y(_y){}
	Complex operator -(const Complex &b)const{
		return Complex(x-b.x, y-b.y);
	}
	Complex operator +(const Complex &b)const{
		return Complex(x+b.x, y+b.y);
	}
	Complex operator *(const Complex &b)const{
		return Complex(x*b.x-y*b.y, x*b.y+y*b.x);
	}
};

Complex x1[maxn], x2[maxn];
int sum[maxn];
void change(Complex y[], int len){
	for(int i = 1, j = len/2; i < len-1; i++){
		if (i < j) swap(y[i], y[j]);
		int k = len/2;
		while(j >= k){
			j -= k;
			k /= 2;
		}
		if (j < k) j += k;
	}
}

void fft(Complex y[], int len, int 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].x /= len;
	}
}

int main () {
	cin >> n >> m;
	int len = 1;
	while(len <= (n+m)) len <<= 1;
	for(int i = 0; i <= n; i++) scanf("%lf", &x1[i].x);
	fft(x1, len, 1);
	for(int i = 0; i <= m; i++) scanf("%lf", &x2[i].x);
	fft(x2, len, 1);
	for(int i = 0; i < len; i++)
		x1[i] = x1[i]*x2[i];
	fft(x1, len, -1);
	for(int i = 0; i <= n+m; i++){
		sum[i] = (int)(x1[i].x+0.5); // sum[] 是最后的答案 
		printf("%d%c", sum[i], i ==n+m?'\n':' ');
	}
	return 0;
}

 

____________________________________________________________________________

int rev[maxl];
void get_rev(int bit)//bit表示二进制位数,计算一个数在二进制翻转之后形成的新数
{
    for(int i=0;i<(1<<bit);i++)
        rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
}
void fft(cd *a,int n,int dft)//n表示我的多项式位数
{
    for(int i=0;i<n;i++) if(i<rev[i]) swap(a[i],a[rev[i]]);
    //中间的那个if保证了每个数做多只被交换了1次
    //如果不写那么会有一些数被交换两次,导致最终的位置没有变
    for(int step=1;step<n;step<<=1)//模拟一个合并的过程
    {
        cd wn=exp(cd(0,dft*PI/step));//计算当前单位复根
        for(int j=0;j<n;j+=step<<1)
        {
            cd wnk(1,0);//计算当前单位复根
            for(int k=j;k<j+step;k++)
            {//蝴蝶操作
                cd x=a[k];
                cd y=wnk*a[k+step];
                a[k]=x+y;//这就是上文中F(x)=G(x)+ωH(x)的体现
                a[k+step]=x-y;
                    //后半个“step”中的ω一定和“前半个”中的成相反数
                    //“红圈”上的点转一整圈“转回来”,转半圈正好转成相反数
                    //一个数相反数的平方与这个数自身的平方相等..
                wnk*=wn;
            }
        }
    }
    if(dft==-1) for(int i=0;i<n;i++) a[i]/=n;
    //考虑到如果是IDFT操作,整个矩阵中的内容还要乘上1/n  
}

 

posted @ 2018-08-18 18:53  楼主好菜啊  阅读(228)  评论(0编辑  收藏  举报