洛谷.3803.[模板]多项式乘法(FFT)

最新推荐文章于 2022-12-22 15:00:56 发布

weixin_30593261

最新推荐文章于 2022-12-22 15:00:56 发布

阅读量92

点赞数

原文链接：http://www.cnblogs.com/SovietPower/p/8444779.html

版权

题目链接：洛谷、LOJ.
FFT相关：快速傅里叶变换(FFT)详解、FFT总结、从多项式乘法到快速傅里叶变换.

5.4 又看了一遍，这个也不错。
2019.3.7 叕看了一遍，推荐这个。

#include <cmath>
#include <cctype>
#include <cstdio>
#include <algorithm>
#define gc() getchar()
const int N=1e6+5;
const double PI=acos(-1);

int n,m;
struct Complex
{
    double x,y;
    Complex(double xx=0,double yy=0) {x=xx, y=yy;}
    Complex operator + (const Complex &a) {return Complex(x+a.x, y+a.y);}
    Complex operator - (const Complex &a) {return Complex(x-a.x, y-a.y);}
    Complex operator * (const Complex &a) {return Complex(x*a.x-y*a.y, x*a.y+y*a.x);}
}A[N*3],B[N*3];//size!

void Fast_Fourier_Transform(Complex *a,int lim,int opt)
{
    for(int j=0,i=0; i<lim; ++i)
    {
        if(i>j) std::swap(a[i],a[j]);
        for(int l=lim>>1; (j^=l)<l; l>>=1);
    }
    for(int i=2; i<=lim; i<<=1)//最后等于lim即整个序列的合并 
    {
        int mid=i>>1;
        Complex Wn(cos(2.0*PI/i),opt*sin(2.0*PI/i)),t;
        for(int j=0; j<lim; j+=i)
        {
            Complex w(1,0);
            for(int k=0; k<mid; ++k,w=w*Wn)
                a[j+mid+k]=a[j+k]-(t=w*a[j+mid+k]),
                a[j+k]=a[j+k]+t;
        }
    }
}
inline int read()
{
    int now=0,f=1;register char c=gc();
    for(;!isdigit(c);c=gc()) if(c=='-') f=-1;
    for(;isdigit(c);now=now*10+c-'0',c=gc());
    return now*f;
}

int main()
{
    n=read(),m=read();
    for(int i=0; i<=n; ++i) A[i].x=(double)read();//scanf("%lf",&A[i].x);
    for(int i=0; i<=m; ++i) B[i].x=(double)read();//scanf("%lf",&B[i].x);
    int lim=1;
    while(lim<=n+m) lim<<=1;
    Fast_Fourier_Transform(A,lim,1);
    Fast_Fourier_Transform(B,lim,1);
    for(int i=0; i<=lim; ++i) A[i]=A[i]*B[i];//size!
    Fast_Fourier_Transform(A,lim,-1);
    for(int i=0; i<=n+m; ++i) printf("%d ",(int)(A[i].x/lim+0.5));
    
    return 0;
}

递归实现：

#include <cmath>
#include <cstdio>
#include <cctype>
#include <cstdlib>
#define gc() getchar()
const int N=2e6+5;
const double PI=acos(-1.0);

int n,m;
struct Complex
{
    double x,y;
    Complex(double xx=0,double yy=0) {x=xx, y=yy;}
    Complex operator + (const Complex &a) {return Complex(x+a.x, y+a.y);}
    Complex operator - (const Complex &a) {return Complex(x-a.x, y-a.y);}
    Complex operator * (const Complex &a) {return Complex(x*a.x-y*a.y, x*a.y+y*a.x);}
}A[N],B[N];

inline int read()
{
    int now=0,f=1;register char c=gc();
    for(;!isdigit(c);c=gc()) if(c=='-') f=-1;
    for(;isdigit(c);now=now*10+c-'0',c=gc());
    return now*f;
}
void Fast_Fourier_Transform(Complex *a,int lim,int type)
{
    if(lim==1) return;
    Complex a1[lim>>1],a2[lim>>1];//爆栈 
    for(int i=0; i<lim; i+=2)
        a1[i>>1]=a[i], a2[i>>1]=a[i+1];
    Fast_Fourier_Transform(a1,lim>>1,type),
    Fast_Fourier_Transform(a2,lim>>1,type);
    Complex Wn(cos(2.0*PI/lim),type*sin(2.0*PI/lim)),w(1,0),t;//Wn:单位根 w:幂
    for(int i=0; i<(lim>>1); ++i,w=w*Wn)
        a[i]=a1[i]+(t=w*a2[i]),
        a[i+(lim>>1)]=a1[i]-t;
}

int main()
{
    n=read(),m=read();
    for(int i=0; i<=n; ++i) A[i].x=read();
    for(int i=0; i<=m; ++i) B[i].x=read();
    int lim=1;
    while(lim<=n+m) lim<<=1;
    Fast_Fourier_Transform(A,lim,1);
    Fast_Fourier_Transform(B,lim,1);
    for(int i=0; i<=lim; ++i) A[i]=A[i]*B[i];
    Fast_Fourier_Transform(A,lim,-1);
    for(int i=0; i<=n+m; ++i) printf("%d ",(int)(A[i].x/lim+0.5));

    return 0;
}

转载于:https://www.cnblogs.com/SovietPower/p/8444779.html