多项式全家桶

多项式全家桶

多项式乘法

多项式乘法$CX=AX\times BX$中的系数满足加法卷积$c_i=\Sigma a_j b_k [i==j+k]$，可以将一个乘法看作卷积，反之也是。

一般来说，多项式直接相乘的时间复杂度都是$O(n^2)$级别的，使用压位等技巧也只能做到常数上的优化，面对更大的数据规模就捉襟见肘了。

既然是卷积，就可以考虑卷积通用的方法进行优化。如果能找到映射$f$,使得在$f$下点乘的性质成立，并且可求逆运算，即

$$A\rightarrow f(A)\\ B\rightarrow f(B)\\ f(C)=f(A)*f(B) \Leftrightarrow f(c_i)=f(a_i)*f(b_i)\\ f(C)\rightarrow C$$

就有可能找到更好的做法。

快速傅里叶变换

傅里叶变换是一个物理上的方法，虽然他们名字相同，但是不明白他和这个过程有什么使用上的联系。得知复数单位根有着优良的性质，自然想到要用他来优化多项式乘法

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

对于多项式

$$f(x)=\Sigma_0^n c_i x_i \\ (n=2^k)$$

可以把他进行奇偶分解，记奇偶子式分别为$h(x),g(x)$

$$f(x)=\Sigma_0^n c_i x_i\\ =\Sigma c_j x_j + x \Sigma c_k x_k\\ (j\&2 = 0 ,k\&2 = 1)\\ =g(x)+xh(x)$$

代入$x_i=\omega_n^i$，

$$f(\omega)=g(\omega^2)+\omega h(\omega^2)$$

这个求$f$的过程是典型的折半递归，得到了一个时间复杂度为$O(nlogn)$的求解法。求出了映射$A\rightarrow f(A)$，可以在线性的时间内实现点乘，接下来考虑逆变换$f(C)\rightarrow C$。

如果把dft过程看作一个矩阵乘法，那么idft就是矩阵的逆。其实用一个更好理解的想法，再次代入$\omega^{-1}$就逆映射。需要注意，直接代入求得的并不是原序列，由线代的知识可知，$A^{-1}=\frac{A^*}{|A|}$，需要再除以序列的长度n。

代码:

//op=1 or -1 
void fft(int n,cpl *a,int op){
    if(n==1)return ;
    //奇偶分类
    for(int i=0;i<n;i+=2)temp[i/2]=a[i],
    temp[i/2+n/2]=a[i+1];
    ffor(i,0,n-1)a[i]=temp[i];
    //分治 
    fft(n/2,a,op),fft(n/2,a+n/2,op);
    //合并 
    cpl Wn={cos(2.0*PI/n),sin(2.0*PI/n)*op},w=1; 
    for(int i=0,nl=n/2;i<nl;++i,w=w*Wn){
        cpl t1=a[i],t2=w*a[i+nl];
        a[i]=t1+t2;
        a[i+nl]=t1-t2;
    }
}
void mult(){
    fft(n,a1,1),fft(n,a2,1);
    ffor(i,0,n-1)a3[i]=a1[i]*a2[i];
    //求逆，补精度
    fft(n,a3,-1);
    ffor(i,0,nout)a3[i]=int(a3[i].real()/n+0.5);
}

迭代fft

考虑在$f(a_i)$中所有有贡献的元素，可以发现他们的下标都等于$i>>k$。在$f(\omega)=g(\omega^2)+\omega h(\omega^2)$划分递归的过程，x也在作二进制下的右移。

如果能提前得到在$f(a_i)$所涉及的位置，就不需要递归，可以剩下栈开销与优化常数。

自底而顶地合并元素，可以发现当$i=rev[i]$时，满足fft的性质，其中

$$rev[i]=rev[i/2]/2+(i\&1)*n/2;$$

代码:

//op=1 or -1 
void fft(cpl *a,int op){
    //奇偶分类
    ffor(i,0,n-1)if(i<rev[i])swap(a[i],a[rev[i]]);
    for(int mid=1;mid<n;mid<<=1){
        cpl Wn={cos(PI/mid),sin(PI/mid)*op};
        for(int p=0;p<n;p+=(mid<<1)){
            cpl w=1;
            for(int k=0;k<mid;++k,w=w*Wn){
                cpl t1=a[p+k],t2=w*a[mid+p+k];
                a[p+k]=t1+t2;
                a[mid+p+k]=t1-t2;
            }
        }
    }
}

快速数论变换

fft使用了复数，有运行常数和精度的劣势，而且许多题目要对998244353取模，一个基于整数的变换算法会更合适。

由数论知识可知，对某些质数$P$，存在原根$g$，使得

$$g_n^k=-g_n^{k+n/2}(mod P)\\ g_n^{2k}=g_{n/2}^{k}(mod P)$$

那么就可以利用类似于上述的过程来得到变换

代码:

void ntt(ll *a,int op){
    //奇偶分类
    ffor(i,0,n-1)if(i<rev[i])swap(a[i],a[rev[i]]);
    for(int mid=1;mid<n;mid<<=1){
        ll Wn=fsp((op==1?G:Gi),(P-1)/(mid<<1));
        for(int p=0;p<n;p+=(mid<<1)){
            ll w=1;
            for(int k=0;k<mid;++k,w=w*Wn%P){
                ll t1=a[p+k],t2=w*a[mid+p+k]%P;
                a[p+k]=(t1+t2)%P;
                a[mid+p+k]=(t1-t2+P)%P;
            }
        }
    }
}
void mult(){
    ntt(a1,1),ntt(a2,1);
    ffor(i,0,n-1)a3[i]=a1[i]*a2[i]%P;
    ntt(a3,-1);
    ffor(i,0,nout)a3[i]=a3[i]*fsp(n)%P;
}

多项式运算

求导，积分

对于一个多项式函数，求导和积分的原理都非常简单。在计算机中的实现就是移动项并变换相应系数。
在ntt的做法中，可以预处理阶乘和逆阶乘。

代码:

ll fac[N],ifac[N];
//取k阶导或求k次积分
int deviate(ll *a,int n,int k=1){
    ffor(i,0,n-k-1)a[i]=a[i+k]*fac[i+k]%P*ifac[i]%P;
    return n-k;
}
int integrate(ll *a,int n,int k=1){
    rfor(i,n-1,0)a[i+k]=a[i]*ifac[i+k]%P*fac[i]%P;
    ffor(i,0,k-1)a[i]=0;
    return n+k;
}
void init(){
    fac[0]=ifac[0]=1;
    ffor(i,1,N-1)fac[i]=fac[i-1]*i%P;
    ifac[N-1]=fsp(fac[N-1]);
    rfor(i,N-1,1)ifac[i-1]=ifac[i]*i%P;
} 

求逆

$F(x)$的多项式逆元定义为$F(x)G(x)\equiv 1 (mod x^n)$。这个逆元只针对原多项式长度范围内的部分，对于乘积超出的部分直接忽略。

考虑递归求解。在递归边界$n=1$，$G_0=f_0^{-1}$。

在模$x^{\frac{n}{2}}$的意义下，定义$F(X)$的逆多项式$H(x)$。$H(x)$通过递归逐次求解。

$$F(x)H(x)\equiv 1(mod x^{\frac{n}{2}})\\ F(x)G(x)\equiv 1(mod x^n)\\ \downarrow \\ F(x)(H(x)-G(x))\equiv 0(mod x^{\frac{n}{2}})\\ H(x)-G(x)\equiv 0(mod x^{\frac{n}{2}})$$

可以看出，在低位$x^{\frac{n}{2}}$以下，$H(X)=G(X)$，即$H(X)$就是$G(X)$的低位部分。至于高位部分，考虑平方

$$G^2(x)\equiv 2H(X)G(x)-H^2(x)(mod x^n)\\ F(x)*G^2(x)\equiv F(x)*2H(x)G(x)-F(x)*H^2(x)(mod x^n)\\ \downarrow \\ G(x)\equiv 2H(x)-F(x)*H^2(x)(mod x^n)$$

这样，通过$\frac{n}{2}$项式$H(X)$，就可以得到$n$项式$G(X)$。

根据fft的经验，二进制分治的递归可以改写成迭代，求逆也可以改用迭代，这样时空开销更小。

代码:

// n:多项式f原长度
ll inv_f[N];
void poly_inv(ll *f,ll *g,int n){
    ffor(i,0,n-1)g[i]=0;g[0]=fsp(f[0]);
    for(ll len=1;len<(n<<1);len<<=1){
        ll lim=len<<1;
        ffor(i,0,lim-1)inv_f[i]=0;
        ffor(i,0,len-1)inv_f[i]=f[i];
        ntt(inv_f,lim),ntt(g,lim);
        ffor(i,0,lim-1)g[i]=(P+2ll-inv_f[i]*g[i]%P)%P*g[i]%P;
        ntt(g,lim,-1);
        ffor(i,len,lim-1)g[i]=0;
    }
}

开根

$F(x)$的根式定义为$G^2(x)\equiv F(x) (mod x^n)$。按照求逆的套路，考虑分治合并。

当$H(x)$是$F(x)$在一半长度下的根式

$$H(x)\equiv G(x)(mod x^{\frac{n}{2}})\\ H(x)-G(x)\equiv 0(mod x^{\frac{n}{2}})\\ (H(x)-G(x))^2\equiv 0(mod x^n)\\ \downarrow \\ G^2(x)+H^2(x)-2G(x)H(x)\equiv 0(mod x^n)\\ F(x)+H^2(x)-2G(x)H(x)\equiv 0(mod x^n)\\ \frac{F(x)+H^2(x)}{2H(x)}\equiv G(x)(mod x^n)\\$$

递归，求$H(x)$的逆元，就可由半长根式得全长根式

代码:

ll sqrt_f[N],sqrt_inv[N];
void poly_sqrt(ll *f,ll *g,int n){
    ffor(i,0,n-1)g[i]=0;g[0]=1;
    for(ll len=1;len<(n<<1);len<<=1){
        ll lim=len<<1;
        ffor(i,0,lim-1)sqrt_f[i]=sqrt_inv[i]=0;
        ffor(i,0,len-1)sqrt_f[i]=f[i];
        poly_inv(g,sqrt_inv,len);
        ntt(sqrt_f,lim),ntt(sqrt_inv,lim);
        ffor(i,0,lim-1)sqrt_f[i]=sqrt_inv[i]*sqrt_f[i]%P;
        ntt(sqrt_f,lim,-1);
        ffor(i,0,len-1)g[i]=(sqrt_f[i]+g[i])%P*inv2%P;
        ffor(i,len,lim-1)g[i]=0;
    }
}