exp等)

(latest updated on 2021.02.23)

前置知识NTT

所有操作均在对$P=\text{998244353}$取模下进行

代码在最下面，由于板子实在有一点长，所以。。。

下文中$\pmod {x^n}$表示求出了多项式的前$n$项

$[x^i]F(x)$表示$F(x)$第$i$项的系数

每个小问题的模板题都可以在洛谷上找到

\[\ \]

1.多项式求乘法逆

~~(为什么叫做乘法逆?因为还有求复合逆和模逆元的)~~

求 $G(x)\equiv \frac{1}{F(x)} \pmod {x^n}$

形象化的理解就是$F(x)\cdot G(x) \pmod {x^n}$只有第一项是$1$，其他项都是$0$

这个由于是第一个操作，很多人还并不是很能理解多项式操作到底是什么东西，所以讲多一点

Part1 O($n^2$)

为了便于理解这个问题，先考虑一个最简单的模拟

$[x^i]F\cdot G(x)=\sum [x^j]F(x)[x^{i-j}]G(x)$

第一项$[x^0]G(x)=\frac{1}{[x^0]F(0)} \pmod P$，因此求逆的前提条件是$[x^0]F(x)\ne 0$

考虑从$1$到$n-1$依次求出每一项，先从前面的项中得到所有$j>0$的和$Sum$，然后带入$j=0$时知道

\[[x^i]G(x)=-\frac{Sum=\sum_{j=1}^{j\leq i}[x^j]F(x)[x^{i-j}]G(x)}{[x^0]F(0)} \]

\[\ \]

Part2 O($n\log^2n$)

上面这个式子是一个类似$dp$转移的东西，可以直接分治NTT优化掉

\[\ \]

Part3 $O(n\log n)$

考虑倍增求解，设已经求出了

\[H(x)\equiv \frac{1}{F(x)},\pmod {x^{\frac{n}{2}}} \]

其中递归边界是$n=1$时,$[x^0]G(x)=\frac{1}{[x^0]F(0)} \pmod P$，因此求逆的前提条件是$[x^0]F(x)\ne 0$

则

\[H(x)\equiv G(x)\pmod {x^{\frac{n}{2}}} \]

\[H(x)-G(x)\equiv 0\pmod {x^{\frac{n}{2}}} \]

我们对于$H(x)-G(x)$平方，结果的前$n$项不可能由两个$\ge \frac{n}{2}$的项相乘得到，而前$\frac{n}{2}$项都是$0$，所以

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

所以通过平方可以扩大模数，这很常用

展开平方的式子

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

两边乘上$F(x)$

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

\[G(x)\equiv 2H(x)-H(x)^2F(x)\pmod {x^n} \]

带入这个式子倍增求解即可

分析复杂度，每次有一个$H(x)^2F(x)$，可以通过$NTT$求出，倍增过程中访问的长度是$O(n+\frac{n}{2}+\frac{n}{4}...)=O(n)$

所以总复杂度就是$O(n\log n)$

\[\ \]

2.多项式开根号

求$G(x)^2\equiv F(x) \pmod {x^n}$

同样的，递归求解，边界是$[x^0]=\sqrt{[x^0]F(x)} \pmod P$

可以发现我们需要求二次剩余。。。但是一般题目保证了$[x^0]F(x)\in\{0,1\}$

设已经求出$H(x)^2\equiv F(x) \pmod{ x^{\lceil \frac{n}{2} \rceil}}$

\[H(x)\equiv G(x) \pmod {x^{\lceil \frac{n}{2}\rceil}} \]

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

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

\[G(x)\equiv \frac{H(x)^2+F(x)}{2H(x)} \pmod {x^n} \]

带入这个式子倍增求解即可

复杂度为$O(n\log n)$，由于需要求逆，实际可能会比较难写

\[\ \]

3.多项式求$\ln$

对$\begin{aligned} G(x)\equiv \ln F(x) \pmod {x^n} \end{aligned} $ 两边求导，注意这里是复合函数求导!!!

$\begin{aligned} G'(x)\equiv F'(x)\frac{1}{F(x)} \pmod {x^n}\end{aligned}$

求出$G'(x)$，然后求原函数即可

通常保证$[x^0]F(x)=1$，否则不好求$\ln 1$，所以求出原函数后首项为0

复杂度为$O(n\log n)$

\[\ \]

4.多项式求exp

多项式求$\text{exp}$即求$G(x)=e^{F(x)} \mod x^n$

多项式求$\text{exp}$常见的解法有两种

CDQ分治+$\text{NTT}$

要求$G(x)=e^{F(x)}$

式子两边求导(右边要复合函数求导)，$G'(x)=F'(x) e^{F(x)}$

也就是说，$G'(x)=F'(x)G(x)$

两边同时积分得到$\begin{aligned} G(x)=\int{F'(x)G(x)}\end{aligned}$

我们知道，$ [x^i] \begin{aligned}\int H(x) =\frac{ [x^{i-1}]H(x)}{i}\end{aligned} $

带入上面的式子得到$\displaystyle i\cdot [x^i]G(x)= \sum_{j=0}^{i-1}[x^j]F'(x)\cdot [x^{i-1-j}]G(x)$

那么对于这个式子，直接使用分治NTT求解，其复杂为$O(n\log n)$

\[\ \]

牛顿迭代

这是一种渐进意义上更优的做法，但实际在$10^6$以下几乎不可能更快，而且代码难写

但是不管平时用不用，牛顿迭代的知识学习一下肯定是最好的

把题目转化为，对于函数$f(G)=\ln G-F$

求出在$\mod x^n$意义下的零点

其中$f(x)=\ln x-c$

考虑迭代求解，设已经求出$H(x)=e^{F(x)}\pmod {x^{\frac{n}{2}}}$

边界条件是$[x^0]H(x)=e^{[x^0]F(x)}$(由于没有办法求$e^x$在模意义下的值，所以通常必须要满足$[x^0]F(x)=0$)

带入牛顿迭代的结果

\[G=H-\frac{f(H)}{f'(H)}=H(F-\ln H+1) \]

每次求$\ln$ 复杂度和$\text{NTT}$相同，所以总复杂度为$O(n\log n)$

事实上这个还有优化的余地，就是在求$\ln$的时候，多项式逆的部分可以同步倍增求出，不需要每次都倍增一下~~(但是好像效果并不是特别明显)~~

\[\ \]

5.多项式$k$次幂

$G(x)\equiv F(x)^k\pmod {x^n}$

$\ln G(x)=k \ln F(x) \pmod {x^n}$

求出$\ln G(x)$之后，$\exp$回来即可

由于要求$\ln$，所以这样求的条件是$[x^0]F(x)=1$ (可以通过平移和系数变换来调整为$[x^0]F(x)=1$)

~~很显然这个方法对于开根号也是适用的~~

复杂度$O(n\log n)$

\[\ \]

6.多项式带余除法

问题：给定$F(x),G(x)$，其次数为$n,m,n>m$

求$F(x)=G(x)P(x)+R(x)$，其中$P(x)$次数为$n-m$，$R(x)$次数为$m-1$

考虑先求解$P(x)$，下面引入一种翻转运算

$F^R(x)=x^nF(\frac{1}{x})$，即将$F(x)$的系数翻转排列

用$\frac{1}{x}$带入问题的式子，得到

$\displaystyle F(\frac{1}{x})=G(\frac{1}{x})P(\frac{1}{x})+R(x)$

$\displaystyle x^nF(\frac{1}{x})=x^mG(\frac{1}{x})\cdot x^{n-m}P(\frac{1}{x})+x^nR(x)$

$\displaystyle F^R(x)=G^R(x)\cdot P^R(x)+x^{n-m+1}R^R(x)$

要求的$P(x)$是$n-m$次的，所以$R^R(x)\cdot x^{n-m+1}$并没有贡献

此时可以认为$\displaystyle P^R(x)=\frac{F^R(x)}{G^R(x)}$，求逆即可得到

得到$P(x)$之后，带入$R(x)=F(x)-G(x)P(x)$即可

\[\ \]

应用：多项式多点求值，常系数线性齐次递推

\[\ \]

以上是基本运算，如果不想继续~~吸多项式~~请直接跳到最下面的代码

多项式与点值式

下降幂多项式初步

\[\ \]

所有的操作均用$\text{vector} $来实现，主要是为了理清思路，并且清零问题上会比较容易解决，同时对于每次计算完多项式的长度的要求会显得更加严格

实际在UOJ/Luogu上会非常慢，在LOJ上不错

稍微整理了一下，没怎么卡过常，所以应该还是比较可读的

代码总览(请使用C++11,O2编译运行)

基本运算的总模板题Loj - 150

#include<bits/stdc++.h>
using namespace std;

typedef long long ll;
typedef double db;
typedef unsigned long long ull;
typedef pair <int,int> Pii;
#define reg register
#define pb push_back
#define mp make_pair
#define Mod1(x) ((x>=P)&&(x-=P))
#define Mod2(x) ((x<0)&&(x+=P))
#define rep(i,a,b) for(int i=a,i##end=b;i<=i##end;++i)
#define drep(i,a,b) for(int i=a,i##end=b;i>=i##end;--i)
template <class T> inline void cmin(T &a,T b){ ((a>b)&&(a=b)); }
template <class T> inline void cmax(T &a,T b){ ((a<b)&&(a=b)); }

char IO;
template <class T=int> T rd(){
    T s=0; int f=0;
    while(!isdigit(IO=getchar())) if(IO=='-') f=1;
    do s=(s<<1)+(s<<3)+(IO^'0');
    while(isdigit(IO=getchar()));
    return f?-s:s;
}

const int N=1<<17,P=998244353;

int n,k;

ll qpow(ll x,ll k=P-2) {
    ll res=1;
    for(;k;k>>=1,x=x*x%P) if(k&1) res=res*x%P;
    return res;
}

/*
void NTT(int n,int *a,int f){
    rep(i,1,n-1) if(i<rev[i]) swap(a[i],a[rev[i]]);
    for(int i=1;i<n;i<<=1) {
        int w=qpow(3,(P-1)/i/2);
        for(int l=0;l<n;l+=i*2) {
            int e=1;
            for(int j=l;j<l+i;++j,e=1ll*e*w%P) {
                int t=1ll*a[j+i]*e%P;
                a[j+i]=a[j]-t,((a[j+i]<0)&&(a[j+i]+=P));
                a[j]+=t,((a[j]>=P)&&(a[j]-=P));
            }
        }
    }
    if(f==-1) {
        reverse(a+1,a+n);
        int Inv=qpow(n);
        rep(i,0,n-1) a[i]=1ll*a[i]*Inv%P;
    }
}

int e[N];
void NTT(int n,int *a,int f){
    rep(i,1,n-1) if(i<rev[i]) swap(a[i],a[rev[i]]);
    e[0]=1;
    for(int i=1;i<n;i<<=1) {
        int w=qpow(3,(P-1)/i/2);
        for(int j=1;j<i;++j) e[j]=1ll*e[j-1]*w%P;
        for(int l=0;l<n;l+=i*2) {
            for(int j=l;j<l+i;++j) {
                int t=1ll*a[j+i]*e[j-l]%P;
                a[j+i]=a[j]-t,((a[j+i]<0)&&(a[j+i]+=P));
                a[j]+=t,((a[j]>=P)&&(a[j]-=P));
            }
        }
    }
    if(f==-1) {
        reverse(a+1,a+n);
        int Inv=qpow(n);
        rep(i,0,n-1) a[i]=1ll*a[i]*Inv%P;
    }
}

int e[N];
void NTT(int n,int *a,int f){
    rep(i,1,n-1) if(i<rev[i]) swap(a[i],a[rev[i]]);
    e[0]=1;
    for(int i=1;i<n;i<<=1) {
        int w=qpow(3,(P-1)/i/2);
        for(int j=i-2;j>=0;j-=2) e[j+1]=1ll*w*(e[j]=e[j>>1])%P;
        for(int l=0;l<n;l+=i*2) {
            for(int j=l;j<l+i;++j) {
                int t=1ll*a[j+i]*e[j-l]%P;
                a[j+i]=a[j]-t,((a[j+i]<0)&&(a[j+i]+=P));
                a[j]+=t,((a[j]>=P)&&(a[j]-=P));
            }
        }
    }
    if(f==-1) {
        reverse(a+1,a+n);
        int Inv=qpow(n);
        rep(i,0,n-1) a[i]=1ll*a[i]*Inv%P;
    }
}

int w[N];
void Init(int N){
    w[N>>1]=1;
    int t=qpow(3,(P-1)/N);
    rep(i,(N>>1)+1,N-1) w[i]=1ll*w[i-1]*t%P;
    drep(i,(N>>1)-1,1) w[i]=w[i<<1];
}
void NTT(int n,int *a,int f){
    rep(i,1,n-1) if(i<rev[i]) swap(a[i],a[rev[i]]);
    for(int i=1;i<n;i<<=1) {
        int *e=w+i;
        for(int l=0;l<n;l+=i*2) {
            for(int j=l;j<l+i;++j) {
                int t=1ll*a[j+i]*e[j-l]%P;
                a[j+i]=a[j]-t,((a[j+i]<0)&&(a[j+i]+=P));
                a[j]+=t,((a[j]>=P)&&(a[j]-=P));
            }
        }
    }
    if(f==-1) {
        reverse(a+1,a+n);
        int Inv=qpow(n);
        rep(i,0,n-1) a[i]=1ll*a[i]*Inv%P;
    }
}
*/

/*
namespace MTT{
    const double PI=acos((double)-1);
    int rev[N];
    struct Cp{
        double x,y;
        Cp(){ ; }
        Cp(double _x,double _y): x(_x),y(_y){ } 
        inline Cp operator + (const Cp &t) const { return (Cp){x+t.x,y+t.y}; }
        inline Cp operator - (const Cp &t) const { return (Cp){x-t.x,y-t.y}; }
        inline Cp operator * (const Cp &t) const { return (Cp){x*t.x-y*t.y,x*t.y+y*t.x}; }
    }A[N],B[N],C[N],w[N/2];
#define E(x) ll(x+0.5)%P
    void FFT(int n,Cp *a,int f){
        rep(i,0,n-1) if(rev[i]<i) swap(a[i],a[rev[i]]);
        w[0]=Cp(1,0);
        for(reg int i=1;i<n;i<<=1) {
            Cp t=Cp(cos(PI/i),f*sin(PI/i));
            for(reg int j=i-2;j>=0;j-=2) w[j+1]=t*(w[j]=w[j>>1]);
            // 上面提到的最优板子
            for(reg int l=0;l<n;l+=2*i) {
                for(reg int j=l;j<l+i;j++) {
                    Cp t=a[j+i]*w[j-l];
                    a[j+i]=a[j]-t;
                    a[j]=a[j]+t;
                }
            }
        }
        if(f==-1) rep(i,0,n-1) a[i].x/=n,a[i].y/=n;
    }
    void Multiply(int n,int m,int *a,int *b,int *res,int P){
        // [0,n-1]*[0,m-1]->[0,n+m-2]
        int S=(1<<15)-1;
        int R=1,cc=-1;
        while(R<=n+m-1) R<<=1,cc++;
        rep(i,1,R) rev[i]=(rev[i>>1]>>1)|((i&1)<<cc);
        rep(i,0,n-1) A[i]=Cp((a[i]&S),(a[i]>>15));
        rep(i,0,m-1) B[i]=Cp((b[i]&S),(b[i]>>15));
        rep(i,n,R-1) A[i]=Cp(0,0);
        rep(i,m,R-1) B[i]=Cp(0,0);
        FFT(R,A,1),FFT(R,B,1);
        rep(i,0,R-1) {
            int j=(R-i)%R;
            C[i]=Cp((A[i].x+A[j].x)/2,(A[i].y-A[j].y)/2)*B[i];
            B[i]=Cp((A[i].y+A[j].y)/2,(A[j].x-A[i].x)/2)*B[i];
        }
        FFT(R,C,-1),FFT(R,B,-1);
        rep(i,0,n+m-2) {
            ll a=E(C[i].x),b=E(C[i].y),c=E(B[i].x),d=E(B[i].y);
            res[i]=(a+((b+c)<<15)+(d<<30))%P;
        }
    }
#undef E
}
*/


namespace Polynomial{

    typedef vector <int> Poly;
    void Show(Poly a,int k=0){ 
        if(!k){ for(int i:a) printf("%d ",i); puts(""); }
        else for(int i:a) printf("%d\n",i);
    }
    int rev[N],w[N];
    int Inv[N+1],Fac[N+1],FInv[N+1];

    void Init_w() { 
        int t=qpow(3,(P-1)/N);
        w[N>>1]=1;
        rep(i,(N>>1)+1,N-1) w[i]=1ll*w[i-1]*t%P;
        drep(i,(N>>1)-1,1) w[i]=w[i<<1];
        Inv[0]=Inv[1]=Fac[0]=Fac[1]=FInv[0]=FInv[1]=1;
        rep(i,2,N) {
            Inv[i]=1ll*(P-P/i)*Inv[P%i]%P; 
            FInv[i]=1ll*FInv[i-1]*Inv[i]%P;
            Fac[i]=1ll*Fac[i-1]*i%P;
        }
    }
    int Init(int n){
        int R=1,c=-1;
        while(R<n) R<<=1,c++;
        rep(i,1,R-1) rev[i]=(rev[i>>1]>>1)|((i&1)<<c);
        return R;
    }

#define NTTVersion1

#ifdef NTTVersion1
    void NTT(int n,Poly &a,int f){
        rep(i,0,n-1) if(rev[i]<i) swap(a[i],a[rev[i]]);
        for(int i=1;i<n;i<<=1) {
            int *e=w+i;
            for(int l=0;l<n;l+=i*2) {
                for(int j=l;j<l+i;++j) {
                    int t=1ll*a[j+i]*e[j-l]%P;
                    a[j+i]=a[j]-t,Mod2(a[j+i]);
                    a[j]+=t,Mod1(a[j]);
                }
            }
        }
        if(f==-1) {
            reverse(a.begin()+1,a.begin()+n);
            ll base=Inv[n];
            rep(i,0,n-1) a[i]=a[i]*base%P;
        }
    }
    void NTT(int n,int *a,int f){
        rep(i,0,n-1) if(rev[i]<i) swap(a[i],a[rev[i]]);
        for(int i=1;i<n;i<<=1) {
            int *e=w+i;
            for(int l=0;l<n;l+=i*2) {
                for(int j=l;j<l+i;++j) {
                    int t=1ll*a[j+i]*e[j-l]%P;
                    a[j+i]=a[j]-t,Mod2(a[j+i]);
                    a[j]+=t,Mod1(a[j]);
                }
            }
        }
        if(f==-1) {
            reverse(a+1,a+n);
            ll base=Inv[n];
            rep(i,0,n-1) a[i]=a[i]*base%P;
        }
    }

#else 

    void NTT(int n,Poly &a,int f){ 
        static int e[N>>1];
        rep(i,0,n-1) if(rev[i]<i) swap(a[i],a[rev[i]]);
        e[0]=1;
        for(int i=1;i<n;i<<=1) {
            int t=qpow(f==1?3:(P+1)/3,(P-1)/i/2);
            for(int j=i-2;j>=0;j-=2) e[j+1]=1ll*t*(e[j]=e[j>>1])%P;
            for(int l=0;l<n;l+=i*2) {
                for(int j=l;j<l+i;++j) {
                    int t=1ll*a[j+i]*e[j-l]%P;
                    a[j+i]=a[j]-t,Mod2(a[j+i]);
                    a[j]+=t,Mod1(a[j]);
                }
            }
        }
        if(f==-1) {
            ll base=Inv[n];
            rep(i,0,n-1) a[i]=a[i]*base%P;
        }
    }
    void NTT(int n,int *a,int f){ 
        static int e[N>>1];
        rep(i,0,n-1) if(rev[i]<i) swap(a[i],a[rev[i]]);
        e[0]=1;
        for(int i=1;i<n;i<<=1) {
            int t=qpow(f==1?3:(P+1)/3,(P-1)/i/2);
            for(int j=i-2;j>=0;j-=2) e[j+1]=1ll*t*(e[j]=e[j>>1])%P;
            for(int l=0;l<n;l+=i*2) {
                for(int j=l;j<l+i;++j) {
                    int t=1ll*a[j+i]*e[j-l]%P;
                    a[j+i]=a[j]-t,Mod2(a[j+i]);
                    a[j]+=t,Mod1(a[j]);
                }
            }
        }
        if(f==-1) {
            ll base=Inv[n];
            rep(i,0,n-1) a[i]=a[i]*base%P;
        }
    }

#endif


    Poly operator * (Poly a,Poly b){
        int n=a.size()+b.size()-1;
        int R=Init(n);
        a.resize(R),b.resize(R);
        NTT(R,a,1),NTT(R,b,1);
        rep(i,0,R-1) a[i]=1ll*a[i]*b[i]%P;
        NTT(R,a,-1);
        a.resize(n);
        return a;
    }

    Poly operator + (Poly a,Poly b) { 
        int n=max(a.size(),b.size());
        a.resize(n),b.resize(n);
        rep(i,0,n-1) a[i]+=b[i],Mod1(a[i]);
        return a; 
    }
    Poly operator - (Poly a,Poly b) { 
        int n=max(a.size(),b.size());
        a.resize(n),b.resize(n);
        rep(i,0,n-1) a[i]-=b[i],Mod2(a[i]);
        return a; 
    }

    Poly Poly_Inv(Poly a) { // 多项式乘法逆，注意这里求出的是前a.size()项
        int n=a.size();
        if(n==1) return Poly{(int)qpow(a[0],P-2)};
        Poly b=a; b.resize((n+1)/2); b=Poly_Inv(b);
        int R=Init(n<<1);
        a.resize(R),b.resize(R);
        NTT(R,a,1),NTT(R,b,1);
        rep(i,0,R-1) a[i]=(2-1ll*a[i]*b[i]%P+P)*b[i]%P;
        NTT(R,a,-1);
        a.resize(n);
        return a;
    }

    Poly operator / (Poly a,Poly b){ // 多项式带余除法
        reverse(a.begin(),a.end()),reverse(b.begin(),b.end());
        int n=a.size(),m=b.size();
        a.resize(n-m+1),b.resize(n-m+1),b=Poly_Inv(b);
        a=a*b,a.resize(n-m+1);
        reverse(a.begin(),a.end());
        return a;
    }
    Poly operator % (Poly a,Poly b) { // 多项式取模
        int n=b.size()-1;
        if((int)a.size()<=n) return a;
        Poly t=a/b;
        if((int)t.size()>n) t.resize(n);
        t=t*b; t.resize(n); a.resize(n);
        return a-t;
    }

    int Quad(int a,int k=0) { // 二次剩余(不是原根法),用于求Sqrt
        if(a<=1) return a;
        ll x;
        while(1) {
            x=1ll*rand()*rand()%P;
            if(qpow((x*x-a+P)%P,(P-1)/2)!=1) break;
        }
        ll w=(x*x-a+P)%P;
        Pii res=mp(1,0),t=mp(x,1);
        auto Mul=[&](Pii a,Pii b){
            int x=(1ll*a.first*b.first+1ll*a.second*b.second%P*w)%P,y=(1ll*a.first*b.second+1ll*a.second*b.first)%P;
            return mp(x,y);
        };
        int d=(P+1)/2;
        while(d) {
            if(d&1) res=Mul(res,t);
            t=Mul(t,t);
            d>>=1;
        }
        ll r=(res.first%P+P)%P;
        if(k) r=min(r,(P-r)%P);
        return r;
    }
    Poly Sqrt(Poly a){ // 多项式开根号
        int n=a.size();
        if(n==1) return Poly{Quad(a[0],1)};
        Poly b=a; b.resize((n+1)/2),b=Sqrt(b),b.resize(n);
        Poly c=Poly_Inv(b);
        int R=Init(n*2);
        a.resize(R),c.resize(R);
        NTT(R,a,1),NTT(R,c,1);
        rep(i,0,R-1) a[i]=1ll*a[i]*c[i]%P;
        NTT(R,a,-1);
        a.resize(n);
        rep(i,0,n-1) a[i]=1ll*(P+1)/2*(a[i]+b[i])%P;
        return a;
    }

    Poly Deri(Poly a){ //求导
        rep(i,1,a.size()-1) a[i-1]=1ll*i*a[i]%P;
        a.pop_back();
        return a;
    }
    Poly IDeri(Poly a) { //原函数
        a.pb(0);
        drep(i,a.size()-1,1) a[i]=1ll*a[i-1]*Inv[i]%P;
        a[0]=0;
        return a;
    }

    Poly Ln(Poly a){ // 多项式求Ln
        int n=a.size();
        a=Poly_Inv(a)*Deri(a),a.resize(n-1);
        return IDeri(a);
    }
    Poly Exp(Poly a){ // 多项式Exp
        int n=a.size();
        if(n==1) return Poly{1};
        Poly b=a; b.resize((n+1)/2),b=Exp(b); b.resize(n);
        Poly c=Ln(b);
        rep(i,0,n-1) c[i]=a[i]-c[i],Mod2(c[i]);
        c[0]++,b=b*c;
        b.resize(n);
        return b;
    }

    void Exp_Solve(Poly &A,Poly &B,int l,int r){
        static int X[N],Y[N];
        if(l==r) {
            B[l]=1ll*B[l]*Inv[l]%P;
            return;
        }
        int mid=(l+r)>>1;
        Exp_Solve(A,B,l,mid);
        int R=Init(r-l+2);
        rep(i,0,R) X[i]=Y[i]=0;
        rep(i,l,mid) X[i-l]=B[i];
        rep(i,0,r-l-1) Y[i+1]=A[i];
        NTT(R,X,1),NTT(R,Y,1);
        rep(i,0,R-1) X[i]=1ll*X[i]*Y[i]%P;
        NTT(R,X,-1);
        rep(i,mid+1,r) B[i]+=X[i-l],Mod1(B[i]);
        Exp_Solve(A,B,mid+1,r);
    }
    Poly CDQ_Exp(Poly F){
        int n=F.size(); F=Deri(F);
        Poly A(n);
        A[0]=1;
        Exp_Solve(F,A,0,n-1);
        return A;
    }


    Poly Pow(Poly x,int k) { // 多项式k次幂
        x=Ln(x);
        rep(i,0,x.size()-1) x[i]=1ll*x[i]*k%P;
        return Exp(x);
    }

    Poly EvaluateTemp[N<<1];
    void EvaluateSolve1(Poly &a,int l,int r,int p=1){
        if(l==r) { EvaluateTemp[p]=Poly{P-a[l],1}; return; } 
        int mid=(l+r)>>1;
        EvaluateSolve1(a,l,mid,p<<1),EvaluateSolve1(a,mid+1,r,p<<1|1);
        EvaluateTemp[p]=EvaluateTemp[p<<1]*EvaluateTemp[p<<1|1];
    }
    void EvaluateSolve2(Poly &res,Poly F,int l,int r,int p=1){
        if(l==r){ res[l]=F[0]; return; }
        int mid=(l+r)>>1;
        EvaluateSolve2(res,F%EvaluateTemp[p<<1],l,mid,p<<1);
        EvaluateSolve2(res,F%EvaluateTemp[p<<1|1],mid+1,r,p<<1|1);
    }
    Poly Evaluate(Poly a,Poly b,int flag=1){ // 多项式多点求值
        Poly res(b.size());
        if(flag) EvaluateSolve1(b,0,b.size()-1);
        EvaluateSolve2(res,a,0,b.size()-1);
        return res;
    }
    Poly InterpolationSolve(Poly &T,int l,int r,int p=1){ 
        if(l==r) return Poly{T[l]};
        int mid=(l+r)>>1;
        return InterpolationSolve(T,l,mid,p<<1)*EvaluateTemp[p<<1|1]+InterpolationSolve(T,mid+1,r,p<<1|1)*EvaluateTemp[p<<1];
    }
    Poly Interpolation(Poly X,Poly Y){ // 多项式快速插值
        int n=X.size();
        EvaluateSolve1(X,0,n-1);
        Poly T=Evaluate(Deri(EvaluateTemp[1]),X,0);
        rep(i,0,n-1) T[i]=Y[i]*qpow(T[i])%P;
        return InterpolationSolve(T,0,n-1);
    }

    void FFPTrans(Poly &a,int f){ // FFP<->EGF
        int n=a.size();
        Poly b(n);
        if(f==1) rep(i,0,n-1) b[i]=FInv[i];
        else rep(i,0,n-1) b[i]=(i&1)?P-FInv[i]:FInv[i];
        a=a*b; a.resize(n);
    }
    Poly FFPMul(Poly a,Poly b){ // FFP卷积
        int n=a.size()+b.size()-1;
        a.resize(n),b.resize(n);
        FFPTrans(a,1),FFPTrans(b,1);
        rep(i,0,n-1) a[i]=1ll*a[i]*b[i]%P*Fac[i]%P;
        FFPTrans(a,-1);
        return a;
    }
    Poly PolyToFFP(Poly F){ // 多项式转FFP
        int n=F.size();
        Poly G(n);
        rep(i,0,n-1) G[i]=i;
        G=Evaluate(F,G);
        rep(i,0,n-1) F[i]=1ll*G[i]*FInv[i]%P;
        FFPTrans(F,-1);
        return F;
    }
    Poly FFPToPoly(Poly F){ // FFP转多项式
        FFPTrans(F,1);
        int n=F.size(); Poly X(n);
        rep(i,0,n-1) X[i]=i,F[i]=1ll*F[i]*Fac[i]%P;
        EvaluateSolve1(X,0,n-1);
        rep(i,0,n-1) {
            F[i]=1ll*F[i]*FInv[i]%P*FInv[n-i-1]%P;
            if((n-i-1)&1) F[i]=(P-F[i])%P;
        }
        return InterpolationSolve(F,0,n-1);
    }
}
using namespace Polynomial;

Poly Lag(int n,Poly X,Poly Y){
    Poly T(n+1),R(n+1),A(n+1);
    T[0]=1;
    rep(i,0,n) drep(j,i+1,0) T[j]=(1ll*T[j]*(P-X[i])+(j?T[j-1]:0))%P;
    rep(i,0,n) {
        ll t=1;
        rep(j,0,n) if(i!=j) t=t*(X[i]-X[j]+P)%P;
        t=qpow(t)*Y[i]%P,R[n+1]=T[n+1];
        drep(j,n,0) A[j]=(A[j]+t*R[j+1])%P,R[j]=(T[j]+1ll*R[j+1]*X[i]%P+P)%P;
    }
    return A;
}

int main(){
    int n=rd();
    Init_w();
    Poly F(n);
    rep(i,0,n-1) F[i]=rd();
    Show(CDQ_Exp(F));
}

\[\ \]

posted @ 2020-04-29 14:01 chasedeath 阅读(1240) 评论(1) 编辑收藏举报

刷新页面返回顶部

orangejuice

多项式运算 (求逆/ln/exp等)

多项式运算 (求逆/ln/exp等)

1.多项式求乘法逆

Part1 O(\(n^2\))

Part2 O(\(n\log^2n\))

Part3 \(O(n\log n)\)

2.多项式开根号

3.多项式求\(\ln\)

4.多项式求exp

CDQ分治+\(\text{NTT}\)

牛顿迭代

5.多项式\(k\)次幂

6.多项式带余除法

公告