bzoj2194 快速傅立叶之二

传送门

上面骗人的

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

fft详解(这次是真的)

fft用于快速计算多项式乘法 满足形如c[n] = sigma(a[i] * b[n-i])的式子都可以有效计算

这题是上面形式反过来 

我们可以再返回去

把a,c数组翻转完发现上面的式子就是c[n-k] = sigma(a[n-i] * b[i-k])

所以fft板子粘上就行

注意lim开到n<<1

Code:

#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
#define rep(i,a,n) for(int i = a;i <= n;++i)
#define per(i,n,a) for(int i = n;i >= a;--i)
#define inf 2147483647
#define ms(a,b) memset(a,b,sizeof a)
using namespace std;
typedef double D;
typedef long long ll;
const D PI = acos(-1);
ll read() {
    ll as = 0,fu = 1;
    char c = getchar();
    while(c < '0' || c > '9') {
        if(c == '-') fu = -1;
        c = getchar();
    }
    while(c >= '0' && c <= '9') {
        as = as * 10 + c- '0';
        c = getchar();
    }
    return as * fu;
}
//head
const int N = 1<<21;
struct cp {
    D x,y;
    cp(D X = 0,D Y = 0) {x = X,y = Y;}
    cp operator + (const cp &o) const {return cp(x+o.x,y+o.y);}
    cp operator - (const cp &o) const {return cp(x-o.x,y-o.y);}
    cp operator * (const cp &o) const {return cp(x*o.x-y*o.y,x*o.y+y*o.x);}
} a[N],b[N];

int r[N],lim = 1,base;
void fft(cp *a,int d) {
    rep(i,0,lim-1) if(i < r[i]) swap(a[i],a[r[i]]);
    for(int i = 1;i < lim;i <<= 1) {
        cp T(cos(PI/i),d * sin(PI/i));
        for(int k = 0;k < lim;k += (i<<1)) {
            cp t(1,0);
            rep(l,0,i-1) {
                cp nx = a[k+l],ny = a[k+l+i] * t;
                a[k+l] = nx+ny,a[k+l+i] = nx-ny;
                t = t * T;
            }
        }
    }
}
int n;
int main() {
    n = read();
    while(lim <= n<<1) lim <<= 1,base++;
    rep(i,0,lim-1) r[i] = (r[i>>1]>>1) | ((i&1) << base-1);

    rep(i,0,n-1) a[i] = read(),b[i] = read();
    reverse(a,a+n);

    fft(a,1),fft(b,1);
    rep(i,0,lim-1) a[i] = a[i] * b[i];
    fft(a,-1);
    per(i,n-1,0) printf("%d
", int(a[i].x / lim + 0.5));
    return 0;
}