BZOJ 2194 [快速傅里叶变换 卷积]

 

---

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

 

---

 

卷积 (f x g)(n)=∑{f(i)*g(n-i):i=0...n}

多项式乘法就是一个系数向量的卷积

可以用FFT快速计算卷积

遇到和不是定值的情况可以反转一个向量

如本题反转a向量后

　　c[k]=∑(a[n-i-1]*b[i-k]) k<=i<=n-1

更换求和指标 i=i-k

　　c[k]=∑(a[n-i-k-1]*b[i]) 0<=i<=n-k-1

把-k-1消去，令t=n-k-1

　　c[n-t-1]=∑(a[t-i]*b[i]) 0<=i<=t

这样就是标准的卷积形式啦

[update 2017-03-30]

重做了一下

反转一个向量，变成和为常数的形式

$ c_k = \sum\limits_{i=k}^{n-1} a_i b_{n-1-i+k} = d_{n+k-1} $

这样计算d是没问题的，因为a和b只有$0...n-1$非0，其他都是0

 

#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
using namespace std;
typedef long long ll;
const int N=(1<<18)+5, INF=1e9;
const double PI=acos(-1);
inline int read(){
    char c=getchar();int x=0,f=1;
    while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();}
    while(c>='0'&&c<='9'){x=x*10+c-'0';c=getchar();}
    return x*f;
}

struct meow{
    double x, y;
    meow(double a=0, double b=0):x(a), y(b){}
};
meow operator +(meow a, meow b) {return meow(a.x+b.x, a.y+b.y);}
meow operator -(meow a, meow b) {return meow(a.x-b.x, a.y-b.y);}
meow operator *(meow a, meow b) {return meow(a.x*b.x-a.y*b.y, a.x*b.y+a.y*b.x);}
meow conj(meow a) {return meow(a.x, -a.y);}
typedef meow cd;

struct FFT{
    int n, rev[N];
    void ini(int lim) {
        n=1; int k=0;
        while(n<lim) n<<=1, k++;
        for(int i=0; i<n; i++) {
            int t=0;
            for(int j=0; j<k; j++) if(i&(1<<j)) t |= (1<<(k-1-j));
            rev[i]=t;
        }
    }
    void dft(cd *a, int flag) {
        for(int i=0; i<n; i++) if(i<rev[i]) swap(a[i], a[rev[i]]);
        for(int l=2; l<=n; l<<=1) {
            int m=l>>1; 
            cd wn = meow(cos(2*PI/l), flag*sin(2*PI/l));
            for(cd *p=a; p!=a+n; p+=l) {
                cd w(1, 0);
                for(int k=0; k<m; k++) {
                    cd t = w*p[k+m];
                    p[k+m] = p[k] - t;
                    p[k] = p[k] + t;
                    w=w*wn;
                }
            }
        }
        if(flag==-1) for(int i=0; i<n; i++) a[i].x/=n;
    }
    void mul(cd *a, cd *b, int lim) {
        ini(lim);
        dft(a, 1); dft(b, 1);
        for(int i=0; i<n; i++) a[i]=a[i]*b[i];
        dft(a, -1);
    }
}f;

int n;
cd a[N], b[N];
int main() {
    freopen("in","r",stdin);
    n=read();
    for(int i=0; i<n; i++) a[i].x=read(), b[n-1-i].x=read();
    f.mul(a, b, n+n-1);
    for(int i=n-1; i<2*n-1; i++) printf("%d\n", int(a[i].x+0.5));
}