[营业日志]萌新的多项式入门

· · 个人记录

同步更新于caijiのBlog

多项式工业入门

一些定义

template<const int mod>
struct modint{
    int x;
    modint<mod>(int o=0){x=o;}
    modint<mod> &operator = (int o){return x=o,*this;}
    modint<mod> &operator +=(modint<mod> o){return x=x+o.x>=mod?x+o.x-mod:x+o.x,*this;}
    modint<mod> &operator -=(modint<mod> o){return x=x-o.x<0?x-o.x+mod:x-o.x,*this;}
    modint<mod> &operator *=(modint<mod> o){return x=1ll*x*o.x%mod,*this;}
    modint<mod> &operator ^=(int b){
        modint<mod> a=*this,c=1;
        for(;b;b>>=1,a*=a)if(b&1)c*=a;
        return x=c.x,*this;
    }
    modint<mod> &operator /=(modint<mod> o){return *this *=o^=mod-2;}
    modint<mod> &operator +=(int o){return x=x+o>=mod?x+o-mod:x+o,*this;}
    modint<mod> &operator -=(int o){return x=x-o<0?x-o+mod:x-o,*this;}
    modint<mod> &operator *=(int o){return x=1ll*x*o%mod,*this;}
    modint<mod> &operator /=(int o){return *this *= ((modint<mod>(o))^=mod-2);}
    template<class I>friend modint<mod> operator +(modint<mod> a,I b){return a+=b;}
    template<class I>friend modint<mod> operator -(modint<mod> a,I b){return a-=b;}
    template<class I>friend modint<mod> operator *(modint<mod> a,I b){return a*=b;}
    template<class I>friend modint<mod> operator /(modint<mod> a,I b){return a/=b;}
    friend modint<mod> operator ^(modint<mod> a,int b){return a^=b;}
    friend bool operator ==(modint<mod> a,int b){return a.x==b;}
    friend bool operator !=(modint<mod> a,int b){return a.x!=b;}
    bool operator ! () {return !x;}
    modint<mod> operator - () {return x?mod-x:0;}
    modint<mod> &operator++(int){return *this+=1;}
};
const int N=4e6+5;

const int mod=998244353;
const modint<mod> GG=3,Ginv=modint<mod>(1)/3,I=86583718;
struct poly{
    vector<modint<mod>>a;
    modint<mod>&operator[](int i){return a[i];}
    int size(){return a.size();}
    void resize(int n){a.resize(n);}
    void reverse(){std::reverse(a.begin(),a.end());}
};
int rev[N];
inline poly one(){poly a;a.a.push_back(1);return a;}

基础快速变幻

inline int ext(int n){int k=0;while((1<<k)<n)k++;return k;}
inline void init(int k){int n=1<<k;for(int i=0;i<n;i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<(k-1));}
inline void ntt(poly&a,int k,int typ){
    int n=1<<k;
    for(int i=0;i<n;i++)if(i<rev[i])swap(a[i],a[rev[i]]);
    for(int mid=1;mid<n;mid<<=1){
        modint<mod> wn=(typ>0?GG:Ginv)^((mod-1)/(mid<<1));
        for(int r=mid<<1,j=0;j<n;j+=r){
            modint<mod> w=1;
            for(int k=0;k<mid;k++,w=w*wn){
                modint<mod> x=a[j+k],y=w*a[j+k+mid];
                a[j+k]=x+y,a[j+k+mid]=x-y;
            }
        }
    }
    if(typ<0){
        modint<mod> inv=modint<mod>(1)/n;
        for(int i=0;i<n;i++)a[i]*=inv;
    }
}

多项式加、减、乘

poly operator +(poly a,poly b){
    int n=max(a.size(),b.size());a.resize(n),b.resize(n);
    for(int i=0;i<n;i++)a[i]+=b[i];return a;
}
poly operator -(poly a,poly b){
    int n=max(a.size(),b.size());a.resize(n),b.resize(n);
    for(int i=0;i<n;i++)a[i]-=b[i];return a;
}
inline poly operator*(poly a,poly b){
    int n=a.size()+b.size()-1,k=ext(n);
    a.resize(1<<k),b.resize(1<<k),init(k);
    ntt(a,k,1);ntt(b,k,1);for(int i=0;i<(1<<k);i++)a[i]*=b[i];
    ntt(a,k,-1),a.resize(n);return a;
}
inline poly operator*(poly a,modint<mod> b){for(int i=0;i<a.size();i++)a[i]*=b;return a; }
inline poly operator/(poly a,modint<mod> b){for(int i=0;i<a.size();i++)a[i]/=b;return a; }
inline poly operator-(poly a){for(int i=0;i<a.size();i++)a[i]=-a[i];return a; }

多项式求逆

如果多项式F只有一项，那么显然G_0就是F_0的逆元。

若有n项，递归求解。

假如我们已知F(x)H(x) \equiv 1 \pmod{x^{\lceil \frac{n}{2} \rceil}}

又显然F(x)G(x) \equiv 1 \pmod{x^{\lceil \frac{n}{2} \rceil}}

那么F(x)(G(x)-H(x))\equiv 0 \pmod{x^{\lceil \frac{n}{2} \rceil}}

即G(x)-H(x)\equiv 0 \pmod{x^{\lceil \frac{n}{2} \rceil}}

两边同时平方。

G(x)^2+H(x)^2-2G(x)H(x) \equiv 0 \pmod{x^n}

两边同时乘F(x)，再由F(x)G(x) \equiv 1 \pmod{x^n}，得出：

G(x) \equiv 2H(x)-F(x)H(x)^2 \pmod{x^n} 
poly inv(poly F,int k){
    int n=1<<k;F.resize(n);
    if(n==1){F[0]=modint<mod>(1)/F[0];return F;}
    poly G,H=inv(F,k-1);
    G.resize(n),H.resize(n<<1),F.resize(n<<1);
    for(int i=0;i<n/2;i++)G[i]=H[i]*2;
    init(k+1),ntt(H,k+1,1),ntt(F,k+1,1);
    for(int i=0;i<(n<<1);i++)H[i]=H[i]*H[i]*F[i];
    ntt(H,k+1,-1),H.resize(n);
    for(int i=0;i<n;i++)G[i]-=H[i];return G;
}
inline poly inv(poly a){
    int n=a.size();
    a=inv(a,ext(n)),a.resize(n);return a;;
}

求导数、原函数

(x^a)^\prime=ax^{a-1} \int x^adx=\frac{1}{a+1}x^{a+1} 
inline poly deriv(poly a){//求导 
    int n=a.size()-1;
    for(int i=0;i<n;i++)a[i]=a[i+1]*(i+1);
    a.resize(n);return a;
}
inline poly inter(poly a){//求原 
    int n=a.size()+1;a.resize(n);
    for(int i=n;i>=1;i--)a[i]=a[i-1]/i;
    a[0]=0;return a;
}

多项式ln

G(x)=F(A(x)),F(x)=ln(x)

两边同时求导，得到：

G(x)^\prime=F^\prime(A(x))A^\prime(x)=\frac{A^\prime(x)}{A(x)}

在求一遍逆即可。

inline poly ln(poly a){
    int n=a.size();
    a=inter(deriv(a)*inv(a));
    a.resize(n);return a;
}

多项式exp

计算F(x)=e^{A(x)}

变形得\ln F(x)-A(x)=0

设G(F(x))=lnF(x)-A(x)，那么就是要求这一个函数的零点。

那么我们把F(x)看做变量，A(x)看做常数，对这个进行求导，得G^\prime(F(x))=\frac{1}{F(x)}

那么代入牛顿迭代的公式得F(x)\equiv F_0(x)-\frac{G(F_0(x))}{G'(F_0(x))}\pmod{x^n}

F(x)\equiv F_0(x)(1-\ln F_0(x)+A(x))\pmod{x^n}

(这里的F_0(x)是指在\bmod x^{\frac n 2}下的答案)

然后因为A(0)=0，所以F(x)的常数项为1

poly exp(poly a,int k){
    int n=1<<k;a.resize(n);
    if(n==1)return one();
    poly f0=exp(a,k-1);f0.resize(n);
    return f0*(one()+a-ln(f0)); 
}
poly exp(poly a){
    int n=a.size();
    a=exp(a,ext(n));a.resize(n);return a;
}

多项式开根

设H^2(x)\equiv F(x)(\pmod\ x^{\frac n2})

G(x)\equiv H(x)(\pmod\ x^{\frac n2}) G(x)-H(x)\equiv 0(\pmod\ x^{\frac n2}) (G(x)-H(x))^2\equiv 0(\pmod\ x^{\frac n2}) G^2(x)-2H(x)*G(x)+H^2(x)\equiv 0(\pmod\ x^n) F(x)-2H(x)*G(x)+H^2(x)\equiv 0(\pmod\ x^n) G(x)=\frac {F(x)+H^2(x)}{2H(x)}

a_0=1

poly sqrt(poly F,int k){
    int n=1<<k;F.resize(n);
    if(n==1){F[0]=1;return F;}
    poly G,H=sqrt(F,k-1);
    H.resize(n);init(k);poly invH=inv(H);
    G.resize(n),H.resize(n<<1),F.resize(n<<1);
    init(k+1),ntt(H,k+1,1);
    for(int i=0;i<(n<<1);i++)H[i]=H[i]*H[i];
    ntt(H,k+1,-1),H.resize(n);
    for(int i=0;i<n;i++)G[i]=H[i]+F[i];
    G=G*invH;G.resize(n);for(int i=0;i<n;i++)G[i]/=2;
    return G;
}
inline poly sqrt(poly a){
    int n=a.size();
    a=sqrt(a,ext(n)),a.resize(n);return a;;
}

a_0≠1

套用模板中的二次剩余即可

namespace residue{
    //求二次剩余
    int I2;
    struct complex{
        long long real,imag;
        complex&operator=(long long x){real=x,imag=0;}
        complex(long long real=0,long long imag=0):real(real),imag(imag){}
        bool operator==(const complex&b)const{return real==b.real&&imag==b.imag;}
        complex operator*(const complex &b)const{return complex((real*b.real+I2*imag%mod*b.imag)%mod,(real*b.imag+imag*b.real)%mod);}
    };
    complex ksm(complex a,int b){
        complex res=1;
        while(b){
            if(b&1)res=res*a;
            a=a*a;b>>=1;
        }return res;
    }
    bool check(int x){return ksm(complex(x,0),(mod-1)/2).real==1;}
    pair<int,int>solve(int n){
        if(n==0)return {0,0};
        long long a=rand()%mod;
        while(!a||check((a*a+mod-n)%mod))a=rand()%mod;
        I2=(a*a+mod-n)%mod;
        int x0=ksm(complex(a,1),(mod+1)/2).real;
        return {min(x0,mod-x0),max(x0,mod-x0)};
    }
}
poly sqrt(poly F,int k){
    int n=1<<k;F.resize(n);
    if(n==1){F[0]=residue::solve(F[0].x).first;return F;}
    poly G,H=sqrt(F,k-1);
    H.resize(n);init(k);poly invH=inv(H);
    G.resize(n),H.resize(n<<1),F.resize(n<<1);
    init(k+1),ntt(H,k+1,1);
    for(int i=0;i<(n<<1);i++)H[i]=H[i]*H[i];
    ntt(H,k+1,-1),H.resize(n);
    for(int i=0;i<n;i++)G[i]=H[i]+F[i];
    G=G*invH;G.resize(n);for(int i=0;i<n;i++)G[i]/=2;
    return G;
}
inline poly sqrt(poly a){
    int n=a.size();
    a=sqrt(a,ext(n)),a.resize(n);return a;;
}

多项式快速幂

a_0=1

G(x)=(F(x))^k\pmod{x^n} \ln G(x)=kF(x)\pmod{x^n} G(x)=e^{kF(x)}\pmod{x^n} 
inline poly pow(poly a,modint<mod> k){//保证a[0]=1 
    a=ln(a);for(int i=0;i<a.size();i++)a[i]*=k.x;
    return exp(a);
}

a_0≠1

把最低项改为1

F(x)^k=\left(\frac{F(x)}{ax^t} \right)^k(ax)^{tk}

实现了一个modpair，因为k既要对mod取模，也要对\varphi(mod)取模

struct modpair{
    //用于快速幂中的次数
    modint<mod>k1;modint<mod-1>k2;
    struct trueint{
        double lg;int x;
        trueint &operator=(int o){return x=o,lg=log10(o),*this;}
        trueint &operator+=(int o){return lg<=8&&(x+=o,lg=log10(x)),*this;}
        trueint &operator*=(int o){return x*=o,lg+=log10(o),*this;}
    }k3;
    modpair(modint<mod>_1,modint<mod-1>_2,trueint _3):k1(_1),k2(_2),k3(_3){}
    modpair(int o=0){k1=o,k2=o,k3=o;}
    modpair &operator = (int o){return k1=o,k2=o,k3=o,*this;}
    modpair &operator +=(int o){return k1+=o,k2+=o,k3+=o,*this;}
    modpair &operator *=(int o){return k1*=o,k2*=o,k3*=o,*this;}
    friend modpair operator +(modpair a,int o){return a+=o;}
    friend modpair operator *(modpair a,int o){return a*=o;}
    modpair operator-(){return modpair(k1,k2,k3);}
};
inline poly pow2(poly a,modpair m){//不保证a[0]=1 
    int k=0;modint<mod> val;
    while(a[k]==0&&k<a.size())k++;
    if(k==a.size()||k!=0&&m.k3.lg>8||1ll*m.k3.x*k>=a.size()){
    for(int i=0;i<a.size();i++)a[i]=0;return a;}//bye~
    val=a[k];poly b;b.resize(a.size()-k);
    for(int i=0;i<b.size();i++)b[i]=a[i+k]/val;
    b=pow(b,m.k1);for(int i=0;i<a.size();i++)a[i]=0;
    for(int i=0;i<b.size()&&i+k*m.k1.x<a.size();i++)
        a[i+k*m.k1.x]=b[i]*(val^m.k2.x);
    return a;
}

多项式三角函数

由欧拉公式得到：

e^{-iF(x)}=\cos(F(x))-i\sin(F(x)) \end{cases}

解方程得到

\cos F(x)=\frac{e^{iF(x)}+e^{-iF(x)}}{2} \sin F(x)=\frac{e^{iF(x)}-e^{-iF(x)}}{2i}

在取模条件下，i=\sqrt{-1}=\sqrt{mod-1}，即mod-1的二次剩余

inline poly sin(poly a){return (exp(a*I)-exp(-a*I))/(I*2);}
inline poly cos(poly a){return (exp(a*I)+exp(-a*I))/2;}

多项式反三角函数

G(x)\equiv\sin^{-1}F(x)\space (\text{mod }x^n) G'(x)\equiv \frac{F'(x)}{\sqrt{1-F(x)^2}}\space ({\text{mod }x^n}) G(x)\equiv \int \frac{F'(x)}{\sqrt{1-F(x)^2}}\text dx\space ({\text{mod }x^n}) H(x)\equiv\tan^{-1}F(x)\space (\text{mod }x^n) H'(x)\equiv\frac{F'(x)}{1+F(x)^2}\space (\text{mod }x^n) H(x)\equiv\int\frac{F'(x)}{1+F(x)^2} \text dx\space (\text{mod }x^n) 
inline poly asin(poly a){
    poly G=-a*a;G[0]++;
    return inter(deriv(a)*inv(sqrt(G)));
}
inline poly atan(poly a){
    poly G=a*a;G[0]++;
    return inter(deriv(a)*inv(G));
}

完整代码

using namespace std;
template<const int mod>
struct modint{
    int x;
    modint<mod>(int o=0){x=o;}
    modint<mod> &operator = (int o){return x=o,*this;}
    modint<mod> &operator +=(modint<mod> o){return x=x+o.x>=mod?x+o.x-mod:x+o.x,*this;}
    modint<mod> &operator -=(modint<mod> o){return x=x-o.x<0?x-o.x+mod:x-o.x,*this;}
    modint<mod> &operator *=(modint<mod> o){return x=1ll*x*o.x%mod,*this;}
    modint<mod> &operator ^=(int b){
        modint<mod> a=*this,c=1;
        for(;b;b>>=1,a*=a)if(b&1)c*=a;
        return x=c.x,*this;
    }
    modint<mod> &operator /=(modint<mod> o){return *this *=o^=mod-2;}
    modint<mod> &operator +=(int o){return x=x+o>=mod?x+o-mod:x+o,*this;}
    modint<mod> &operator -=(int o){return x=x-o<0?x-o+mod:x-o,*this;}
    modint<mod> &operator *=(int o){return x=1ll*x*o%mod,*this;}
    modint<mod> &operator /=(int o){return *this *= ((modint<mod>(o))^=mod-2);}
    template<class I>friend modint<mod> operator +(modint<mod> a,I b){return a+=b;}
    template<class I>friend modint<mod> operator -(modint<mod> a,I b){return a-=b;}
    template<class I>friend modint<mod> operator *(modint<mod> a,I b){return a*=b;}
    template<class I>friend modint<mod> operator /(modint<mod> a,I b){return a/=b;}
    friend modint<mod> operator ^(modint<mod> a,int b){return a^=b;}
    friend bool operator ==(modint<mod> a,int b){return a.x==b;}
    friend bool operator !=(modint<mod> a,int b){return a.x!=b;}
    bool operator ! () {return !x;}
    modint<mod> operator - () {return x?mod-x:0;}
    modint<mod> &operator++(int){return *this+=1;}
};
const int N=4e6+5;

const int mod=998244353;
namespace residue{
    //求二次剩余
    int I2;
    struct complex{
        long long real,imag;
        complex&operator=(long long x){real=x,imag=0;}
        complex(long long real=0,long long imag=0):real(real),imag(imag){}
        bool operator==(const complex&b)const{return real==b.real&&imag==b.imag;}
        complex operator*(const complex &b)const{return complex((real*b.real+I2*imag%mod*b.imag)%mod,(real*b.imag+imag*b.real)%mod);}
    };
    complex ksm(complex a,int b){
        complex res=1;
        while(b){
            if(b&1)res=res*a;
            a=a*a;b>>=1;
        }return res;
    }
    bool check(int x){return ksm(complex(x,0),(mod-1)/2).real==1;}
    pair<int,int>solve(int n){
        if(n==0)return {0,0};
        long long a=rand()%mod;
        while(!a||check((a*a+mod-n)%mod))a=rand()%mod;
        I2=(a*a+mod-n)%mod;
        int x0=ksm(complex(a,1),(mod+1)/2).real;
        return {min(x0,mod-x0),max(x0,mod-x0)};
    }
}
const modint<mod> GG=3,Ginv=modint<mod>(1)/3,I=86583718;
struct poly{
    vector<modint<mod>>a;
    modint<mod>&operator[](int i){return a[i];}
    int size(){return a.size();}
    void resize(int n){a.resize(n);}
    void reverse(){std::reverse(a.begin(),a.end());}
};
int rev[N];
inline int ext(int n){int k=0;while((1<<k)<n)k++;return k;}
inline void init(int k){int n=1<<k;for(int i=0;i<n;i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<(k-1));}
inline void ntt(poly&a,int k,int typ){
    int n=1<<k;
    for(int i=0;i<n;i++)if(i<rev[i])swap(a[i],a[rev[i]]);
    for(int mid=1;mid<n;mid<<=1){
        modint<mod> wn=(typ>0?GG:Ginv)^((mod-1)/(mid<<1));
        for(int r=mid<<1,j=0;j<n;j+=r){
            modint<mod> w=1;
            for(int k=0;k<mid;k++,w=w*wn){
                modint<mod> x=a[j+k],y=w*a[j+k+mid];
                a[j+k]=x+y,a[j+k+mid]=x-y;
            }
        }
    }
    if(typ<0){
        modint<mod> inv=modint<mod>(1)/n;
        for(int i=0;i<n-1;i++)a[i]*=inv;
    }
}
inline poly one(){poly a;a.a.push_back(1);return a;}
poly operator +(poly a,poly b){
    int n=max(a.size(),b.size());a.resize(n),b.resize(n);
    for(int i=0;i<n;i++)a[i]+=b[i];return a;
}
poly operator -(poly a,poly b){
    int n=max(a.size(),b.size());a.resize(n),b.resize(n);
    for(int i=0;i<n;i++)a[i]-=b[i];return a;
}
inline poly operator*(poly a,poly b){
    int n=a.size()+b.size()-1,k=ext(n);
    a.resize(1<<k),b.resize(1<<k),init(k);
    ntt(a,k,1);ntt(b,k,1);for(int i=0;i<(1<<k);i++)a[i]*=b[i];
    ntt(a,k,-1),a.resize(n);return a;
}
inline poly operator*(poly a,modint<mod> b){for(int i=0;i<a.size();i++)a[i]*=b;return a; }
inline poly operator/(poly a,modint<mod> b){for(int i=0;i<a.size();i++)a[i]/=b;return a; }
inline poly operator-(poly a){for(int i=0;i<a.size();i++)a[i]=-a[i];return a; }
poly inv(poly F,int k){
    int n=1<<k;F.resize(n);
    if(n==1){F[0]=modint<mod>(1)/F[0];return F;}
    poly G,H=inv(F,k-1);
    G.resize(n),H.resize(n<<1),F.resize(n<<1);
    for(int i=0;i<n/2;i++)G[i]=H[i]*2;
    init(k+1),ntt(H,k+1,1),ntt(F,k+1,1);
    for(int i=0;i<(n<<1);i++)H[i]=H[i]*H[i]*F[i];
    ntt(H,k+1,-1),H.resize(n);
    for(int i=0;i<n;i++)G[i]-=H[i];return G;
}
inline poly inv(poly a){
    int n=a.size();
    a=inv(a,ext(n)),a.resize(n);return a;;
}
inline poly deriv(poly a){//求导 
    int n=a.size()-1;
    for(int i=0;i<n;i++)a[i]=a[i+1]*(i+1);
    a.resize(n);return a;
}
inline poly inter(poly a){//求原 
    int n=a.size()+1;a.resize(n);
    for(int i=n;i>=1;i--)a[i]=a[i-1]/i;
    a[0]=0;return a;
}
inline poly ln(poly a){
    int n=a.size();
    a=inter(deriv(a)*inv(a));
    a.resize(n);return a;
}
poly exp(poly a,int k){
    int n=1<<k;a.resize(n);
    if(n==1)return one();
    poly f0=exp(a,k-1);f0.resize(n);
    return f0*(one()+a-ln(f0)); 
}
poly exp(poly a){
    int n=a.size();
    a=exp(a,ext(n));a.resize(n);return a;
}
poly sqrt(poly F,int k){
    int n=1<<k;F.resize(n);
    if(n==1){F[0]=residue::solve(F[0].x).first;return F;}
    poly G,H=sqrt(F,k-1);
    H.resize(n);init(k);poly invH=inv(H);
    G.resize(n),H.resize(n<<1),F.resize(n<<1);
    init(k+1),ntt(H,k+1,1);
    for(int i=0;i<(n<<1);i++)H[i]=H[i]*H[i];
    ntt(H,k+1,-1),H.resize(n);
    for(int i=0;i<n;i++)G[i]=H[i]+F[i];
    G=G*invH;G.resize(n);for(int i=0;i<n;i++)G[i]/=2;
    return G;
}
inline poly sqrt(poly a){
    int n=a.size();
    a=sqrt(a,ext(n)),a.resize(n);return a;;
}
inline poly pow(poly a,modint<mod> k){//保证a[0]=1 
    a=ln(a);for(int i=0;i<a.size();i++)a[i]*=k.x;
    return exp(a);
}

struct modpair{
    //用于快速幂中的次数
    modint<mod>k1;modint<mod-1>k2;
    struct trueint{
        double lg;int x;
        trueint &operator=(int o){return x=o,lg=log10(o),*this;}
        trueint &operator+=(int o){return lg<=8&&(x+=o,lg=log10(x)),*this;}
        trueint &operator*=(int o){return x*=o,lg+=log10(o),*this;}
    }k3;
    modpair(modint<mod>_1,modint<mod-1>_2,trueint _3):k1(_1),k2(_2),k3(_3){}
    modpair(int o=0){k1=o,k2=o,k3=o;}
    modpair &operator = (int o){return k1=o,k2=o,k3=o,*this;}
    modpair &operator +=(int o){return k1+=o,k2+=o,k3+=o,*this;}
    modpair &operator *=(int o){return k1*=o,k2*=o,k3*=o,*this;}
    friend modpair operator +(modpair a,int o){return a+=o;}
    friend modpair operator *(modpair a,int o){return a*=o;}
    modpair operator-(){return modpair(k1,k2,k3);}
};
inline poly pow2(poly a,modpair m){//不保证a[0]=1 
    int k=0;modint<mod> val;
    while(a[k]==0&&k<a.size())k++;
    if(k==a.size()||k!=0&&m.k3.lg>8||1ll*m.k3.x*k>=a.size()){
    for(int i=0;i<a.size();i++)a[i]=0;return a;}//bye~
    val=a[k];poly b;b.resize(a.size()-k);
    for(int i=0;i<b.size();i++)b[i]=a[i+k]/val;
    b=pow(b,m.k1);for(int i=0;i<a.size();i++)a[i]=0;
    for(int i=0;i<b.size()&&i+k*m.k1.x<a.size();i++)
        a[i+k*m.k1.x]=b[i]*(val^m.k2.x);
    return a;
}
inline poly sin(poly a){return (exp(a*I)-exp(-a*I))/(I*2);}
inline poly cos(poly a){return (exp(a*I)+exp(-a*I))/2;}
inline poly asin(poly a){
    poly G=-a*a;G[0]++;
    return inter(deriv(a)*inv(sqrt(G)));
}
inline poly atan(poly a){
    poly G=a*a;G[0]++;
    return inter(deriv(a)*inv(G));
}

多项式黑科技

任意模数NTT(MTT)

把两个多项式 A(x),B(x) 分别拆成 A_0(x),A_1(x),B_0(x),B_1(x)后，求出它们的点值表示。

构造两个多项式：

P(x)=A(x)+iB(x) Q(x)=A(x)+iB(x)

由于A,B的虚部都为0，P,Q的每一项系数都互为共轭。对P做一遍DFT,可以std::conj在O(n)求出Q的点值表示。然后A_0,A_1,B_0,B_1的点值就是有手就行了

接下求A,B之间的两两乘积，乘出来对A_0B_0,A_1B_0,A_0B_0,A_0B_1做IDFT。

由于虚部不为0，但仍然可以借用上述的方法。构造两个多项式：

P(x)=A_0(x)B_0(x)+iA_1(x)B_0(x) Q(X)=A_0(x)B_1(x)+iA_1(x)B_1(x)

分别做IDFT，由于 A_0B_0,A_0B_1,A_1B_0,A_1B_1这四个多项式卷起来后的系数表示中虚部一定为0，那么P的实部与虚部就分别为A_0(x)B_0(x)和A_1(x)B_0(x)，而Q就位A_0(x)B_1(x)与A_1(x)B_1(x)

给出部分代码：

typedef std::complex<double>complex;
const int N=4e6+10;const double PI=acos(-1);const complex I=complex(0,1);
int rev[N];complex Wn[N];int M,mod;
int ksm(int x,int y) {
    int re=1;
    for(;y;y>>=1,x=1LL*x*x%mod)if(y&1)re=1LL*re*x%mod;
    return re;
}
inline long long num(complex x){double d=x.real();return d<0?(long long)(d-0.5)%mod:(long long)(d+0.5)%mod;}
struct poly{
    std::vector<complex>a0,a1;
    int size(){return a0.size();}
    void resize(int n){a0.resize(n);a1.resize(n);}
    void set(int x,long long y){
        y%=mod;
        a0[x]=y/M;
        a1[x]=y%M;
    }
    long long get(int x){return (M*M*num(a0[x].real())%mod +
                M*(num(a0[x].imag())+num(a1[x].real()))%mod+num(a1[x].imag()))%mod;}
    long long val(int x){return (long long)(M*a0[x].real()+a1[x].real()+mod)%mod;}
};
poly operator+(poly a,poly b){
    int n=std::max(a.size(),b.size());a.resize(n),b.resize(n);
    for(int i=0;i<n;i++)a.set(i,a.val(i)+b.val(i));return a;
}
poly operator-(poly a,poly b){
    int n=std::max(a.size(),b.size());a.resize(n),b.resize(n);
    for(int i=0;i<n;i++)a.set(i,a.val(i)-b.val(i));return a;
}
inline poly one(){poly a;a.resize(1);a.set(0,1);return a;}
inline int ext(int n){int k=0;while((1<<k)<n)k++;return k;}
inline void init(int k){
    int n=1<<k;
    for(int i=0;i<n;i++)
        rev[i]=(rev[i>>1]>>1)|((i&1)<<(k-1));
    for(int i=0;i<n;i++)
        Wn[i]={cos(PI/n*i),sin(PI/n*i)};
}
void FFT(std::vector<complex>&A,int n,int t){
    if(t<0)for(int i=1;i<n;i++)if(i<(n-i))std::swap(A[i],A[n-i]);
    for(int i=0;i<n;i++)
        if(i<rev[i])std::swap(A[i],A[rev[i]]);
    for(int m=1;m<n;m<<=1)
        for(int i=0;i<n;i+=m<<1)
            for(int k=i;k<i+m;k++){
                complex W=Wn[1ll*(k-i)*n/m];
                complex a0=A[k],a1=A[k+m]*W;
                A[k]=a0+a1;A[k+m]=a0-a1; 
            }
    if(t<0)for(int i=0;i<n;i++)A[i]/=n;
}
void MTT(poly &A,int n,int t){
    for(int i=0;i<n;i++)A.a0[i]=A.a0[i]+I*A.a1[i];
    FFT(A.a0,n,t);
    for(int i=0;i<n;i++)A.a1[i]=std::conj(A.a0[i?n-i:0]);
    for(int i=0;i<n;i++){
        complex p=A.a0[i],q=A.a1[i];
        A.a0[i]=(p+q)*0.5;A.a1[i]=(q-p)*0.5*I;
    }
}
inline poly operator*(poly a,poly b){
    int n=a.size()+b.size()-1,k=ext(n);
    a.resize(1<<k),b.resize(1<<k),init(k);
    MTT(a,1<<k,1);MTT(b,1<<k,1);
    for(int i=0;i<(1<<k);i++){
        complex p=a.a0[i]*b.a0[i]+I*a.a1[i]*b.a0[i];
        complex q=a.a0[i]*b.a1[i]+I*a.a1[i]*b.a1[i];
        a.a0[i]=p,a.a1[i]=q;
    }
    FFT(a.a0,1<<k,-1);FFT(a.a1,1<<k,-1);a.resize(n);
    long long tmp;for(int i=0;i<n;i++)
        tmp=a.get(i),a.set(i,tmp);
    return a;
}
inline poly deriv(poly a){//求导 
    int n=a.size()-1;
    for(int i=0;i<n;i++)a.set(i,a.val(i+1)*(i+1));
    a.resize(n);return a;
}
inline poly inter(poly a){//求原 
    int n=a.size()+1;a.resize(n);
    for(int i=n;i>=1;i--)a.set(i,a.val(i-1)*ksm(i,mod-2));
    a.set(0,0);return a;
}
poly inv(poly F,int k){
    int n=1<<k;F.resize(n);
    if(n==1){F.set(0,ksm(F.val(0),mod-2));return F;}
    poly G,H=inv(F,k-1);
    G.resize(n),H.resize(n<<1),F.resize(n<<1);
    for(int i=0;i<n/2;i++)G.set(i,H.val(i)*2);
    H=H*H;H.resize(n);H=H*F;H.resize(n);
    G=G-H;return G;
}
inline poly inv(poly a){
    int n=a.size();
    a=inv(a,ext(n)),a.resize(n);return a;;
}
inline poly ln(poly a){
    int n=a.size();
    a=inter(deriv(a)*inv(a));
    a.resize(n);return a;
}

别问我vector太慢怎么办下次一定改

多项式除法

太神奇了蒟蒻只能膜拜

设多项式A为n次多项式，考虑一种操作R，使得

A_R(x)=x^nA(\frac 1x)

不难发现A_R就是对A进行系数翻转。

F(x)=Q(x)\times G(x)+R(x) F(\frac 1x)=Q(\frac 1x)\times G(\frac 1x)+R(\frac 1x) x^nF(\frac 1x)=x^{n-m}Q(\frac 1x)\times x^mG(\frac 1x)+x^{n-m+1}\times x^{m-1}R(\frac 1x) F_R(x)=Q_R(x)\times G_R(x)+x^{n-m+1}\times R_R(x) F_R(x)\equiv Q_R(x)\times G_R(x) \pmod {x^{n-m+1}} Q_R(x)\equiv F_R(x)\times G_R^{-1}(x) \pmod {x^{n-m+1}}

然后R就有手就行了

R(x)=F(x)-G(x)\times Q(x) 
inline pair<poly,poly>div(poly rF,poly rG){
    poly Q,R;int q=rF.size()-rG.size()+1,n=rF.size(),m=rG.size();poly F=rF,G=rG;
    rF.reverse();rG.reverse();rF.resize(q);rG.resize(q);
    Q=rF*inv(rG);Q.resize(q);Q.reverse();
    R=F-(G*Q);R.resize(m-1);
    return {Q,R};
}