bzoj2194 快速傅立叶之二

Posted 2021-03-14 wxyww

题目链接

problem

给出两个长度为n的数列a，b。求一个数列c满足：[c[k] = sumlimits_{i = k} ^ na[i]b[i - k]]

(nle 10^5)

solution

长得和卷积很像，观察一下卷积的形式：(c[k]=sumlimits_{i=0}^ia[i]b[k-i])

所以先把b数组翻转过来。

然后所求的式子就变成了(c[k]=sumlimits_{i=k}^na[i]b[n+k-i])

将(n+k)看做整体，那么形式就和卷积比较像了,令(d[n+k]=c[k]=sumlimits_{i=k}^na[i]b[n + k - i])。因为a的(n+1)到(n+k)位置均为0，b的(0)到(k-1)位置均为0.

所以上式也可以写为(d[n+k]=sumlimits_{i=0}^{n+k}a[i]b[n+k-i])。

这样就可以借助FFT计算了。

注意最后的答案是问c[k]，也就是d[k+n]

code

/*
* @Author: wxyww
* @Date:   2020-01-22 08:04:15
* @Last Modified time: 2020-01-22 08:24:04
*/
#include<cstdio>
#include<iostream>
#include<cstdlib>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<queue>
#include<vector>
#include<ctime>
using namespace std;
typedef long long ll;
const int N = 400010;
const double pi = 3.1415926535897931;
ll read() {
    ll x=0,f=1;char c=getchar();
    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 complex {
    double x,y;
}A[N],B[N];

complex operator * (const complex &A,const complex &B) {
    return (complex){A.x * B.x - A.y * B.y,A.x * B.y + A.y * B.x};
}
complex operator + (const complex &A,const complex &B) {
    return (complex){A.x + B.x,A.y + B.y};
}
complex operator - (const complex &A,const complex &B) {
    return (complex){A.x - B.x,A.y - B.y};
}
int rev[N];
void FFT(complex *a,int n,int xs) {
    for(int i = 0;i <= n;++i) if(rev[i] > i) swap(a[rev[i]],a[i]);

    for(int m = 2;m <= n;m <<= 1) {
        complex w1 = (complex){cos(2 * pi / m),xs * sin(2 * pi / m)};
        for(int i = 0;i < n;i += m) {
            complex w = (complex){1,0};
            for(int k = 0;k < (m >> 1);++k) {
                complex u = a[i + k];
                complex t = w * a[i + k + (m >> 1)];
                a[i + k] = u + t;
                a[i + k + (m >> 1)] = u - t;
                w = w * w1;
            }
        }
    }

}

int main() {
    int n = read();
    for(int i = 0;i < n;++i) {
        A[i].x = read();B[n - i - 1].x = read();
    }
    --n;
    int tot = 1;
    while(tot <= (n << 1)) tot <<= 1;

    for(int i = 0;i <= tot;++i) rev[i] = (rev[i >> 1] >> 1) | (i & 1 ? tot >> 1 : 0);


    FFT(A,tot,1);FFT(B,tot,1);

    for(int i = 0;i <= tot;++i) A[i] = A[i] * B[i];
    FFT(A,tot,-1);
    
    for(int i = 0;i <= n;++i) printf("%d
",(int)round(A[i + n].x / tot));



    return 0;
}