简单的 FFT 变形

「BZOJ2194」快速傅立叶之二

2015年4月29日3,8300

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

思路分析 ：

　　初看题目所要求的式子，很像卷积， f(x) * g(x) = sigma(f(x) g(t-x))  那么我们只要将 b数组变换一下即可， 另 d[i] = b[n-i-1] ， 则a[i]*b[k-i] = a[i]*b[n-1-(n+k-i-1)] = a[i]*d[n+k-1-i] ( k-1 < i < n) 这不就是一个标准的卷积了吗，fft 即可

代码示例 ：（未测试）

#include <bits/stdc++.h>
using namespace std;
#define ll long long
const int maxn = 3e5+5;
const double pi = acos(-1.0);

int n;
struct Complex{
    double x, y;
    Complex (double _x=0, double _y=0):x(_x), y(_y){}
    Complex operator -(const Complex &b)const{
        return Complex(x-b.x, y-b.y);
    }
    Complex operator +(const Complex &b)const{
        return Complex(x+b.x, y+b.y);
    }
    Complex operator *(const Complex &b)const{
        return Complex(x*b.x-y*b.y, x*b.y+y*b.x);
    }
};
Complex x1[maxn], x2[maxn];
void change(Complex y[], int len){
    for(int i = 1, j = len/2; i < len-1; i++){
        if (i < j) swap(y[i], y[j]);
        int k = len/2;
        while(j >= k){
            j -= k;
            k /= 2;
        }
        if (j < k) j += k;
    }
}
 
void fft(Complex y[], int len, int on){
    change(y, len);
    for(int h = 2; h <= len; h <<= 1){
        Complex wn(cos(-on*2*pi/h), sin(-on*2*pi/h));
        for(int j = 0; j < len; j += h){
            Complex w(1, 0);
            for(int k = j; k < j+h/2; k++){
                Complex u = y[k];
                Complex t = w*y[k+h/2];
                y[k] = u+t;
                y[k+h/2] = u-t;
                w = w*wn;
            }
        }
    }
    if (on == -1){
        for(int i = 0; i < len; i++)
            y[i].x /= len;
    }
}

int main () {
	
	cin >> n;
	for(int i = 0; i < n; i++) scanf("%lf%lf", &x1[i].x, &x2[n-i-1].x);	
	int len = 1;
	while(len < 2*n) len <<= 1;
	
	fft(x1, len, 1); fft(x2, len, 1);
	for(int i = 0; i < len; i++) x1[i] = x1[i]*x2[i];
	fft(x1, len, -1);
	
	for(int i = n-1; i < 2*n-1; i++){
		int x = (int)(x1[i].x+0.5);
		printf("%d
", x);
	}	
	return 0;
}