多项式牛顿迭代

tNuxvs5t · 2023-06-10 11:03:02 · 个人记录

以下仅为作者的个人理解，有错轻喷。

泰勒级数

定义 多项式拟合，作用于函数 f(x)，有 f(x)=F(x)=\sum_{i=0}^\infty a_ix^i

考虑在点 x_0 泰勒展开的意义：在 x_0 处 0\to n 阶导数相同。

有如下性质：

代入（3）

化简得

因此证明了一个满足 g(F_0(x))=0\pmod {x^{n/2}} 的多项式通过一次迭代可以求出 g(F(x))=0\pmod {x^n}

并给出相关公式。

应用

常见的迭代方式主要有递归版和幂次版，由于常数关系，选择使用后者。

具体表现为：先求出 \pmod x 的解，然后层层迭代直到求出 \pmod {x^{lim}} 的解。

对多项式 exp：

求 B(x)，使 B(x)\equiv e^{A(x)}\pmod {x^n} 。

取对数

移项

根据多项式牛顿迭代法则

使用 \ln 和 卷积 模板，迭代即可。

代码：

#include <bits/stdc++.h>
using namespace std;
#define int long long
#define rep(i,a,b) for(int i=a;i<=b;i++)
const int N=(1<<18)+1,mod=998244353,yg=3;
int rot[N];
int qpow(int x,int y){
    int ans=1;
    while(y){
        if(y&1)ans=ans*x%mod;
        y>>=1;x=x*x%mod;
    }return ans;
}
const int nyg=qpow(yg,mod-2);
struct Poly{
    int *a;
    int& operator[](int pos){return a[pos];}
    int lg=0;
    bool operator==(Poly& b){
        if(lg!=b.lg)return false;
        rep(i,0,lg){
            if(a[i]!=b[i])return false;
        }
        return true;
    }
    Poly(int SZ){a=new int[SZ+2];rep(i,0,SZ+1)a[i]=0;}
    void ntt(int lim,int mode=1){
        rep(i,1,lim-1){
            rot[i]=(rot[i>>1]>>1)|((i&1)*(lim>>1));
        }
        rep(i,0,lim-1){
            if(i<rot[i])swap(a[i],a[rot[i]]);
        }
        for(int mid=1;mid<lim;mid<<=1){
            int w0=qpow((mode==1)?yg:nyg,(mod-1)/(mid<<1));
            for(int i=0;i<lim;i+=(mid<<1)){
                int wn=1;
                rep(j,0,mid-1){
                    int X=a[i+j],Y=a[i+j+mid]*wn%mod;
                    a[i+j]=(X+Y)%mod;
                    a[i+j+mid]=(X-Y+mod)%mod;
                    wn=wn*w0%mod;
                }
            }
        }
        if(mode!=1){
            int ny=qpow(lim,mod-2);
            rep(i,0,lim-1)a[i]=a[i]*ny%mod;
        }
    }
    static Poly mul(Poly a,Poly b){
        int yz=1;
        while(yz<a.lg+b.lg)yz<<=1;
        Poly _a(yz),_b(yz),c(yz);
        rep(i,0,a.lg)_a[i]=a[i];
        rep(i,0,b.lg)_b[i]=b[i];
        _a.ntt(yz,1);_b.ntt(yz,1);
        rep(i,0,yz-1){
            c[i]=_a[i]*_b[i]%mod;
        }
        c.ntt(yz,-1);
        c.lg=a.lg+b.lg;
        return c;
    }
    Poly inv(){
        int lim=1;
        while(lim<lg)lim<<=1;
        Poly u(lim);
        u[0]=qpow(a[0],mod-2);
        int tar=1;
        while(tar<lim){
            u.lg=tar*2;
            Poly w=mul(u,u);
            Poly tp(tar*2);tp.lg=tar*2;
            rep(i,0,tar*2)tp[i]=a[i];
            w=mul(w,tp);
            rep(i,0,tar*2)u[i]=((2*u[i]-w[i])%mod+mod)%mod;
            tar*=2;
        }
        u.lg=lg;
        return u;
    }
    Poly calc(){
        Poly u(lg+1);
        rep(i,1,lg+1){
            u[i]=a[i-1]*qpow(i,mod-2)%mod;
        }
        u.lg=lg+1;
        return u;
    }
    Poly deri(){
        Poly u(lg-1);
        rep(i,0,lg-1){
            u[i]=a[i+1]*(i+1)%mod;
        }
        u.lg=lg-1;
        return u;
    }
    Poly ln(){
        Poly u=mul(deri(),inv());
        u.lg=lg;
        return u.calc();
    }
    Poly exp(){
        Poly u(1);
        u[0]=1;
        u.lg=1;
        while(u.lg<2*lg){
            Poly u0=u.ln();
            rep(i,0,u0.lg){
                u0[i]=(-u0[i]+(i==0)+a[i])%mod+mod;
                u0[i]%=mod;
            }
            u=mul(u0,u);
        }
        u.lg=lg;
        return u;
    }
};
int n;
Poly g(N);
signed main(){
    ios::sync_with_stdio(0);
    cin>>n;
    rep(i,0,n-1)cin>>g[i];
    g.lg=n;
    Poly f=g.exp();
    rep(i,0,n-1)cout<<f[i]<<' ';
}