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;
}

转载于:https://www.cnblogs.com/zbtrs/p/8639793.html