FFT [TPLY]

题目链接
https://www.luogu.org/problemnew/show/1919
https://www.luogu.org/problemnew/show/3803

资料推荐

orz 大佬博客
http://www.cnblogs.com/cjoieryl/p/8206721.html (orz YL 大佬)
http://blog.csdn.net/iamzky/article/details/22712347 (超级易懂)

知识点
复数:
https://baike.baidu.com/item/%E5%A4%8D%E6%95%B0/254365?fr=aladdin
单位根 (单位复数根)
https://baike.baidu.com/item/%E5%8D%95%E4%BD%8D%E6%A0%B9
多项式表示
http://blog.csdn.net/acdreamers/article/details/39005227
卷积
http://blog.csdn.net/bitcarmanlee/article/details/54729807

书籍推荐
ACM/ICPC 算法基础训练教程 7.4 快速傅里叶变换
算法导论 30 章多项式与快速傅里叶变换

觉得看不懂很正常, 知识点很抽象需要逐步理解
而对于本文, 本文作用为对两位大佬内容解读并且加入自己看法
读者可以先花时间阅读以上推荐的文字, 产生一定理解再读本文. 由于作者很弱, 可能会产生错误, 还请大神帮忙指出.

正文

Part1. 初识 FFT

作用: 求多项式乘法的系数 (就是初中学的多项式)
FFT 多项式乘法与普通形式有差异

一般多项式乘多项式法则:
多项式与多项式相乘，先用一个多项式的每一项与另一个多项式的每一项相乘，再把所得的积相加。举一个例子由多项式乘多项式法则可以得到 (a+b)(c+d)=a(c+d)+b(c+d)=ac+ad+bc+bd

FFT‘ 绕圈’ 法
FFT 处理多项式相乘, 先从多项式最普通的表达形式–-系数表达形式转化成点值表达形式再进行点值乘法. 最后, 再将点值乘法的结果转回系数表达形式.
而从-系数表达形式转化成点值表达形式的过程, 我们叫做 DFT, 即离散傅里叶变换.
把点值乘法的结果转回系数表达形式的步骤, 我们称之为离散傅里叶变换的逆运算, 也叫插值.

下面再对一些可能不理解的概念进行解释
１～对于点值表达式我的理解 :
点值表达形式就是用若干函数上的点的坐标来表示该函数
比如 f(x)=$x^2$-x+2={ ( 0 , 2 ) , ( 1 , 2 ) , ( 2 , 4 ) }
但是这些点的数量一定不小于系数的个数 n
这是一个常识
就像求解 n 元方程一定要有 n 组方程才能得解
２～点值表达式的计算:
假设我们现在有用点乘表示的多项式
A(x)={(0,1),(1,2),(2,7)}
B(x)={(0,2),(1,2),(2,4)}
[P.S.: 点乘表达式能运算的条件是两个点要有相同的 x 值]
C(x)={(0,2×1),(1,2×2),(2,7×4)}={(0,2),(1,4),(2,28)}

FFT 流程图如下:(借鉴他人博客)

part2. 复数

0. 为什么要学习复数?
YL 大佬原话:” 学习复数, 所有这些七七八八的定理, 都是为了 FFT 的奇妙变换!”

所以, 为了 FFT, 好好学习复数吧!(再一次 ORZ YL 大佬)

1. 复数的定义: 我们把形如 a+bi（a,b 均为实数）的数称为复数，其中 a 称为实部，b 称为虚部，i 称为虚数单位,i 的值为 $\sqrt{-1}$
2. 单位根的定义: 数学上，n 次单位根是 n 次幂为 1 的复数。
单位根的概念不理解很正常, 下面举个例子 (来自百度百科)
单位的一次根有一个：1
单位的二次根有两个：+1 和-1。
单位的三次根是:
{$ 1 $, $\frac{-1+\sqrt{3}i}{2}$,$\frac{-1-\sqrt{3}i}{2}$}
咱们证明一下
$(\frac{-1+\sqrt{3}i}{2})^3$=$(\frac{-1+\sqrt{3}i}{2})^2$×$(\frac{-1+\sqrt{3}i}{2})$
………………..=$(\frac{1+3i^2-2\sqrt{3}i}{4})$×$(\frac{-1+\sqrt{3}i}{2})$
………………..=$(\frac{-2-2\sqrt{3}i}{4})$×$(\frac{-1+\sqrt{3}i}{2})$
………………..=$(\frac{(-2)×(1+\sqrt{3}i)×(\sqrt{3}i-1)}{8})$
………………..=-$(\frac{1}{4})$×($3i^2-1$)
………………..=$-(\frac{1}{4})$*(-4)
………………..=1
所以 $\frac{-1+\sqrt{3}i}{2}$是单位的三次根
单位的四次根是 {1,+i,-1,-i}

希望你已经理解了什么是单位复数根 (单位根), 咱们继续

强烈建议学习一下著名的欧拉公式，后面会用到。
单位复数根定义：对于 w$_n$=$e^{2πi/n}$我们称 w$_n$为 n 次单位根 (前面讲了 n 次单位跟哦)。

下面有一个定理需要记住 (YL 大佬的忠告)
PS 由于太弱不会证明,YL 大佬的博客上由详细证明可以参考
(w 上标表示指数)
消去引理

$$\omega^{dk} _ {dn}=\omega^{k}_{n}$$

这个引理非常非常重要哦，你后面一定一定会见到！
非常，非常重要！记住！

这些复杂的数学公式可能会使你觉得枯燥无味, 但为后面的运算奠定了重要的基础, 了解一下下啦!

Part3.DFT 与-DFT 进阶

在本章内于括号里面的数为什么是 $\omega^{k}{n}$而不是 x, 你要到 FFT 才能明白. 但本章内容中的 $\omega^{k}{n}$你可以当成 x 看 (下一章就不行了)

DFT
对于 k=0~n-1，定义：
$$y _ k=A(\omega^{k} _ {n}) = \Sigma^{n-1} _ {j=0} a _ j(\omega^{k} _ {j})^j$$
对于这个式子, 你这么看 $y _ k=A(x)=\Sigma^{n-1} _ {j=0} a _ j(x)^j$
我们知道 DFT 求的是点乘
点乘就是坐标
平时你是怎么求坐标的?
选定一个 x 值, 把 x 值带到里面求出 y
上面这个式子就是这么做的啦 (还不懂就再看看多项式的系数表达式)

所以我们继续分析, 对于每一个 x, 我们都要 O(n) 地将值带到多项式里求值, 我们最少要 n 个坐标, 所以总复杂度是 O(n * n)

最后我们将得到的 y 称为 a 的离散傅里叶变换，记作 $y=DFT_n(a)$ (这里的 y，a 指的是所有的 y_k,a_k(k 有 n 种, 也就是要取 k 个))

这个定义先记着, 还要注意的是, 这一步的作用是将系数表达形式变成点乘表达形式 $O(n^2)$.
再用之前的点乘算法算出点乘的结果 O(n)(Part1. 讲了点乘, 复杂度证明很显然).
最后用下面的离散傅里叶逆变换(我叫它-DFT), 再换回去 $O(n^2)$所以总复杂度为 $O(n^2)$

-DFT
我这里直接摆出结果，有兴趣的可以自己证明一下
对于 $y_k = \Sigma^{n-1} _ {i=0}a _ i(x)^i$
有 $a_k = \frac{1}{n}\Sigma^{n-1} _ {i=0}y _ i(x)^i$
再对比 DFT 公式
$y _ k = \Sigma^{n-1} _ {j=0} a _ j(x)^j$
可以发现,DFT 和-DFT 的差别很小, 基本上只是 y,a 调换了一下罢了 (也有不同的…)

Part4 FFT 水落石出

终于来到了最后, 前面所有令人头大的背景知识, 都是为了这一章而服务.
所以这一章是至关重要的.
如果前面还没有完全理解, 请重新巩固扎实.

网上都说,FFT 的复杂度为 O(nlgn)
那是怎么来的呢?
答案是: 神奇的分治
先把 A(x) 这个多项式拆分成奇数项与偶数项
$$A^{[0]}(x) = a _ 0 + a _ 2x + a_4x^2 + … + a _ {n-2}x^{\frac{n}{2} – 1}$$
$$A^{[1]}(x) = a _ 1 + a _ 3x + a_5x^2 + … + a _ {n-1}x^{\frac{n}{2} – 1}$$
$$A(x) = A^{[0]}(x^2) + xA^{[1]}(x^2) $$
这个是正确的, 希望读者拿起草稿纸自己演算一遍
咱们把这个多项式拆分成这样子以后
复数闪亮登场
开始用复数的公式了哦
咱们不妨把 $\omega^{k}$带入到 x 中去
$$A(\omega^k _ n)=A^{[0]}((\omega^k _ n)^2) + \omega^k _ n A^{[1]}((\omega^k _ n)^2)$$
然后, 通过之前提到过的单位复数根消去引理我们可以得到 （说了消去引理非常重要）
咱们用消去引理把 $(\omega^k _ n)^2$化一下
$（（\omega^k _ n)^2$=$\omega^{2k} _ n$

………….=$\omega^{\frac{2k}{2}} _ {\frac{n}{2}}$

………….=$\omega^{k} _ {\frac{n}{2}}$
又因为 $$有名的欧拉公式\omega^{\frac{n}{2}} _ {n}=\omega _ {2}=e^{k\pi i}= cos k\pi + i sin k\pi = -1(有名的欧拉公式)$$
所以

$$\omega^{k+\frac{n}{2}} _ n=\omega^{k} _ n*\omega^{\frac{n}{2}} _ n=-\omega^{k} _ n$$

所以
$$A(\omega^{k+\frac{n}{2}} _ n)=A^{[0]}((\omega^{k+\frac{n}{2}} _ n)^2) + \omega^{k+\frac{n}{2}} _ nA^{[1]}((\omega^{k+\frac{n}{2}} _ n)^2)$$
$$=A^{[0]}((-\omega^k _ {n})^2)-\omega^k _ {n}A^{[1]}((-\omega^k _ {n})^2)$$
$$=A^{[0]}(\omega^{2k} _ {n})-\omega^k _ {n}A^{[1]}(\omega^{2k} _ {n})$$
$$再一次用到咱们的消去引理 再一次用到咱们的消去引理 $$
$$。。。（式）=A^{[0]}(\omega^k _ {\frac{n}{2}})-\omega^k _ {n}A^{[1]}(\omega^k _ {\frac{n}{2}})。。。（式 1）$$
然而我们同时也可以暴力把原式化一下：

$$A(\omega^k _ n)=A^{[0]}((\omega^k _ n)^2) + \omega^k _ n A^{[1]}((\omega^k _ n)^2)$$
$$=A^{[0]}(\omega^{2k} _ n) + \omega^k _ n A^{[1]}(\omega^{2k} _ n)$$
$$又因为消去引理 又因为消去引理 $$
$$。。。（式）=A^{[0]}(\omega^{k} _ {\frac{n}{2}}) + \omega^k _ n A^{[1]}(\omega^k _ {\frac{n}{2}})。。。（式 2）$$

蝴蝶操作

所以我们把我们推了这么久的（式 1）和（式 2）拿出来并排看
$$A(\omega^k _ n)=A^{[0]}(\omega^{k} _ {\frac{n}{2}}) + \omega^k _ n A^{[1]}(\omega^k _ {\frac{n}{2}}) $$
$$A(\omega^{k+\frac{n}{2}} _ n)=A^{[0]}(\omega^k _ {\frac{n}{2}})-\omega^k _ {n}A^{[1]}(\omega^k _ {\frac{n}{2}})$$
然后我们发现
用 $$和和 A^{[0]}(\omega^{k} _ {\frac{n}{2}}) 和 A^{[1]}(\omega^k _ {\frac{n}{2}}) 和\omega^k _ {n}$$
可以同时算出 $$和 A(\omega^k _ n) 和 A(\omega^{k+\frac{n}{2}} _ n)$$
计算量便减小了一半
我们管这叫蝴蝶操作

蝴蝶操作流程如下

然后我们还可以发现
蝴蝶操作无非就是把奇数项和偶数项合并算出新项再往上继续往上合并直到我们得到最终的式子

也就是这样

是可以递归着做
不断分离奇数项偶数项最后再回溯合并
这是由上而下

那么能不能由下而上呢？
我们看最下层叶子节点数字的编号，通过找规律我们发现每个数和他二进制相反的位置互换
举个例子
4 -> 100
1 -> 001
所以 4 和 1 交换
（注意标号从 0 开始）

接下来。。。

我承认我蝴蝶操作讲的不清楚
由于知识有限只能讲到这里
网上关于蝴蝶操作的好的文章有很多。
大家可以去看一看
下面是一些例题
（链接再最上面）
可以把它当成板子背一背

代码

//P3803【模板】多项式乘法（FFT）
#include <algorithm>
#include <iostream>
#include <cstring>
#include <cstdio>
#include <queue>
#include <cmath>

#define rg register int
#define ll long long
#define db double
#define il inline

#define INF 2147483647
#define SZ 10000001

using namespace std;

const double pi=acos(-1);

struct Complex
{
    db r,i;
    Complex(){}
    Complex(db a,db b):r(a),i(b) {}
    Complex operator+(const Complex B)const
    {return Complex(r+B.r,i+B.i);}
    Complex operator-(const Complex B)const
    {return Complex(r-B.r,i-B.i);}
    Complex operator*(const Complex B)const
    {return Complex(r*B.r-i*B.i,r*B.i+i*B.r);}
};
Complex X,Y,a[SZ],b[SZ];

int r[SZ],n,m,l;
il void FFT(Complex *a,int x)
{
    for(rg i=0;i<n;++i) if(i<r[i]) swap(a[i],a[r[i]]);
    for(rg i=1;i<n;i<<=1)
    {
        Complex wn(cos(pi/i),x*sin(pi/i));
        for(rg j=0;j<n;j+=(i<<1))
        {
            Complex w(1,0);
            for(rg k=0;k<i;++k,w=w*wn)
            {
                X=a[j+k],Y=a[i+j+k]*w;
                a[j+k]=X+Y,a[i+j+k]=X-Y;
            }
        }
    }
    if(x==-1) for(rg i=0;i<n;++i) a[i].r=a[i].r/n;
}

int main()
{
    scanf("%d%d",&n,&m);
    for(rg i=0;i<=n;++i) scanf("%lf",&a[i].r);
    for(rg i=0;i<=m;++i) scanf("%lf",&b[i].r);
    m+=n;for(n=1;n<=m;n<<=1) ++l;
    for(rg i=0;i<n;++i) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
    FFT(a,1);FFT(b,1);
    for(rg i=0;i<=n;++i) a[i]=a[i]*b[i];
    FFT(a,-1);
    for(rg i=0;i<=m;++i) printf("%d ",(int)(a[i].r+0.5));
    return 0;
}

P1919【模板】A*B Problem 升级版（FFT 快速傅里叶）
#include <algorithm>
#include <iostream>
#include <cstring>
#include <cstdio>
#include <queue>
#include <cmath>

#define rg register int
#define ll long long
#define il inline
#define ldb long double

#define INF 2147483647
#define SZ 1000001

using namespace std;
const ldb pi=acos(-1);

struct Complex
{
    ldb r,i;
    Complex(){}
    Complex(ldb a,ldb b):r(a),i(b){}
    il Complex operator+(const Complex B)const
    {return Complex(r+B.r,i+B.i);}
    il Complex operator-(const Complex B)const
    {return Complex(r-B.r,i-B.i);}
    il Complex operator*(const Complex B)const
    {return Complex(r*B.r-i*B.i,r*B.i+i*B.r);}
};
Complex X,Y,a[SZ],b[SZ];
int n,m,l,r[SZ];
il void FFT(Complex *a,int x)
{
    for(rg i=0;i<n;++i) if(i<r[i]) swap(a[i],a[r[i]]);
    for(rg i=1;i<n;i<<=1)
    {
        Complex wn(cos(pi/i),x*sin(pi/i));
        for(rg j=0;j<n;j+=(i<<1))
        {
            Complex w(1,0);
            for(rg k=0;k<i;++k,w=w*wn)
            {
                X=a[j+k],Y=a[i+j+k]*w;
                a[j+k]=X+Y,a[i+j+k]=X-Y;
            }
        }
    }
    if(x==-1) for(rg i=0;i<n;++i) a[i].r=a[i].r/n;
}

char ch1[SZ],ch2[SZ];
int res[SZ];
int main()
{
    scanf("%d",&n);
    cin>>ch1;
    cin>>ch2;
    for(rg i=0;i<n;++i)
    {
        a[i].r=ch1[n-i-1]-'0';
        b[i].r=ch2[n-i-1]-'0';
    }
    m=n+n-2;for(n=1;n<=m;n<<=1) ++l;
    for(rg i=0;i<n;++i) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
    FFT(a,1);FFT(b,1);
    for(rg i=0;i<=n;++i) a[i]=a[i]*b[i];
    FFT(a,-1);
    for(rg i=0;i<=m;++i) res[i]=(int)(a[i].r+0.5);
    for(rg i=0;i<=m;++i) res[i+1]+=res[i]/10,res[i]%=10;
    rg top=m+1;
    while(res[top]>9)
    {
        res[top+1]=res[top]/10;
        res[top]%=10;
        ++top;
    }
    while(!res[top]) --top;
    for(rg i=top;i>=0;--i)
     printf("%d",res[i]);
    return 0;
}

几点注意
pi 最好不要手写（当然如果你能背到小数点后十几位我也不栏你）
不然很可能出精度问题
最好手写 Complex 速度快很多
注意一些从 0 开始的地方
加油哦!

END

$\times$