bzoj2194 快速傅立叶之二

2194: 快速傅立叶之二

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

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

分析：很显然，这道题要用FFT来做，但是题目所给的式子并不是一个标准的卷积的形式，需要进行变形.

　　　一般而言，能用FFT求的式子c[k] = a[i] * b[j], i + j是一个只与k有关的式子，其它的都是常数项.  对于这道题，把bi变成b_n-i-1

那么ck = ai*b_i - k = ai * b_n - (i - k) - 1 = ai * b_n -i - 1 + k.  a和b的下标加起来恰好是n + k,n是常数，满足做fft的要求.

　　　最后c数组的值要整体向后移k位. 

#include <cstdio>
#include <cmath>
#include <cstring>
#include <iostream>
#include <algorithm>

using namespace std;

const int maxn = (1 << 18) + 5;
const double pai = acos(-1.0);
int n,a[maxn],b[maxn],len;

struct node
{
    double real, imag;
    node(double real = 0.0, double imag = 0.0)
    {
        this->real = real, this->imag = imag;
    }
    node operator - (const node&elem) const
    {
        return node(this->real - elem.real, this->imag - elem.imag);
    }
    node operator + (const node&elem) const
    {
        return node(this->real + elem.real, this->imag + elem.imag);
    }
    node operator * (const node&elem) const
    {
        return node(this->real * elem.real - this->imag * elem.imag, this->real * elem.imag + this->imag * elem.real);
    }
    void set(double real = 0.0, double imag = 0.0)
    {
        this->real = real, this->imag = imag;
    }
} A[maxn],B[maxn];

void pre()
{
    for (int i = 0; i < n; i++)
        scanf("%d%d",&a[i],&b[i]);
    len = 1;
    while (len < (n << 1))
        len <<= 1;
    for (int i = 0; i < n; i++)
        A[i].set(a[i],0),B[i].set(b[n - i - 1],0);
    for (int i = n; i < len; i++)
        A[i].set(),B[i].set();
}

void Swap(node &a,node &b)
{
    node temp = a;
    a = b;
    b = temp;
}

void zhuan(node *y)
{
    for (int i = 1,j = len >> 1,k; i < len - 1; i++)
    {
        if (i < j)
            Swap(y[i],y[j]);
        k = len >> 1;
        while (j >= k)
        {
            j -= k;
            k >>= 1;
        }
        if (j < k)
            j += k;
    }
}

void FFT(node *y,int op)
{
    zhuan(y);
    for (int h = 2; h <= len; h <<= 1)
    {
        node temp(cos(op * pai * 2 / h),sin(op * pai * 2 / h));
        for (int i = 0; i < len; i += h)
        {
            node W(1,0);
            for (int j = i; j < i + h / 2; j++)
            {
                node u = y[j];
                node v = W * y[j + h / 2];
                y[j] = u + v;
                y[j + h / 2] = u - v;
                W = W * temp;
            }
        }
    }
    if (op == -1)
        for (int i = 0; i < len; i++)
            y[i].real /= len;
}

void solve(node *A,node *B)
{
    FFT(A,1);
    FFT(B,1);
    for (int i = 0; i < len; i++)
        A[i] = A[i] * B[i];
    FFT(A,-1);
}

int main()
{
    scanf("%d",&n);
    pre();
    solve(A,B);
    for (int i = n - 1; i < 2 * n - 1; i++)
        printf("%d\n",(int)(A[i].real + 0.5));

    return 0;
}