解题报告（二）C、（darkBZOJ 2194） 快速傅立叶之二（FFT、卷积的概念、常用变换）

繁凡出品的全新系列：解题报告系列 —— 超高质量算法题单，配套我写的超高质量题解和代码，题目难度不一定按照题号排序，我会在每道题后面加上题目难度指数（$1\sim 5$），以模板题难度 $1$ 为基准。

---

这样大家在学习算法的时候就可以执行这样的流程：

阅读我的【学习笔记】 / 【算法全家桶】学习算法 $\Rightarrow$ 阅读我的相应算法的【解题报告】获得高质量题单 $\Rightarrow$ 根据我的一句话题解的提示尝试自己解决问题 $\Rightarrow$ 点开我的详细题解链接学习巩固（好耶）
%

题单链接：【解题报告】多项式问题（多项式乘法及其各种运算）（ICPC / CCPC / NOIP / CF / AT / NC / P）超高质量题解

C、（BZOJ 2194） 快速傅立叶之二（FFT、卷积的概念、常用变换）

Wbelink

https://darkbzoj.tk/problem/2194

Problem

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

$n\le 10^5$

Solution

$$C(x)=A(x)\ast B(x)=\sum _{k=0}^{2n-2}(\sum _{k=i+j}{a}_i{b}_j)x^k$$

我们知道多项式乘法：
$$C(x)=A(x)\ast B(x)=\sum _{i=0}^n(\sum _{j=0}^i{a}_j{b}_{i-j})x^{j+i-j}=\sum _{i=0}^n(\sum _{j=0}^i{a}_j{b}_{i-j})x^i$$

也就是对于 $C(x)$ 的第 $k$ 项的系数：

$$[k]C(x)=\sum _{i=0}^k{a}_i{b}_{k-i}$$

可以理解为只有小于 $k$ 的系数才有可能给 $x^k$ 做贡献，即 $x^j\times x^{k-j}=x^k$。而大的次幂乘起来就大于 $k$ 了。

而本题中要求的式子实际上可以写成：

$$[k]C(x)=\sum _{i=k}^{n-1}{a}_i{b}_{i-k}$$

发现本题中要求的式子和标准的多项式卷积不太一样，所以考虑如何转化一下。

很明显这两个式子的唯一差别就是卷积循环的是前半部分，题中所求的式子循环的是后半部分，所以我们可以想办法把他们交换一下：我们设 $a^{\prime }$ 是 $a$ 数组的翻转，即： $a[i]=a^{\prime }[n-1-i]$ 或 $a^{\prime }[i]=a[n-1-i]$ ，二者等价。 （减一是因为 $a$ 数组是从 $0\sim n-1$）。

故原式等于：

$$[k]D(x)=\sum _{i=k}^{n-1}{a}_{n-1-i}^{\prime }{b}_{i-k}=[n-1-k]C(x)$$

因为： $i:k\sim n-1$

$j=n-1-i:n-1-k\sim n-1-(n-1)\Rightarrow 0\sim n-1-k$

$$[k]D(x)=\sum _{j=n-1-i=0}^{n-1-k}{a}_{n-1-i}^{\prime }{b}_{i-k}$$

$n-1-i+i-k=n-1-k$ 是一个常数。

这样就变成了一个标准的卷积的形式，即可使用 FFT 求解。

Code

#include <cstdio>
#include <iostream>
#include <algorithm>
#include <cstring>
#include <cmath>
#include <map>
#include <queue>
using namespace std;
typedef long long ll;
typedef int itn;
typedef pair<int, int>PII;
const int N = 2e5 + 7, mod = 1e9 + 7;

const double PI = acos(-1.0);

int n, m, RR[N], limit, L;

struct Complex
{
    double x, y;
    Complex (double x = 0, double y = 0) : x(x), y(y) { }
    Complex operator + (const Complex &t)
    {
        return Complex (x + t.x, y + t.y);
    }
    Complex operator - (const Complex &t)
    {
        return Complex (x - t.x, y - t.y);
    }
    Complex operator * (const Complex &t)
    {
        return Complex (x * t.x - y * t.y, x * t.y + y * t.x);
    }
}A[N], B[N], C[N];


void FFT(Complex *A, double type)
{
    for(int i = 0; i < limit; ++ i)
        if(i < RR[i])
            swap(A[i], A[RR[i]]);

    for(int mid = 1; mid < limit; mid <<= 1) {
        Complex wn(cos(PI / mid), type * sin(PI / mid));

        for(int len = mid << 1, pos = 0; pos < limit; pos += len) {
            Complex w(1, 0);
            for(int k = 0; k < mid; ++ k, w = w * wn) {
                Complex x = A[pos + k];
                Complex y = w * A[pos + mid + k];
                A[pos + k] = x + y;
                A[pos + mid + k] = x - y;
            }
        }
    }
    if(type == -1) {
        for(int i = 0; i <= limit; ++ i)
            A[i].x /= limit;
    }
}

//a'[i] = a[n - 1 - i];

int main()
{
    scanf("%d", &n);
    for(int i = 0; i < n; ++ i) {
        int x, y;
        scanf("%d%d", &x, &y);
        A[n - 1 - i].x = x;
        B[i].x = y;
    }
    L = 0, limit = 1;
    while(limit <= n * 2) limit <<= 1, L ++ ;
    for(int i = 0; i < limit; ++ i) {
        RR[i] = (RR[i >> 1] >> 1) | ((i & 1) << (L - 1));
    }

    FFT(A, 1), FFT(B, 1);
    for(int i = 0; i < limit; ++ i) {
        C[i] = A[i] * B[i];
    }
    FFT(C, -1);
    for(int i = n - 1; i >= 0; -- i) {
        printf("%d\n", (int)(C[i].x + 0.5));
    }
    return 0;
}