多项式学习笔记

lizhous · 2023-02-10 13:54:09 · 个人记录

乘法

概念

多项式乘法时，我们发现暴力乘十分缓慢，但是点值乘十分快速。考虑求 A 和 B 的卷积。

一个 n 次多项式可以被 n+1 个点确定。

设多项式 A(x) 的系数为 (a_0,a_1,\cdots,a_n)

对其奇偶分类得 A(x)=\sum\limits a_{2i}*x^{2i}+\sum a_{2i+1}*x^{2i+1}

提取 n 得 A(x)=\sum\limits a_{2i}*x^{2i}+x*\sum a_{2i+1}*x^{2i}

化简得 A(x)=\sum\limits a_{2i}*(x^2)^i+x*\sum a_{2i+1}*(x^2)^i

设 (\frac{n}{2}-1) 次多项式 A_0=\sum\limits a_{2i}*x^{2i},A_1=\sum\limits a_{2i+1}*x^{2i}

则 A(x)=A_0(x^2)+x*A_1(x^2)

我们发现可以递归

考虑继续优化。

复数

定义 i^2=-1

则所有形如 z=a+b*i 的数 z 组成的集合为复数集，记为 C。而 a 称为 z 的实部， b 为虚部。

复平面

类似坐标系，纵虚横实。复数 z 对应从原点指向 (a,b) 的向量。幅角为实轴横方向与向量的夹角，记为 θ。

加法：与向量加法一样

乘法：复数相乘，模长相乘，幅角相加。(a+bi)*(c+di)=(ac-bd)+(bc+ad)i

单位根

在复平面上，以原点为圆心，1 为半径作圆，所得的圆叫单位圆。以圆点为起点，圆的 n 等分点为终点，做 n 个向量，设幅角为正且最小的向量对应的复数为 \omega_n，称为 n 次单位根。剩余 n-1 个点则为 \omega_n^2,\omega_n^3,\omega_n^4,\cdots。

单位根的幅角为周角的 \dfrac{1}{n}。

欧拉公式有：w_n^k=\cos k\frac{2π}{n}+i\sin k\frac{2π}{n}

性质 1: \omega_{n}^{k}=\omega_{2n}^{2k}

性质 2: \omega _{n}^{k+\frac{n}{2}}=-\omega _{n}^{k}

性质 3: w_n^n=1

图上自证不难

FFT

带入 x=w_n^k(k<\frac{1}{2}n) 得

带入 x=w_n^{k+\frac{n}{2}}(k<\frac{1}{2}n) 得

不难发现，两式的值可以互相转换求得。那么当 k=[0,\frac{n}{2}-1] 时，k+\frac{n}{2}=[\frac{n}{2},n] 两式分别处理了一半区间。递归处理，复杂度 O(n\log n)。

IFFT

函数转点并相乘后，需要转化回系数。此时使用傅里叶逆变换。

设向量 c_1,c_2,\cdots c_{n-1} 满足 c_k=\sum\limits^{n}_{j=1} y_jx^j。当 x=w_n^{-k} 时

所以和之前一样递归，带入 x=w_n^{-k} 即可。其实就是最后反过来。

代码

递归实现即可

使用STL：complex实现复数

#include<iostream>
#include<cstdio>
#include<complex>
#include<cmath>
#include<algorithm>
using namespace std;
double pi=acos(-1);
void FFT(complex<double> *A,int N)
{
    if(N==1) return;
    complex<double> A1[N>>1],A2[N>>1];
    for(int i=0;i<N;i+=2)
    {
        A1[(i)>>1]=A[i];
        A2[(i)>>1]=A[i+1];
    }
    int N2=N>>1;
    FFT(A1,N2);
    FFT(A2,N2);
    complex<double> W=complex<double>(cos(pi/N2),sin(pi/N2));
    complex<double> now=complex<double>(1,0);
    for(int i=0;i<N2;i++)
    {
        A[i]=A1[i]+now*A2[i];
        A[i+N2]=A1[i]-now*A2[i];
        now=now*W;
    }
}
void IFFT(complex<double> *A,int N)
{
    FFT(A,N);
    reverse(A+1,A+N);
}

但是，递归常数很大。

改成迭代就可以了。

迭代优化

观察 FFT 组成的递归树。

0 1 2 3 4 5 6 7
0 2 4 6|1 3 5 7
0 4|2 6|1 5|3 7
0|4|2|6|1|5|3|7

实质上是二进制从低到高位逐位排序，所以把每个数二进制翻转后变成正常排序，就是

0 1 2 3 4 5 6 7

那么直接翻转即可。

之后就是模拟递归的流程，写成循环。

#include<iostream>
#include<cstdio>
#include<complex>
#include<cmath>
#include<algorithm>
using namespace std;
double pi=acos(-1);
int tax[5000001];
void Rev(int n,int cn)
{
    for(int i=0;i<n;i++)
    {
        int rw=0,ii=i;
        for(int z=1;z<=cn;z++)
        {
            rw<<=1;
            rw^=ii&1;
            ii>>=1;
        }
        tax[i]=rw;
    }
}
void FFT(complex<double> *A,int N)
{
    complex<double> B[N+2];
    for(int i=0;i<N;i++)
    {
        B[i]=A[i];
    }
    for(int i=1;i<N;i++)
    {
        A[i]=B[tax[i]];
    }
    for(int len=2,M=1;len<=N;M=len,len<<=1)
    {
        complex<double> W(cos(pi/M),sin(pi/M));
        for(int l=0,r=len-1;r<=N;l+=len,r+=len)
        {
            complex<double> w(1.0,0.0);
            for(int p=l;p<l+M;p++)
            {
                complex<double> x=A[p]+w*A[p+M];
                complex<double> y=A[p]-w*A[p+M];
                A[p]=x;
                A[p+M]=y;
                w*=W;
            }
        }
    }
}
void IFFT(complex<double> *A,int N)
{
    FFT(A,N);
    reverse(A+1,A+N);
}
int n,m;
complex<double> F[5000001],G[5000001];
int main()
{
    scanf("%d%d",&n,&m);
    for(int i=0,_;i<=n;i++)
    {
        scanf("%d",&_);
        F[i]=_;
    }
    for(int i=0,_;i<=m;i++)
    {
        scanf("%d",&_);
        G[i]=_;
    }
    int p2=1,cn=0;
    while(p2<=m+n) p2<<=1,cn++;
    Rev(p2,cn);
    FFT(F,p2);
    FFT(G,p2);
    for(int i=0;i<p2;i++)
    {
        F[i]*=G[i];
    }
    IFFT(F,p2);
    for(int i=0;i<=n+m;i++)
    {
        printf("%d ",int(F[i].real()/p2+0.5));
    }
}

NTT

模意义下，把 FFT 中的单位根换成原根即可。

其中 w_{n}^{1}≡g^{\frac{mod-1}{n}}\pmod {mod}

NTT板子

#include<iostream>
#include<cstdio>
#include<complex>
#include<cmath>
#include<algorithm>
#define mod 998244353
#define int long long
using namespace std;
int tax[5000001];
const int Gen=3;//原根 
int n,m;
int F[5000001],G[5000001];
int poww(int a,int b)
{
    int ret=1;
    while(b)
    {
        if(b&1)
        {
            ret*=a;
            ret%=mod;
        }
        a*=a;
        a%=mod;
        b>>=1;
    }
    return ret;
}
void Rev(int n,int cn)
{
    for(int i=0;i<n;i++)
    {
        int rw=0,ii=i;
        for(int z=1;z<=cn;z++)
        {
            rw<<=1;
            rw^=ii&1;
            ii>>=1;
        }
        tax[i]=rw;
    }
}
void FFT(int *A,int N)
{
    int B[N+2];
    for(int i=0;i<N;i++)
    {
        B[i]=A[i];
    }
    for(int i=1;i<N;i++)
    {
        A[i]=B[tax[i]];
    }
    for(int len=2,M=1;len<=N;M=len,len<<=1)
    {
        int W=poww(Gen,(mod-1)/len);
        for(int l=0,r=len-1;r<=N;l+=len,r+=len)
        {
            int w=1;
            for(int p=l;p<l+M;p++)
            {
                int x=(A[p]+w*A[p+M]%mod)%mod;
                int y=(A[p]-w*A[p+M]%mod)%mod;
                A[p]=x;
                A[p+M]=y;
                w=w*W%mod;
            }
        }
    }
}

void IFFT(int *A,int N)
{
    FFT(A,N);
    reverse(A+1,A+N);
    int invn=poww(N,mod-2);
    for(int i=0;i<=N;i++)
    {
        A[i]=((A[i]*invn)%mod+mod)%mod;
    }
}
signed main()
{
    scanf("%lld%lld",&n,&m);
    for(int i=0;i<=n;i++)
    {
        scanf("%lld",&F[i]);
    }
    for(int i=0;i<=m;i++)
    {
        scanf("%lld",&G[i]);
    }
    int p2=1,cn=0;
    while(p2<=m+n) p2<<=1,cn++;
    Rev(p2,cn);
    FFT(F,p2);
    FFT(G,p2);
    for(int i=0;i<p2;i++)
    {
        F[i]=F[i]*G[i]%mod;
    }
    IFFT(F,p2);
    int invn=poww(p2,mod-2);
    for(int i=0;i<=n+m;i++)
    {
        printf("%lld ",F[i]);
    }
}

封装版本

#include<iostream>
#include<cstdio>
#include<complex>
#include<cmath>
#include<algorithm>
#define mod 998244353
#define int long long
using namespace std;
int tax[5000001];
const int Gen=3;//原根 
int n,m;
int F[5000001],G[5000001];
int poww(int a,int b)
{
    int ret=1;
    while(b)
    {
        if(b&1)
        {
            ret*=a;
            ret%=mod;
        }
        a*=a;
        a%=mod;
        b>>=1;
    }
    return ret;
}
void Rev(int n,int cn)
{
    for(int i=0;i<n;i++)
    {
        int rw=0,ii=i;
        for(int z=1;z<=cn;z++)
        {
            rw<<=1;
            rw^=ii&1;
            ii>>=1;
        }
        tax[i]=rw;
    }
}
void FFT(int *A,int N)
{
    int B[N+2];
    for(int i=0;i<N;i++)
    {
        B[i]=A[i];
    }
    for(int i=1;i<N;i++)
    {
        A[i]=B[tax[i]];
    }
    for(int len=2,M=1;len<=N;M=len,len<<=1)
    {
        int W=poww(Gen,(mod-1)/len);
        for(int l=0,r=len-1;r<=N;l+=len,r+=len)
        {
            int w=1;
            for(int p=l;p<l+M;p++)
            {
                int x=(A[p]+w*A[p+M]%mod)%mod;
                int y=(A[p]-w*A[p+M]%mod)%mod;
                A[p]=x;
                A[p+M]=y;
                w=w*W%mod;
            }
        }
    }
}
void IFFT(int *A,int N)
{
    FFT(A,N);
    reverse(A+1,A+N);
    int invn=poww(N,mod-2);
    for(int i=0;i<=N;i++)
    {
        A[i]=((A[i]*invn)%mod+mod)%mod;
    }
}
void NTT(int *FF,int *GG,int n,int m,int *R)
{
    int p2=1,cn=0;
    while(p2<=m+n) p2<<=1,cn++;
    int F[p2+2],G[p2+2];
    for(int i=0;i<=p2;i++) F[i]=G[i]=0;
    for(int i=0;i<=n;i++) F[i]=FF[i];
    for(int i=0;i<=m;i++) G[i]=GG[i];
    Rev(p2,cn);
    FFT(F,p2);
    FFT(G,p2);
    for(int i=0;i<p2;i++)
    {
        F[i]=F[i]*G[i]%mod;
    }
    IFFT(F,p2);
    for(int i=0;i<=n+m;i++)
    {
        R[i]=F[i];
    }
}
signed main()
{
    scanf("%lld%lld",&n,&m);
    for(int i=0;i<=n;i++)
    {
        scanf("%lld",&F[i]);
    }
    for(int i=0;i<=m;i++)
    {
        scanf("%lld",&G[i]);
    }
    NTT(F,G,n,m,F);
    for(int i=0;i<=n+m;i++)
    {
        printf("%lld ",F[i]);
    }
}

多项式求逆

给多项式 F，求 G 满足 F(x) * G(x) \equiv 1 \pmod{x^n}。

容易知道当 n=1 时，G(x)=([0]F)^{-1}。

考虑已知 G\pmod{x^n}，如何求得 G'\pmod{x^{2n}}。

两边乘 F

NTT 优化即可。

注意清空。

#include<iostream>
#include<cstdio>
#include<complex>
#include<cmath>
#include<algorithm>
#define mod 998244353
#define int long long
using namespace std;
int tax[5000001];
const int Gen=3;//原根 
int n,m;
int F[5000001],G[5000001];
int poww(int a,int b)
{
    int ret=1;
    while(b)
    {
        if(b&1)
        {
            ret=ret*a%mod;
        }
        a=a*a%mod;
        b>>=1;
    }
    return ret;
}
void Rev(int n,int cn)
{
    for(int i=1;i<=n;i++)
    {
        tax[i]=0;
        int it=i;
        for(int z=1;z<=cn;z++)
        {
            tax[i]<<=1;
            tax[i]|=(it&1);
            it>>=1;
        }
    }
}
void FFT(int *F,int n)
{
    for(int i=0;i<=n;i++)
    {
        if(tax[i]>i)
        {
            swap(F[i],F[tax[i]]);
        }
    }
    for(int len=2;len<=n;len*=2)
    {
        int W=poww(Gen,(mod-1)/len);
        for(int l=0,r=len-1;r<n;l+=len,r+=len)
        {
            int w=1;
            for(int j=l;j<l+(len>>1);j++)
            {
                int x=F[j],y=F[j+(len>>1)];
                F[j]=(x+y*w)%mod;
                F[j+(len>>1)]=(x-y*w)%mod;
                w=w*W%mod;
            }
        }
    }
}
void IFFT(int *F,int n)
{
    FFT(F,n);
    reverse(F+1,F+n);
    int invn=poww(n,mod-2);
    for(int i=0;i<=n;i++)
    {
        F[i]=F[i]*invn%mod;
    }
}
void INV(int *F,int *G,int n)
{
    if(n==1)
    {
        G[0]=poww(F[0],mod-2);
        return ;
    }
    INV(F,G,n/2+(n%2==1));
    int p2=1,cn=0;
    while(p2<=n+n) p2<<=1,cn++;
    int H[p2+3];
    for(int i=0;i<=p2;i++)
    {
        if(i<n) H[i]=F[i];
        else H[i]=0;
    }
    Rev(p2,cn);
    FFT(H,p2);
    FFT(G,p2);
    for(int i=0;i<=p2;i++)
    {
        G[i]=G[i]*(2-H[i]*G[i]%mod)%mod;    
    }
    IFFT(G,p2);
    for(int i=n;i<=p2;i++)
    {
        G[i]=0; 
    }
}
signed main()
{
    scanf("%lld",&n);
    for(int i=0;i<n;i++)
    {
        scanf("%lld",&F[i]);
    }
    int p2=1;
    while(p2<=n) p2<<=1;
    INV(F,G,n);
    for(int i=0;i<n;i++)
    {
        printf("%lld ",(G[i]+mod)%mod);
    }
}

多项式求导

拆开每一项求导。

合起来。

积反过来即可。

void D(int *F,int n)
{
    for(int i=1;i<=n;i++)
    {
        F[i-1]=F[i]*i%mod;
    }
    F[n]=0;
}
void J(int *F,int n)
{
    for(int i=n;i>0;i--)
    {
        F[i]=F[i-1]*poww(i,mod-2)%mod;
    }
    F[0]=0;
}

多项式 ln

给多项式 F，求 G 满足 G(x) \equiv ln(F(x)) \pmod{x^n}。

算出来再积回去即可。

#include<iostream>
#include<cstdio>
#include<complex>
#include<cmath>
#include<algorithm>
#define int long long
#define mod 998244353
using namespace std;
double pi=acos(-1);
int tax[5000001];
const int Gen=3;//原根 
int n,m;
void Rev(int n,int cn)
{
    for(int i=0;i<n;i++)
    {
        int rw=0,ii=i;
        for(int z=1;z<=cn;z++)
        {
            rw<<=1;
            rw^=ii&1;
            ii>>=1;
        }
        tax[i]=rw;
    }
}
long long poww(int a,int b)
{
    int re=1;
    while(b)
    {
        if(b&1)
        {
            re*=a;
            re%=mod;
        }
        a*=a;
        a%=mod;
        b>>=1;
    }
    return re;
}
void FFT(int *A,int N)
{
    int B[N+2];
    for(int i=0;i<N;i++)
    {
        B[i]=A[i];
    }
    for(int i=1;i<N;i++)
    {
        A[i]=B[tax[i]];
    }
    for(int len=2,M=1;len<=N;M=len,len<<=1)
    {
        int W=poww(Gen,(mod-1)/len);
        for(int l=0,r=len-1;r<=N;l+=len,r+=len)
        {
            int w=1;
            for(int p=l;p<l+M;p++)
            {
                int x=(A[p]+w*A[p+M]%mod)%mod;
                int y=(A[p]-w*A[p+M]%mod)%mod;
                A[p]=x;
                A[p+M]=y;
                w=w*W%mod;
            }
        }
    }
}
void IFFT(int *A,int N)
{
    FFT(A,N);
    reverse(A+1,A+N);
    int invn=poww(N,mod-2);
    for(int i=0;i<=N;i++)
    {
        A[i]=((A[i]*invn)%mod+mod)%mod;
    }
}
long long F[5000001],G[5000001],H[5000001],p2;
void INV(int ci,int *F,int *G)
{
    if(ci==1)
    {
        G[0]=poww(F[0],mod-2);
        return;
    }
    INV((ci+1)>>1,F,G);
    int p2=1,cn=0;
    while(p2<=(ci<<1)) p2<<=1,cn++;
    Rev(p2,cn);
    for(int i=0;i<ci;i++)
    {
        H[i]=F[i];
    }
    for(int i=ci;i<p2;i++)
    {
        H[i]=0;
    }
    FFT(H,p2);
    FFT(G,p2);
    for(int i=0;i<p2;i++)
    {
        G[i]*=(2-H[i]*G[i]%mod+mod)%mod;
        G[i]%=mod; 
    }
    IFFT(G,p2);
    for(int i=ci;i<p2;i++)
    {
        G[i]=0;
    }
}
void D(int *F,int n)
{
    for(int i=1;i<=n;i++)
    {
        F[i-1]=F[i]*i%mod;
    }
    F[n]=0;
}
void J(int *F,int n)
{
    for(int i=n;i>0;i--)
    {
        F[i]=F[i-1]*poww(i,mod-2)%mod;
    }
    F[0]=0;
}
void LN(int *F,int n)
{
    p2=1;
    INV(n,F,G);
    D(F,n-1);
    int p2=1,cn=0;
    while(p2<=n+n) p2*=2,cn++;
    FFT(F,p2);
    FFT(G,p2);
    for(int i=0;i<=p2;i++)
    {
        F[i]=F[i]*G[i]%mod;
    }
    IFFT(F,p2);
    J(F,p2);
}
signed main()
{
    scanf("%lld",&n);
    for(int i=0,_;i<n;i++)
    {
        scanf("%lld",&_);
        F[i]=_;
    }
    LN(F,n);
    for(int i=0;i<n;i++)
    {
        printf("%lld ",(F[i]+mod)%mod);
    }
}

多项式exp

给多项式 F，求 G 满足 G(x) \equiv e^{F(x)} \pmod{x^n}。

牛顿迭代

求函数 G(x) 交 x 轴的点的 x，进行一个迭代。

假如已经求得一个接近点 x_0，我们要让他更接近正确的值。

在 G(x_0) 做切线，交 x 轴于 (x,y) 则

效率至少时 log 的，所以可以递归。

exp

迭代求解。

分治FFT

给 F 求 G 满足 g_i=\sum f_jg_{i-j}

即对于每个 i，他前面的每个数都会有一个贡献。

分治。

• 先计算左边。
• 计算左边对右边的贡献。
• 扔掉左边，计算右边。

难点就在计算左边对右边的贡献。

假如在计算 [l,r]，对于右边一个点 i 它受到的贡献就是 \sum\limits_{j=[1,mid-l+1]} f_{j}g_{i-j}，正好是左边和 f 卷起来，NTT优化即可。