【BZOJ2194】快速傅立叶之二 [FFT]

快速傅立叶之二

Time Limit: 10 Sec  Memory Limit: 259 MB
[Submit][Status][Discuss]

Description

　　请计算C[k]=sigma(a[i]*b[i-k]) 其中 k < = i < n。 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

HINT

　　 n < = 10 ^ 5

Solution

　　显然是运用FFT，看到题目里 b 的下标为 i-k，于是乎我们就要想一个办法，把它弄成卷积的形式。
　　然后翻转一下，下标就变成了(n-1)-(i-k)。那 Ans[n-1+k]=Σa[i]*b[(n-1)-(i-k)] 啦。

　　至于卷积这个东西，当然是让我们膜拜AireenYe啦！(←_←戳名字有惊喜)

Code

#include<iostream>    
#include<string>    
#include<algorithm>    
#include<cstdio>    
#include<cstring>    
#include<cstdlib>
#include<cmath>
#include<map>
using namespace std;  
typedef long long s64;
   
const int ONE = 5e5+5;
const double Pi = acos(-1.0);
  
int n, m;
int tn, tl;
int bitRev[ONE];
 
struct Complex
{
        double r, i;
        Complex() {}
        Complex(double _r, double _i)
            : r(_r), i(_i) {}
        friend Complex operator +(Complex a, Complex b)
        {
            return Complex(a.r + b.r, a.i + b.i);
        }
        friend Complex operator -(Complex a, Complex b)
        {
            return Complex(a.r - b.r, a.i - b.i);
        }
        friend Complex operator *(Complex a, Complex b)
        {
            return Complex(a.r*b.r - a.i*b.i, a.r*b.i + a.i*b.r);
        }
}a[ONE], b[ONE];
 
int get()
{    
        int res=1,Q=1;char c;    
        while( (c=getchar())<48 || c>57 ) 
        if(c=='-')Q=-1; 
        res=c-48;     
        while( (c=getchar())>=48 && c<=57 )    
        res=res*10+c-48;    
        return res*Q;
}
 
 
void FFT_init(int n)
{
        tn = 1, tl = 0;
        while(tn < n)
            tn <<= 1, tl ++;
        for(int i = 0; i < tn; i++)
        {
            int l = bitRev[i >> 1] >> 1;
            int r = (i & 1) << (tl - 1);
            bitRev[i] = l | r;
        }
}
 
void FFT(Complex *a, int rev)
{
        for(int i = 0; i < tn; i++)
        {
            int r = bitRev[i];
            if(i < r) swap(a[i], a[r]);
        }
         
        for(int k = 1; k < tn; k <<= 1)
        {
            Complex wn(cos(Pi / k), rev * sin(Pi / k));
            for(int s = 0; s < tn; s += k<<1)
            {
                Complex w(1, 0);
                for(int i = s; i < s + k; i++)
                {
                    Complex f1 = a[i], f2 = w * a[i + k];
                    a[i] = f1 + f2;
                    a[i + k] = f1 - f2;
                    w = w * wn; 
                }
            }
        }
}
 
int main()
{
        n = get();
        for(int i=0; i<n; i++)
            a[i].r = get(), b[n-1-i].r = get();
         
        FFT_init(n + n + 1);
         
        FFT(a, 1);  FFT(b, 1);
        for(int i=0; i<tn; i++)
            a[i] = a[i] * b[i];
        FFT(a, -1);
         
        for(int i=n-1; i<n+n-1; i++)
            printf("%d\n", (int)(a[i].r / tn + 0.5));
         
}