BZOJ 2194 快速傅立叶之二

最新推荐文章于 2020-11-24 12:47:06 发布

linjiayang2016

最新推荐文章于 2020-11-24 12:47:06 发布

阅读量192

点赞数

CC 4.0 BY-SA版权

分类专栏： FFT 文章标签： FFT

本文链接：https://blog.youkuaiyun.com/linjiayang2016/article/details/81590205

FFT 专栏收录该内容

7 篇文章

订阅专栏

本文介绍了一种使用快速傅立叶变换（FFT）解决特定数学问题的方法。通过将序列进行翻转和构造，使得原本难以直接应用FFT的问题得以解决。文章详细展示了如何利用FFT计算两个序列的卷积，并给出了具体的实现代码。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

题目大意

已知 $a,b$ 序列，计算 $c_k=\sum{a_i\times b_{i-k}}$ 。

题解

观察题目名称，可以想到FFT……
FFT能解决的是形如下面的式子：
$h_k=\sum f_i\times g_{k-i}$
可以发现， $f$ 数组的下标和 $g$ 数组的下标和不变，且恒为 $k$ 。
小小地转换一下，把 $b$ 数组前后翻转，得：
$c_k=\sum{a_i\times b_{n-1-(i-k)}}=\sum{a_i\times b_{n-1-i+k}}$
这样，和不变，恒为 $n-1+k$ 。然后再构造 $f$ 序列：
$c_k=\sum{a_i\times b_{n-1-i+k}}=f_{n-1+k}$
用FFT计算 $f$ 序列，然后再输出相应的 $c$ 序列就好了。

代码

#include<bits/stdc++.h>
using namespace std;
const int maxn=2148576;
const double pi=acos(-1.0);
struct comp{
    double x,y;
    comp(double xx=0,double yy=0):x(xx),y(yy) {}
    friend comp operator+(const comp &x,const comp &y) {return comp(x.x+y.x,x.y+y.y);}
    friend comp operator-(const comp &x,const comp &y) {return comp(x.x-y.x,x.y-y.y);}
    friend comp operator*(const comp &a,const comp &b) {return comp(a.x*b.x-a.y*b.y,a.x*b.y+b.x*a.y);}
}a[maxn],b[maxn];
int limit=1,r[maxn];
void fft(comp *t,int ty){
    for(int i=0;i<limit;i++)
        if(i<r[i])
            swap(t[i],t[r[i]]);
    for(int mid=1;mid<limit;mid<<=1){
        comp wn(cos(pi/mid),ty*sin(pi/mid));
        for(int j=0,R=(mid<<1);j<limit;j+=R){
            comp w(1,0);
            for(int k=0;k<mid;k++,w=w*wn){
                comp x=t[j+k],y=w*t[j+k+mid];
                t[j+k]=x+y;
                t[j+k+mid]=x-y;
            }
        }
    }
}
int main(void)
{
    int n,l=0;
    scanf("%d",&n);
    while(limit<=(n<<1))
        limit<<=1,++l;
    for(int i=1;i<limit;i++)
        r[i]=((r[i>>1]>>1)|((i&1)<<(l-1)));
    for(int i=0;i<n;++i)
        scanf("%lf%lf",&a[i].x,&b[n-i-1].x);
    fft(a,1);
    fft(b,1);
    for(int i=0;i<limit;i++)
        a[i]=a[i]*b[i];
    fft(a,-1);
    for(int i=n-1;i<(n<<1)-1;i++)
        printf("%d\n",(int)(a[i].x/limit+0.5));
    return 0;
}