bzoj2194 快速傅里叶之二 [快速傅里叶变换]

in 题目 with 0 comment

Description

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


Input

第一行一个整数N,接下来N行,第i+2..i+N-1行,每行两个数,依次表示a[i],b[i] (0 < = i < N)。


Output

输出N行,每行一个整数,第i行输出C[i-1]。


Sample Input

5
3 1
2 4
1 1
2 4
1 4


Sample Output

24
12
10
6
1


Solution

分析题目,将题中给出的式子简单变换,将b数组反转:
$$ C_k=\sum_{i=k}^n(a_i\cdot b_{i-k}) \rightarrow C_k=\sum_{i=k}^n(a_i\cdot b_{n-i-k})$$
将c数组反转:
$$ C_{n-k}=\sum_{i=k}^n(a_i\cdot b_{k-i}) $$
即转化成了卷积的形式,直接FFT求解即可


Code

#include<bits/stdc++.h> 
#define pi acos(-1.0)
#define maxn 300010
//#define DEBUG
using namespace std;
int n;
complex<double> a[maxn], b[maxn];

inline int read()
{
    char ch;
    int read = 0, sign = 1;
    do
        ch = getchar();
    while ((ch<'0'||ch>'9')&&ch!='-');
    if (ch == '-') sign = -1, ch = getchar();
    while (ch >= '0' && ch <= '9')
    {
        read = read * 10 + ch - '0';
        ch = getchar();
    }
    return read*sign;
}

inline int Power2(int x)
{
    int x0;
    for (x0 = 1; x0 < x; x0 <<= 1);
    return x0;
}

inline int lg(int n)
{
    int l = 0;
    if (n == 0) return l;
    for (int x0 = 1; x0 <= n; x0 <<= 1) l++;
    return l;
}

inline int rev(int x, int n)
{
    int out = 0;
    while (n--) out = (out + (x & 1)) << 1, x >>= 1;
    return out>>1;
}

void FFT(complex<double> a[],int n, int flag)
{
    complex<double> A[n+1];
    for (int i = 0, l = lg(n - 1); i < n; ++i) A[rev(i, l)] = a[i];
#ifdef DEBUG
    int l=lg(n-1);
    cerr<<"l="<<l<<endl;
    for(int i=0;i<n;++i) cerr<<rev(i,l)<<" ";
    cerr<<endl;
#endif 
    for (int i = 2; i <= n; i <<= 1)
    {
        complex<double> dw(cos(2*pi/i),sin(flag*2*pi/i));
        for (int j = 0; j < n; j += i)
        {
            complex<double> w(1.0, 0);
            for (int k = 0; k < (i >> 1); k++, w *= dw)
            {
                complex<double> u = A[j + k];
                complex<double> t = w*A[j + k + (i >> 1)];
                A[j + k] = u + t;
                A[j + k + (i >> 1)] = u - t;
            }
        } 
        if (flag == -1)
            for (int i = 0; i < n; i++) a[i].real() = int(A[i].real() / n + 0.5);
        else
            for (int i = 0; i < n; i++) a[i] = A[i];
    }
}

int main()
{
    n = read();
    for (int i = 0; i < n; ++i)
        a[n-i-1] = read(), b[i] = read();
#ifdef DEBUG
    for(int i=0;i<n;++i) cerr<<a[i].real()<<" ";
    cerr<<endl;
    for(int i=0;i<n;++i) cerr<<b[i].real()<<" "; 
    cerr<<endl;
#endif 
    int length = Power2(n);
#ifdef DEBUG
    cerr<<"length="<<length<<endl;
#endif
    FFT(a, length, 1);
    FFT(b, length, 1);
    for (int i = 0; i < length; ++i) a[i] *= b[i];
    FFT(a, length, -1);
    for (int i = n-1; i >= 0; --i) cout << a[i].real() << endl;
    return 0;
}

扫描二维码,在手机上阅读!
Responses