【BZOJ】2194: 快速傅立叶之二

最新推荐文章于 2024-10-29 21:24:58 发布

weixin_33725722

最新推荐文章于 2024-10-29 21:24:58 发布

阅读量105

点赞数

文章标签：数据结构与算法

原文链接：http://www.cnblogs.com/onioncyc/p/8418113.html

版权

本文介绍了一种使用快速傅里叶变换（FFT）处理卷积的方法，以解决特定的数学问题。具体而言，该方法应用于计算C[k]=sigma(a[i]*b[i-k])的场景，其中a和b为两个数组，通过将B数组反转并扩展上下界，最终转化为适用于FFT的标准形式。

摘要由CSDN通过智能技术生成

【题意】请计算C[k]=sigma(a[i]*b[i-k]) 其中 k < = i < n ，并且有 n < = 10 ^ 5。 a,b中的元素均为小于等于100的非负整数。

【算法】快速傅里叶变换（FFT）处理卷积

【题解】题目要求：

$$C_k=\sum_{i=k}^{n-1}A_i*B_{i-k}$$

令B'为B的数组反转，转化成和为定值的形式，即：

$$C_k=\sum_{i=k}^{n-1}A_i*B'_{n+k-1-i}$$

因为：

$$A_i=0,i\in[n,n+k-1]$$

$$B_{n+k-1-i}=0,i\in[0,k-1]$$

所以可以扩展式子的上下界，即：

$$C_k=D_{n+k-1}=\sum_{i=0}^{n+k-1}A_i*B'_{n+k-1-i}$$

这就可以用FFT处理卷积了。

复杂度O(n log n)。

#include<cstdio>
#include<cmath>
#include<cstring>
#include<algorithm>
using namespace std;
const int maxn=300010;
const double PI=acos(-1);
int n;
struct cp{
    double x,y;
    cp(double a,double b){x=a;y=b;}
    cp(){x=y=0;}
    cp operator + (cp a){return (cp){x+a.x,y+a.y};}
    cp operator - (cp a){return (cp){x-a.x,y-a.y};}
    cp operator * (cp a){return (cp){x*a.x-y*a.y,x*a.y+y*a.x};}
}a[maxn],b[maxn];
void fft(cp *a,int n,int f){
    int k=0;
    for(int i=0;i<n;i++){
        if(i>k)swap(a[i],a[k]);
        for(int j=n>>1;(k^=j)<j;j>>=1);
    }
    for(int l=2;l<=n;l<<=1){
        int m=l>>1;
        cp wn=(cp){cos(2*PI*f/l),sin(2*PI*f/l)};
        for(cp *p=a;p!=a+n;p+=l){
            cp w=(cp){1,0};
            for(int i=0;i<m;i++){
                cp t=p[i+m]*w;
                p[i+m]=p[i]-t;
                p[i]=p[i]+t;
                w=w*wn;
            }
        }
    }
    if(f==-1){for(int i=0;i<n;i++)a[i].x/=n;}
}
int main(){
    scanf("%d",&n);
    for(int i=0;i<n;i++)scanf("%lf %lf",&a[i].x,&b[n-i-1].x);
    int N=1;while(N<n*2)N*=2;
    fft(a,N,1);fft(b,N,1);
    for(int i=0;i<N;i++)a[i]=a[i]*b[i];//2n
    fft(a,N,-1);
    for(int i=n-1;i<2*n-1;i++)printf("%d\n",(int)(a[i].x+0.1));
    return 0;
}