FFT(MTT)常数优化

最近闲着无聊研究了下$FFT$的常数优化，大概就是各种$3$次变$2or1.5$次之类的，不过没见过啥题卡这个的吧。

关于$FFT$可以看这里：浅谈FFT&NTT。

关于复数

设$x=a+bi$，其中$i$是虚数单位，那么我们用$\bar x$表示$x$的共轭复数，即$\bar x=a-bi$。

共轭复数有一个这样的性质：

$$\overline{ab}=\bar a \cdot \bar b$$

证明展开就好了，这个是下面优化的关键。

设$\omega_n$为$n$阶单位根，则$\overline{\omega _n^{x}}=\omega_{n}^{-x}$。

idft变dft

设$f(x)=\sum_{i=0}^{n-1}a_ix^i$，注意到：

$$n\cdot[j]{\rm idft}(f)=\sum_{i=0}^{n-1}\omega_{n}^{-ij}a_i=a_0+\sum_{i=1}^{n-1}\omega_{n}^{ij}a_{n-i}$$

也就是说我们如果先dft(a)，再进行一次std::reverse(a+1,a+n)，再除以$n$，就完成了一次$idft$。

多项式乘法优化 1

给出多项式$f(x)=\sum_{i=0}^{n}a_ix^i,g(x)=\sum_{i=0}^{m}b_ix^i$，求其卷积。

这里最开始介绍一种非常简洁的优化方法，构造多项式$h(x)$：

$$h(x)=f(x)+ig(x)$$ $$h^2(x)=f^2(x)-g^2(x)+2if(x)g(x)$$

那么我们只需要取$h^2(x)$的虚部除以$2$就是答案，这只需要做两次$FFT$。

多项式乘法优化 2

这个和上面的关联不大，设$X_i$表示多项式$F(x)$$dft$之后的系数，$a_i$表示$dft$之前的系数，设$F(x)$为$n$项的多项式，且$n=2^k$，注意到：

$$X_i=\sum_{j=0}^{n-1}a_j\omega_{n}^{ij},X_{n-i}=\sum_{j=0}^{n-1}a_j\omega_{n}^{-ij}$$

即：$X_i=\overline{X_{n-i}}$。

这实质上是因为$F$没有虚部的原因，我们换一个有虚部的多项式试试：

$$X_i=\sum_{j=0}^{n-1}(a_j+ib_j)\omega_{n}^{ij}\\ X_{n-i}=\sum_{j=0}^{n-1}(a_j+ib_j)\omega_{n}^{-ij}\\ \overline{X_{n-i}}=\sum_{j=0}^{n-1}(a_j-ib_j)\omega_{n}^{ij}\\$$

等等，我们发现第一个式子和第三个式子很像，两式相加减可以得到：

$$X_i+\overline{X_{n-i}}=2\sum_{j=0}^{n-1}a_j\omega_{n}^{ij}\\ X_i-\overline{X_{n-i}}=2i\sum_{j=0}^{n-1}b_j\omega_{n}^{ij}$$

注意到等式右边就是$a$ $dft$完之后的结果，那么对于多项式$F(x),G(x)$，我们可以构造一个函数然后$dft$一次，然后$O(n)$得到两个多项式$dft$之后的结果，总共只用了一次$FFT$。

---

当然这个玩意也可以这样用：假设我们现在想求$dft(F(x))$，我们把$F(x)$奇偶分类，构造多项式：

$$g(x)=\sum_{i=0}^{n/2-1}(a_{2i}+ia_{2i+1})x^i$$

然后相当于是$0.5$次$FFT$来完成这个事，设$dft(g(x))$每一项为$X_i$，$dft(F(x))$每一项为$Y_i$，那么推一下可以得到：

$$Y_i=\frac{X_i+\overline{X_{n/2-i}}}{2}-2\omega_{n}^i(X_i-\overline{X_{n/2-i}})$$

注意这里只有$i\in [0,n/2)$的值，$Y_{n/2}$特殊处理一下，后面的可以通过前面得到。

MTT常数优化

$\rm MTT$就是拆系数$\rm FFT$，设多项式$s(x),t(x)$，我们要算$s(x)t(x)$，模数任意。

我们拆系数，设拆完了之后是$s(x)=a(x)+b(x)\cdot p,t(x)=c(x)+d(x)\cdot p$。

构造$F(x)=a(x)+i\cdot b(x)$，$G(x)=c(x)+i\cdot d(x)$。

那么有：

$$\begin{align} &F(\omega_n^j)=\sum_{i=0}^{n-1}(a_i+ib_i)\omega_n^{ij}\\ &F(\omega_n^{-j})=\sum_{i=0}^{n-1}(a_i+ib_i)\omega_n^{-ij}\\ &\overline{F(\omega_n^{-j})}=\sum_{i=0}^{n-1}(a_i-ib_i)\omega_n^{ij}\\ \end{align}$$

那么相加减可得$a(x),b(x)$的$dft$。

令$h(x)={\rm dft}(a(x))\cdot {\rm dft}(G(x))={\rm dft}(a(x)\cdot G(x))={\rm dft}(a(x)c(x)+i\cdot a(x)d(x))$。

那么我们$idft$一次$h(x)$就可以得到$a(x)c(x),a(x)d(x)$。

同理可以得到$b(x)c(x),b(x)d(x)$，一共$4$次$dft$。

代码长这样：

void mul(int *r,int *s,int *t,int len) {
    for(N=1,bit=0;N<len;N<<=1,bit++);
    for(int i=1;i<N;i++) pos[i]=pos[i>>1]>>1|((i&1)<<(bit-1));
    for(int i=0;i<N;i++) g[0][i]=cp(r[i]&all,r[i]>>15),g[1][i]=cp(s[i]&all,s[i]>>15);
    fft(g[0]),fft(g[1]);
    for(int i=0;i<N;i++) {
        int j=(N-i)&(N-1);
        g[2][j]=(g[0][i]+conj(g[0][j]))*cp(0.5,0)*g[1][i];
        g[3][j]=(g[0][i]-conj(g[0][j]))*cp(0,-0.5)*g[1][i];
    }fft(g[2]),fft(g[3]);
    for(int i=0;i<N;i++) g[2][i]=g[2][i]/N,g[3][i]=g[3][i]/N;
    for(int i=0;i<N;i++) {
        ll pp=g[2][i].r+0.5,x=g[2][i].i+0.5,y=g[3][i].r+0.5,z=g[3][i].i+0.5;
        t[i]=(pp%p+(((x+y)%p)<<15)+((z%p)<<30))%p;
    }
}