FFT

· · 个人记录

FFT

简单记一下FFT，方便自己理解

两个多项式做乘法，暴力的想法是系数依次考虑，复杂度O(n^2)

然而如果我们观察出来了一个性质，不妨设我们现在想求解：

F(x)*G(x)

那么显然有乘出来的多项式P(x)中第k项的系数c_k满足：

c_k=\sum_{i+j=k}a_i*b_j

这是一个加法卷积。

所以FFT说他用来做加法卷积也可以，说他用来搞多项式乘法也可以

接下来是FFT的核心思想：

将两个多项式转成点值表达。

将点值两两乘起来。

将点值换成系数表达。

前面那一部分我们叫做DFT，后面那一部分叫做IDFT（逆）

前方高能：

当我们在做求点值表达的时候，会发现复杂度很高，原因是因为我们没有代入合适的值。

所以我们考虑代入合适的值进去。

单位根\omega_{n}^k表示单位根

他是定义在复平面的单位圆上的。

n$次单位根会将此圆$n$等分，然后取k块作为$\omega_{n}^k

显然\omega_n^{0}=1

以及n次单位根的n次方=1

复数乘法有个好玩的，就是它的乘法遵循模长相乘幅角相加。

因为单位根定义在单位圆上所以模长是1，而他们的幅角就加起来好了。

我们尝试给一个n-1次多项式代入n个单位根进去（以下n均视为2的整数次幂）

会由于下列两条性质而发生翻天覆地的变化：

第一： $$\omega_{n}^{2k} = \omega_{n/2}^{k}$$ 这个还好吧，显然有$\omega_{n}^k=(\omega_n^1)^k

且有：\omega_{n/2}^1=\omega_{n}^2

第二：

\omega_{n}^{k+n/2}=-\omega_{n}^k

原因很简单

这里我们定义的负数的关于原点对称了的。

且因为：\omega_{n}^{n/2}=-1

所以\omega_{n}^{n/2}*\omega_{n}^k=-\omega_{n}^k

$$\sum_{j=0}^{n-1}(\omega_{n}^k)^j=0/n$$ 简单证明一下： 如果$k=0$，那么显然答案为$n

若k\ne0，那么我们可以记s_n=\sum_{j=0}^{n-1}(\omega_{n}^k)^j

于是有：

s_n*\omega_{n}^k=s_n - 1+(\omega_{n}^k)^{n-1+1}

稍微化简一下得：s_n=\dfrac{(\omega_{n}^k)^n-1}{\omega_{n}^k-1}

显然有分母不为0，而分子=0

于是我们得到了求和性质。

前方真高能！

我们有了折半性质后怎么做DFT?

对于一个多项式F(x)

不妨设它为F(x)=a_0+a_1x+a_2x^2+...a_{n-1}x^{n-1}

不妨设：

Fl(x)=a_0+a_2x+a_4x^2+...a_{n-2}x^{\frac{n}{2}-1} Fr(x)=a_1+a_3x+a_5x^2+...a_{n-1}x^{\frac{n}{2}-1}

于是有：

F(x)=Fl(x^2)+xFr(x^2)

接下来我们代入\omega_n^k

得到：

F(\omega_n^k)=Fl((\omega_n^k)^2)+\omega_n^k*Fr((\omega_n^k)^2)

于是：

F(\omega_n^k)=Fl(\omega_{n/2}^k)+\omega_n^k*Fr(\omega_{n/2}^k)

这里假设k< n/2

所以我们想要求出k取[0,n/2)的所有值的单位根的点值都只要递归处理Fl,Fr即可。

接下来我们考虑k\ge n/2

不妨设它为k+n/2(k<n/2)

于是有：

F(\omega_n^{k+n/2})=Fl((\omega_n^{k+n/2})^2)+\omega_n^{k+n/2}*Fr((\omega_n^{k+n/2})^2)

得到：

F(\omega_n^{k+n/2})=Fl((-\omega_n^{k})^2)-\omega_n^{k}*Fr((-\omega_n^{k})^2)

即：

F(\omega_n^k)=Fl(\omega_{n/2}^k)-\omega_n^k*Fr(\omega_{n/2}^k)

你会惊人的发现这个东东和前面的一样只是一个符号不同而已

这同时表明我们只要知道Fl,Fr的n/2个点值就可以求出F的n个点值。

于是复杂度就是T(x)=x+2*T(\dfrac{x}{2})=x\log x

即n\log n

接下来我们考虑插值，即知道n个点值会逆推得到系数。

这里设我们得到的点值为A

我们发现实际上我们要做的就是求解这样的一组方程组：

(\omega_{n}^0)^0*a_0+(\omega_{n}^0)^1*a_1+...+(\omega_{n}^0)^{n-1}*a_{n-1}=A(\omega_{n}^0) (\omega_{n}^1)^0*a_0+(\omega_{n}^1)^1*a_1+...+(\omega_{n}^1)^{n-1}*a_{n-1}=A(\omega_{n}^1) ...... (\omega_{n}^{n-1})^0*a_0+(\omega_{n}^{n-1})^1*a_1+...+(\omega_{n}^{n-1})^{n-1}*a_{n-1}=A(\omega_{n}^{n-1})

我们可以将之写作矩阵的形式。

于是我们就只要对左边的矩阵求逆了。

假设现在的式子是：V*X=F

那么我们需要求出V'使得V'*V*X=X=F*V'

这个矩阵可以构造出来，我们利用一下求和性质。

于是新矩阵的第i,j项为：

c_{i,j}=\sum_{k=0}^{n-1}a_{i,k}*b_{k,j}

如果我们令V'中的每一项都恰好为V的相反数（雾）

即(i,j)为(\omega_{n}^{-i})^j

于是可以得到：

c_{i,j}=\sum_{k=0}^{n-1}(\omega_{n}^{-i})^k*(\omega_{n}^k)^j

也就是：

c_{i,j}=\sum_{k=0}^{n-1}(\omega_{n}^{-i})^k*(\omega_{n}^j)^k c_{i,j}=\sum_{k=0}^{n-1}(\omega_{n}^{j-i})^k

由求和性质

当i=j时，其为n

否则其为0

所以我们构造出来的V'还要除一个n才是单位矩阵

也就有：V'*F/n=X

我们考虑求解V'*F

因为F固定，所以我们尝试将之看作一个新的系数

首先V'*F肯定是一个大小为1*n的矩阵

所以有：V'*F后得到的第i项其实就是：

A_0*(\omega_{n}^{-i})^0+A_1*(\omega_{n}^{-i})^1+...+A_{n-1}*(\omega_{n}^{-i})^{n-1}

你惊人的发现它变成了我们想要求出当x取\omega_{n}^{-0},\omega_{n}^{-1}...\omega_{n}^{-(n-1)}时的点值了

于是就再做一遍FFT，不过是每次取\omega_{n}^{-1}为基本增量。

复杂度O(n\log n)

下面补充NTT好了

我们尝试利用整数int类型来代替浮点数进行计算。

首先是阶的定义：（下面的等于(=)是同余）

如果有a^n=1(mod~p)

则称满足此条件最小的n为a模p的阶

若一个数g的阶为\phi(p)，则称g为p的原根

我们尝试使用原根代替单位根进行计算。

首先要有单位根的如下几条性质：

其次是折半性质$(DFT)$和求和性质$(IDFT)

我们尝试证明原根亦满足如下三条性质。

1.$若$g$为$p$的原根，则$g^i\%p$互不相等$(0\le i <\phi(p))

证明：

若存在两者相等，不妨设k>j且g^k=g^j(mod~p)

于是有：g^{k-j}=1(mod~p)

又k<\phi(p),j<\phi(p)

故k-j<\phi(p)

与定义矛盾，得证。

2.$我们定义原根$g_{n}^1$为$g^{\frac{p-1}{n}}

那么g_{n}^k=(g_n^1)^k，由性质1，可以得到原根互不相等。

不难得到：g_{n}^{2k}=g_{n/2}^k(mod~p)

把两边拆开，代入g_{n}^1即可。

其次：

我们接下来考虑将之指数看作一个角度为p-1的大圆，那么里面的n就是将之n等分，其余与单位根相似。

那么其乘法就是幅角不断累加。

于是可以得到：-g_n^k=p-g_n^k=g_n^{k+n/2}

这里就感性理解吧，我实在是不会证了，按照单位根的方法理解。

性质3:求和性质

求和性质利用了等比数列求和的定理，我们套用即可，因为是在\%p意义下，乘法操作仍然满足。

故g满足求和性质。

有了这三个性质，我们就能愉快的把原根代入单位根运算了。

然后因为FFT中的DFT和IDFT的n均为2^n

所以我们取的p要p-1整除足够大的n

比较常用的p是998244353，它的原根为3

且998244353=119*8388608(2^{23})+1

非递归版FFT和NTT

我们尝试优化它，从递归变成循环 仔细划分$FFT$会发现其呈现出树的结构 我们假设我们通过某种手段快速的得到了其最后的序列（叶子节点） 那么就可以两两往上合并来求答案 比如原序列： $0~ 1~ 2~ 3~ 4~ 5~ 6~ 7

最后得到的序列应该是:

0~4~2~6~2~5~3~7

我们打个表，发现：

000,001,010,011,100,101,110,111

变成了：

000,100,010,110,001,101,011,111

其实就是二进制反转了

考虑求出这个东东

貌似有一个递推式，不妨记R[i]为i反转后的值

于是有：R[i]=(R[i>>1]>>1)|((i\&1)<<(L-1))

简单讲讲 我们假设$[0,x)$的反转值都已经求出，现在想求出$x$的反转值 如果$x$为奇数，那么$x$可以写作$2*(x>>1)|1

那么有我们将x>>1反转后，可以将之右移一位，然后x还要在其最高位补上一个1

偶数类似，不用补1而已

这样就O(n)

然后因为这个关系双向，所以我们做swap即可。

#include<bits/stdc++.h>
using namespace std;
#define rep( i, s, t ) for( register int i = s; i <= t; ++ i )
#define re register
#define LL long long
int read() {
    char cc = getchar(); int cn = 0, flus = 1;
    while(cc < '0' || cc > '9') {  if( cc == '-' ) flus = -flus;  cc = getchar();  }
    while(cc >= '0' && cc <= '9')  cn = cn * 10 + cc - '0', cc = getchar();
    return cn * flus;
}
const int P = 998244353 ;
const int G = 3 ; 
const int Gi = 332748118 ;
const int N = 4e6 + 5 ; 
int n, m, R[N], L, limit ;
LL A[N], B[N] ;
LL fpow( LL x, int k ) {
    LL ans = 1, base = x ;
    while( k ) {
        if( k & 1 ) ans = ans * base % P ;
        k >>= 1, base *= base, base %= P ; 
    }
    return ans ; 
}
void NTT( LL *a, int type ) {
    for( int i = 0; i < limit; ++ i )
    if( i < R[i] ) swap( a[i], a[R[i]] ) ;
    for( int k = 1; k < limit; k <<= 1 ) {
        LL dg = fpow( ( type == 1 ) ? G : Gi, ( P - 1 ) / ( k << 1 ) ); //这里要注意
        for( int i = 0; i < limit; i += ( k << 1 ) ) //注意*2
        for( LL j = i, g = 1; j < i + k; ++ j, g = ( g * dg ) % P ) {
            LL Nx = a[j], Ny = ( a[j + k] * g ) % P;
            a[j] = ( Nx + Ny ) % P, a[j + k] = ( Nx - Ny + P ) % P ; 
        }
    }
}
void init() {
    while( limit <= n + m ) limit <<= 1, ++ L ; 
    for( int i = 0; i < limit; ++ i ) R[i] = ( R[i >> 1] >> 1 ) | ( ( i & 1 ) << ( L - 1 ) ) ;
}
signed main()
{
    n = read(), m = read(), limit = 1 ; 
    rep( i, 0, n ) A[i] = read() ; 
    rep( i, 0, m ) B[i] = read() ; 
    init(), NTT( A, 1 ), NTT( B, 1 ) ;
    rep( i, 0, limit ) A[i] = ( A[i] * B[i] ) % P ; 
    NTT( A, -1 ) ; int inv = fpow( limit, P - 2 ) ;
    rep( i, 0, n + m ) printf("%lld ", ( A[i] * inv ) % P ) ; 
    return 0;
}

【模板】多项式求逆

给定n-1次多项式F(x)，求解G(x)使得：

G(x)*F(x)=1(mod~x^n)

首先显然那些什么n次幂及以上的部分对答案没有影响

于是我们考虑递推

不妨设已知：H(x)*F(x)=1(mod~x^{n/2})

则显然：G(x)*F(x)=1(mod~x^{n/2})

故做差，得到：

G(x)-H(x)=0(mod~x^{n/2})

平方得到：

(G(x)-H(x))^2=0(mod~x^n)

于是有：

G^2+H^2-2G*H=0(mod~x^n)

两边同时乘F得到：

G+F*H^2-2H=0(mod ~x^n)

于是有：

G=2H-F*H^2

如此递推即可。

从1递推到2^k>n比较合适

我也不知道为什么复杂度是O(n\log n)的，我感觉是O(n \log^2n)的...

这道题真的是令人调得丧心病狂，心态爆炸

按着题解得来打还是挂了我也不知道为什么。

最后有一种常数比较小的做法，就是不做乘法，每次算点值乘起来就好了

这样只要跑3次NTT

以及，以后写代码先自己想，别照着题解打，今天不是照着题解打早调出来了...

P4721 【模板】分治 FFT

好像有分治的做法

但我们搞求逆算了

考虑构造一个函数F(x)=\sum_{i=0}^{\infty}f_i*x^{i}

当然这样的f我们无限的生成下去

类似的构造一个G(x)，没输入的地方都取0

于是显然有：F(x)*G(x)+f_0=F(x)

然后就有：F(x)=f_0*(1-G(x))^{-1}

因为我们只要求出n-1项，所以有：

F(x)=f_0*(1-G(x))^{-1}(mod~x^n)

因为f_0=1，所以我们直接跑多项式求逆就好了...

【模板】任意模数NTT

貌似有很多做法，先记着一个我会的，最慢的，7FFT

将每一项的系数拆成k1*m+b1=a,k2*m+b2=b

于是显然有：m*K1(x)+B1(x)=A(x)，另一个同理。

我们就分别对K1,B1,K2,B2做一遍DFT

然后考虑合并

其在x处的取值，我们做乘法本来是A(x)*B(x)

变成(m*K1+B1)*(m*K2+B2)

然后暴力展开，发现还要对K1*K2,K1*B2+K2*B1,B1*B2做一遍IDFT

最后合并系数即可，7FFT

以后学了别的会慢慢补充.......

MTT

考虑对于两个多项式A,B

构造P=A+iB,Q=A-iB

显然有，记P_k表示P在\omega_{n}^k处的取值

记录X=\dfrac{2\pi*k}{n}

P_k=A(\omega_{n}^k)+i*B(\omega_n^k) P_k=\sum_{j=0}^kA_j(\omega_n^{jk})+i*B_j(\omega_n^{jk}) P_k=\sum_{j=0}^k(A_j+i*B_j)*\omega_{n}^{jk} P_k=\sum_{j=0}^k(A_j+i*B_j)*(cosX+isinX) Q_k=\sum_{j=0}^kA_j(\omega_n^{jk})-i*B_j(\omega_n^{jk}) Q_k=\sum_{j=0}^k(A_j-i*B_j)*\omega_{n}^{jk} Q_k=\sum_{j=0}^k(A_j-i*B_j)(cosX+isinX) Q_k=\sum_{j=0}^k(A_jcosX+B_jsinX)+i*(A_jsinX-B_jcosX) Q_k=conj(\sum_{j=0}^k(A_jcosX+B_jsinX)-i*(A_jsinX-B_jcosX)) Q_k=conj(\sum_{j=0}^kA_j(cosX-i*sinX)+i*B_j(cosX-i*sinX)) Q_k=conj(\sum_{j=0}^k(A_j+i*B_j)(cosX-i*sinX)) Q_k=conj(\sum_{j=0}^k(A_j+i*B_j)\omega_{n}^{-jk}) Q_k=conj(\sum_{j=0}^k(A_j+i*B_j)(\omega_{n}^{n-k})^j) Q_k=conj(P_{n-k})

所以我们只要DFT出P，就能还原Q了

然后求出了P和Q后

就有：

A=\dfrac{P+Q}{2},B=\dfrac{P-Q}{2}*(-i)

然后对A*B做一遍IDFT即可

3$次$FFT->2FFT

上述优化过程后得到的FFT叫做MTT

接下来考虑优化IDFT

假设我们有A的点值和B的点值

根据求和性质，我们只需要求出其在\omega_{n}^{-k}的点值就可以get到系数

不妨构造P=\dfrac{A+iB}{2},Q=\dfrac{A-iB}{2}

P_k=\sum_{j=0}^k(A_j+i*B_j)*(\omega_{n}^{n-k})^j Q_k=\sum_{j=0}^kA_j(\omega_n^{n-k})^j-i*B_j(\omega_n^{n-k})^j

令X=\dfrac{2\pi*(-kj)}{n}

同上可得：

Q_k=conj(P_{n-k})

貌似还有MTT2，可以做到1.5FFT搞出结果，不会告辞

【模板】多项式除法

给定F,G

求出Q,R使得：

F=G*Q+R

注意到如下性质：

f(\dfrac{1}{x})*x^n=f_r(x) 于是有： $$F(x)=G*Q(x)+R(x)$$ $$F(\dfrac{1}{x})*x^n=G*Q(\dfrac{1}{x})*x^n+R(\dfrac{1}{x})*x^n$$ $$F(\dfrac{1}{x})*x^n=G(\dfrac{1}{x})*x^m*Q(\dfrac{1}{x})*x^{n-m}+R(\frac{1}{x})*x^{m-1}*x^{n-m+1}$$ $$F_r(x)=G_r(x)*Q_r(x)+R_r(x)*x^{n-m+1}$$ 同时取模得到： $$F_r(x)=G_r(x)*Q_r(x)(~mod~ x^{n-m+1}~)$$ $$F_r(x)*G_r^{-1}(x)=Q_r(x)(~mod~ x^{n-m+1}~)$$ 于是就求个逆，然后做乘法$......

注意用NTT

$F^{-1}(mod~x^n)$与$F^{-1}(mod~x^{n+1})$是不同的啊QAQ -------------------- 例题： 力 仔细观察$Ei$的定义，可以得出： $$E1=0*q1+\dfrac{q2}{-1}+\dfrac{q3}{-2^2}+...+\dfrac{q_n}{-(n-1)^2}$$ $$E2=\dfrac{q1}{1}+0*q2+...+\dfrac{q_n}{-(n-2)^2}$$ 你惊人的发现它仅仅是系数平移了一位而已。 就是一个加法卷积的形式。 之所以平移一位可以将$E$看作卷积乘法得到的多项式，那么有$E1,E2...En$往后依次$+1

即

E1=\sum_{i+j=k}q_i*b_j E2=\sum_{i+j=k+1}q_i*b_j

于是构造一组b和原来q卷起来就好了，第i+n位为Ei

或许一直以来对FFT的理解都不够吧...

考虑一种变换，应该满足：

FFT(C,i)=FFT(A,i)*FFT(B,i)

我们好像实际上C_i=\sum_{j+k=i}A_j*B_k

接下来就不会了qwq

仔细看了myy的博客貌似搞懂了，补充一下：

不过貌似这篇文章太长了现在很卡

简单来讲从加法卷积入手，将x换成单位根，然后单位根反演，就会发现只需要DFT后再乘起来再IDFT

qwq