多项式全家桶

Soulist · 2019-10-08 19:59:02 · 个人记录

实现卷积的算法

First.FFT

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

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

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

这是一个加法卷积。

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

接下来是FFT的核心思想：

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

将点值两两乘起来。

将点值换成系数表达。

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

前方高能：

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

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

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

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

显然\omega_n^{0}=1

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

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

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

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

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

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

第二：

原因很简单

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

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

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

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

于是有：

稍微化简一下得：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}

不妨设：

于是有：

接下来我们代入\omega_n^k

得到：

于是：

这里假设k< n/2

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

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

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

于是有：

得到：

即：

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

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

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

即n\log n

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

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

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

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

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

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

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

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

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

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

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

于是可以得到：

也就是：

由求和性质

当i=j时，其为n

否则其为0

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

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

我们考虑求解V'*F

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

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

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

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

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

复杂度O(n\log n)

Second.NTT

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

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

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

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

若一个数g的阶为\phi(p)，则称g为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)

与定义矛盾，得证。

那么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

Third.非递归版FFT和NTT

最后得到的序列应该是:

我们打个表，发现：

变成了：

其实就是二进制反转了

考虑求出这个东东

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

于是有：R[i]=(R[i>>1]>>1)|((i\&1)<<(L-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;
}

Fourth.MTT/任意模数FFT

貌似有很多做法，先记着一个我会的，最慢的，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 jk}{n}

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

然后求出了P和Q后

就有：

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

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

接下来考虑优化\rm IDFT

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

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

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

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

同上可得：

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

这样大概可以将两次算点值做到一次同时get

然后就可以将7FFT这种丢人的东西变成4FFT...

我寻思大家为什么都用鱼写的那个奇怪的东西叫做三次变两次啊...真正的三次变两次不应该是这玩意儿吗...这是真正意义上同时DFT出两个点值啊...

然后就要开始上全家桶了qwq

华丽的分割线

多项式全家桶

First.多项式求逆

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

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

于是我们考虑递推

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

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

故做差，得到：

平方得到：

于是有：

两边同时乘F得到：

于是有：

如此递推即可。

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

复杂度O(n\log n)

关于证明：

于是T'(x)=\log x*(O'(x))>T(x)

所以T'(x)=\log x *O'(x)=n\log n

Second. 【模板】分治 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_0=1，所以我们直接跑多项式求逆就好了...

Third.【模板】多项式除法

给定F,G

求出Q,R使得：

注意到如下性质：

Sixth.多项式\ln/多项式\exp(4\&5为泰勒展开+牛顿迭代)

这个东西要前置芝士...

1.求导

定义式：

关于求导运算的一些结论：

证明：

于是就有：

因为\lim_{dx\to 0}f(x+dx)=f(x)

所以就有：

原式=

于是：

证明：

可以构造dg=\dfrac{}{dx}

咕咕咕...

注意到 f^k 是一个复合导，所以 (f^k)'=f'\times k\times f^{k-1}

所以：

一些比较常用的结论：

3.式比较有用

2.积分

积分和求导为逆运算，定义式：

然后积分可加...(面积可加)

于是同理求导可加

然后积分可以直接乘以倍数/系数，反过来求导也可以qwq

Fourth.泰勒展开

首先是一个函数的泰勒展开：

其中f^i表示i阶导数

然而在精度允许的条件下由于阶乘非常的大所以只需要把n开到12,13左右即可

然后由于我们可以随便选择一个x_0进行泰勒展开，所以一般选择0

于是我们就可以将原式转化为

Firth.牛顿迭代

用于求解0点

与二分不同的是牛顿迭代的求法是对一个函数求导，降次，求解降次后的0点以逼近原函数的0点

多项式牛顿迭代即求解满足G(F(x))=0(mod~~x^n)的解F(x)

假设现在已经求出了：

考虑拓展到x^n次意义下，可以在G(F_0(x))处将G(F(x))泰勒展开变成：

因为我们有：G(F(x))=0(mod~~x^n)

所以就有：

又G(F_0(x))=0(mod~~x^{n/2})

然后又由于理论上应该有：F(x)和F_0(x)的n/2项应该相等

所以那个奇怪的减法会导致这个函数没有最后n/2项

于是就发生了神奇的结论！

没有最后n/2项

所以(F(x)-F_0(x))^2没有最后n-1项

于是令人惊喜的是我们泰勒展开的第二项就是(F(x)-F_0(x))^2

所以我们居然只需要展开两项即可！！！woc

所以就得到了神奇的结论：

于是就有：

于是就可以倍增处理了...

问题的关键主要是对G'(F_0(x))求导才对吧...

下面给出两个对于G'(F_0(x))求导的例子：

Sixth.多项式\ln/多项式\exp

1.多项式\ln

然后接下来推一下多项式\rm ln

两边同时求导得：

令u=A(x)

又\ln'(u)=\dfrac{1}{u}

所以得：

所以就只需要计算出A'(x)即可

因为求导与积分互逆，所以求导可加=积分可加，个人感觉

首先有对一个数求导是0

所以我们知道

然后在求一下A(x)的逆A^{-1}(x)

然后将A'(x)和A^{-1}(x)乘起来就好了qwq

然后还原回去就套用上面求导那里的公式

2.多项式\exp

现在想要求解：

于是两边同时取\ln得：

可以设G(B(x))=\ln B(x)-A(x)于是问题变成求解它的0点

这个多项式里面A(x)为定值，而B(x)为变量

我们假设有：

问题是求解B(x)

于是代入牛顿迭代的公式得到：

因为G'(H(x))=\ln H(x)-A(x)

而其中A(x)为常量，由于求导的可加性，A'(x)=0，简单求和可以得到：

所以我们可以得到：

然而因为求解\ln H(x)的复杂度是O(n\log n)的

于是求解\exp的复杂度也是T(x)=x\log x+T(x/2)

即O(n\log n)

不过常数比较大就是了......

Seventh.多项式开根

与上面类似

假设有

两边同时平方得到：

可以设一个函数G(B(x))=B^2(x)-A(x)

代入牛顿迭代的式子就是：

然后类似的，将A(x)看作常量，那么可以得到：G'(B(x))=2B(x)

所以得到：

于是要求逆，复杂度O(n\log n)

#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 gi() {
    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 Gi = 332748118 ;
const int G = 3 ; 
const int N = 5e5 + 5 ; 
int n, R[N], L, limit ;
LL A[N] ; 
LL fpow( LL x, LL k ) {
    LL ans = 1, base = x ;
    while( k ) {
        if( k & 1 ) ans = ( ans * base ) % P ;
        base = ( 1ll * base * base ) % P, k >>= 1 ; 
    }
    return ans % P ; 
} 
LL Gx[N], H[N], Inv ;
void Init( int x ) {
    limit = 1, L = 0 ;
    while( limit < x ) limit <<= 1, ++ L ; //因为都是2的整数次幂，实际上并不需要 <=
    //这样会造成8倍长度的数组，最后TLE 
    rep( i, 0, limit - 1 ) R[i] = ( R[i >> 1] >> 1 ) | ( ( i & 1 ) << ( L - 1 ) ) ;
    Inv = fpow( limit, P - 2 ) ;
}
void NTT( LL *a, int type ) {
    rep( i, 0, limit - 1 ) if( R[i] > i ) swap( a[i], a[R[i]] ) ;
    for( re int k = 1; k < limit; k <<= 1 ) {
        int d = fpow( ( type == 1 ) ? G : Gi, ( P - 1 ) / ( k << 1 ) ) ;
        for( re int i = 0; i < limit; i += ( k << 1 ) )
            for( re LL j = i, g = 1; j < i + k; ++ j, g = ( d * g ) % P ) {
                LL Nx = a[j], Ny = ( a[j + k] * g ) % P ;
                a[j] = ( Nx + Ny ) % P, a[j + k] = ( Nx - Ny + P ) % P ; 
            }
    }
    if( type != 1 ) rep( i, 0, limit - 1 ) a[i] = ( a[i] * Inv ) % P ;
}
void GetInv( LL *a, int x ) {
    int k = 1 ; memset( H, 0, sizeof(H) ), memset( Gx, 0, sizeof(Gx) ), H[0] = fpow( A[0], P - 2 ) ;
    while( k < x ) { //这里同理 
        k <<= 1, Init( k * 2 ) ;
        rep( i, 0, k - 1 ) Gx[i] = a[i] ;
        NTT( Gx, 1 ), NTT( H, 1 ) ;
        rep( i, 0, limit ) H[i] = ( 2ll - Gx[i] * H[i] % P + P ) * H[i] % P ; 
        NTT( H, -1 ) ; rep( i, k, limit ) H[i] = 0 ;
    } 
}
LL Gt[N], B[N] ; 
void Sqrt( LL *a, int x ) {
    int k = 1 ; LL Iv = 499122177 ; B[0] = 1 ; 
    while( k <= x ) {
        k <<= 1 ; rep( i, 0, k - 1 ) Gt[i] = a[i] ; 
        GetInv( B, k ), Init( k * 2 ) ;
        NTT( Gt, 1 ), NTT( H, 1 ), NTT( B, 1 ) ;
        rep( i, 0, limit ) B[i] = 1ll * Iv * ( B[i] % P + Gt[i] * H[i] % P ) % P ;
        NTT( B, -1 ) ; rep( i, k, limit ) B[i] = 0 ; 
    }
}
signed main()
{
    n = gi() ;
    rep( i, 0, n - 1 ) A[i] = gi() ; 
    Sqrt( A, n ) ;
    rep( i, 0, n - 1 ) printf("%d ", B[i] ) ;
    return 0 ;
}