[UOJ 34]多项式乘法(FFT)

题目链接

http://uoj.ac/problem/34

思路

裸的FFT模板题。
具体的FFT算法请看我的另一篇闲扯http://blog.csdn.net/qpswwww/article/details/44223331

代码

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <algorithm>
#include <cmath>
#include <complex>

#define MAXN 2100000

using namespace std;

double PI=3.1415926535897384626;

struct Complex //复数
{
    double real,imag; //实部与虚部
    Complex(){}
    Complex(double _real,double _imag):real(_real),imag(_imag){}
}a[MAXN],b[MAXN];

int bitrev[MAXN]; //设w为N次意义下的单位复根，w[i]=w^i

Complex operator+(Complex a,Complex b)
{
    return Complex(a.real+b.real,a.imag+b.imag);
}

Complex operator-(Complex a,Complex b)
{
    return Complex(a.real-b.real,a.imag-b.imag);
}

Complex operator*(Complex a,Complex b)
{
    return Complex(a.real*b.real-a.imag*b.imag,a.real*b.imag+a.imag*b.real);
}

void FFT(Complex a[],int n,int f) //对长度为n的序列(系数或点值序列)a进行快速傅立叶变换，f为1是求值操作(系数->点值)，f为2是插值操作(点值->系数)
{
    for(int i=0;i<n;i++) if(i<bitrev[i]) swap(a[i],a[bitrev[i]]);
    for(int i=1;i<n;i<<=1) //log(n)次迭代模拟递归
    {
        Complex wn=Complex(cos(PI/i),f*sin(PI/i)); //w^i=w*wn^i
        for(int j=0;j<n;j+=(i<<1))
        {
            Complex w=Complex(1,0);
            for(int k=0;k<i;k++,w=w*wn)
            {
                Complex x=a[j+k],y=w*a[j+k+i]; //x,y即为f0(x)和f1(x)
                a[j+k]=x+y,a[j+k+i]=x-y;
            }
        }
    }
    if(f==-1) //如果是插值的话，那么就需要在结果的前面乘上1/n
        for(int i=0;i<n;i++)
            a[i].real/=n;
}

void mul(Complex a[],int n,Complex b[],int m) //长度为n的多项式a乘上长度为m的多项式b
{
    int L=n+m-1; //最终答案的多项式的长度
    int N=1,k=0; //N=2^k，上面的w就是N次单位复根
    for(;N<L;N<<=1,k++);
    for(int i=0;i<N;i++) bitrev[i]=(bitrev[i>>1]>>1)|((i&1)<<(k-1)); //!!!!!!
    FFT(a,N,1);
    FFT(b,N,1);
    for(int i=0;i<N;i++) a[i]=a[i]*b[i]; //!!!!对两个多项式的系数表示法做完FFT后，就得到了它们的点值表示，直接乘上去即可(不是只乘实部！！)
    FFT(a,N,-1); //然后再从点值表示转回到系数表示
}

int main()
{
    int n,m;
    scanf("%d%d",&n,&m); //此时n和m是最高次项的次数
    n++,m++; //此时n和m是多项式a和b的系数的个数
    for(int i=0;i<n;i++)
        scanf("%lf",&a[i].real);
    for(int i=0;i<m;i++)
        scanf("%lf",&b[i].real);
    mul(a,n,b,m);
    for(int i=0;i<n+m-1;i++)
        printf("%d ",(int)(a[i].real+0.5));
    puts("");
    return 0;
}

改成我现在一直用的模板再交了一发

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <algorithm>
#include <cmath>
#include <complex>
#define MAXN 2100000
using namespace std;
double PI=3.1415926535897384626;
struct Complex
{
    double real,imag;
    Complex(){}
    Complex(double _real,double _imag):real(_real),imag(_imag){}
}a[MAXN],b[MAXN];
inline Complex operator+(Complex a,Complex b)
{
    return Complex(a.real+b.real,a.imag+b.imag);
}
inline Complex operator-(Complex a,Complex b)
{
    return Complex(a.real-b.real,a.imag-b.imag);
}
inline Complex operator*(Complex a,Complex b)
{
    return Complex(a.real*b.real-a.imag*b.imag,a.real*b.imag+a.imag*b.real);
}
inline void reverse(Complex a[],int n)
{
    for(int i=1,j=n/2,k;i<n-1;i++)
    {
        if(i<j) swap(a[i],a[j]);
        k=n/2;
        while(j>=k)
        {
            j-=k;
            k>>=1;
        }
        if(j<k) j+=k;
    }
}
inline void FFT(Complex a[],int n,int flag)
{
    reverse(a,n);
    for(int i=1;i<n;i<<=1)
    {
        Complex wn=Complex(cos(PI/i),flag*sin(PI/i));
        for(int j=0;j<n;j+=(i<<1))
        {
            Complex w=Complex(1,0);
            for(int k=0;k<i;k++,w=w*wn)
            {
                Complex x=a[j+k],y=w*a[j+k+i];
                a[j+k]=x+y;
                a[j+k+i]=x-y;
            }
        }
    }
    if(flag==-1)
        for(int i=0;i<n;i++)
            a[i]=Complex(a[i].real/n,a[i].imag);
}
inline void mul(Complex a[],int n,Complex b[],int m)
{
    int L=n+m-1;
    int N=1,k=0;
    for(;N<L;N<<=1,k++);
    FFT(a,N,1);
    FFT(b,N,1);
    for(int i=0;i<N;i++) a[i]=a[i]*b[i];
    FFT(a,N,-1);
}
int main()
{
    int n,m;
    scanf("%d%d",&n,&m);
    n++,m++;
    for(int i=0;i<n;i++)
        scanf("%lf",&a[i].real);
    for(int i=0;i<m;i++)
        scanf("%lf",&b[i].real);
    mul(a,n,b,m);
    for(int i=0;i<n+m-1;i++)
        printf("%d ",(int)(a[i].real+0.5));
    puts("");
    return 0;
}