多项式ln、牛顿迭代学习笔记

这篇博客写得很好,本文大部分参考此博客。

一些证明的细节在本文中不会提及。

前置知识:多项式求逆,多项式求导、积分 , 泰勒展开。

多项式ln

给定\(G(x)\),求

\[F(x)=ln\ G(x)\]

两边求导,得

\[F'(x)=\frac{G'(x)}{G(x)}\]

\[F(x)=\int \frac{G'(x)}{G(x)}\]

求逆,计算即可。

牛顿迭代

以下我们把一个多项式当做一个数来理解。

\(\delta(F(x))\)是一个关于多项式\(F(x)\)的一个函数,现在我们要求一个\(F(x)\),满足

\[\delta(F(x))=0\]

假设我们已经求出了\(F(x)\)的前\(\frac{n}{2}\)项,记为\(F_0(x)\),现在我们要求它的前\(n\)

\(\delta(F(x))\)进行泰勒展开,有

\[\delta(F(x))=\delta(F_0(x))+\frac{\delta '(F_0(x))}{1!}(F(x)-F_0(x))+\frac{\delta ''(F_0(x))}{2!}(F(x)-F_0(x))^2+...\]

其中\(\delta '(F(x))\)\(\delta (F(x))\)关于\(F(x)\)的导数

因为\(F_0(x)\)\(F(x)\)\(\frac{n}{2}\)项相同,所以\((F(x)-F_0(x))^2\equiv 0 \ \ (mod \ x^n)\),所以上式从第三项开始的值全为0。

也就是

\[\delta(F(x))=\delta(F_0(x))+\delta '(F_0(x))(F(x)-F_0(x))=0\]

整理得

\[F(x)=F_0(x)-\frac{\delta(F_0(x))}{\delta '(F_0(x))}\]

如果没搞懂的话可以参照下面的例子来理解。

牛顿迭代——多项式求逆

给定\(G(x)\),求\(F(x)\)满足

\[F(x)G(x)\equiv1\ (mod\ x^n)\]

就是要使得\(F(x)G(x)-1=0\)

所以设 \(\delta(F(x))=F(x)G(x)- 1\)

就有

\[F(x)=F_0(x)-\frac{\delta(F_0(x))}{\delta '(F_0(x))}\]

\[=F_0(x)-\frac{F_0(x)G(x)-1}{G(x)}\]

\[=F_0(x)-F_0(F_0(x)G(x)-1)\]

\[=2F_0(x)-F_0^2(x)G(x)\ \ (mod\ x^n)\]

牛顿迭代——多项式开根

即求\(F^2(x)-G(x)=0\)

\[F(x)=F_0(x)-\frac{F_0^2(x)-G(x)}{2F_0(x)}\]

\[=\frac{F_0^2(x)+G(x)}{2F_0(x)}\]

多项式Exp

\[F(x)=e^ {G(x)}\]

\[ln\ F(x)=G(x)\]

\[ln\ F(x)-G(x)=0\]

\[ F(x)=F_0(x)-\frac{ln\ F_0(x)-G(x)}{\frac{1}{F_0(x)}}\]

\[=F_0(x)-F_0(x)(ln\ F_0(x)-G(x))\]

\[=F_0(x)(1-ln\ F_0(x)+G(x))\]

最后奉上我丑陋的代码模板:

#include<bits/stdc++.h>
using namespace std;
#define N 600007
#define M 2100007
#define ll long long
const ll mod=998244353;
const int lim=4e5;
ll inv[N];
int rev[M],len;
ll p2(ll x){return x*x%mod;}
ll pw(ll x,ll p)
{
    return p?p2(pw(x,p/2))*(p&1?x:1)%mod:1;
}
void getrev()
{
    for(int i=0;i<len;i++)
        rev[i]=rev[i>>1]>>1|(i&1?len>>1:0);
}
void getlen(int n)
{
    for(len=1;len<=n;len<<=1);
    getrev();
}
void NTT(ll *a,int op)
{
    for(int i=0;i<len;i++)if(i<rev[i])swap(a[i],a[rev[i]]);
    for(int i=1;i<len;i<<=1)
    {
        ll nw=pw(3,(mod-1)/(i<<1));
        for(int j=0;j<len;j+=i<<1)
        {
            ll w=1;
            for(int k=j;k<j+i;k++)
            {
                ll x=a[k],y=a[k+i]*w%mod;
                a[k]=(x+y)%mod,a[k+i]=(x-y+mod)%mod;
                w=w*nw%mod;
            }
        }
    }
    if(op<0)
    {
        reverse(a+1,a+len);
        ll Inv=pw(len,mod-2);
        for(int i=0;i<len;i++)a[i]=a[i]*Inv%mod;
    }
}
ll inv_c[N];
void getinv(int p,ll *a,ll *b)
{
    if(p==1)return a[0]=pw(b[0],mod-2),(void)1;
    getinv((p+1)/2,a,b);
    getlen(p<<1);
    ll *c=inv_c;
    for(int i=0;i<p;i++)c[i]=b[i];
    for(int i=p;i<len;i++)c[i]=0;
    NTT(a,1),NTT(c,1);
    for(int i=0;i<len;i++)a[i]=a[i]*(2-a[i]*c[i]%mod+mod)%mod;
    NTT(a,-1);
    for(int i=p;i<len;i++)a[i]=0;
}
ll mul_c[M],mul_d[M];
void mul(ll *t,ll *a,ll *b)
{
    ll *c=mul_c,*d=mul_d;
    for(int i=0;i<len;i++)c[i]=a[i],d[i]=b[i];
    NTT(c,1),NTT(d,1);
    for(int i=0;i<len;i++)c[i]=c[i]*d[i]%mod;
    NTT(c,-1);
    for(int i=0;i<len;i++)t[i]=c[i];
}
void devir(ll *a,int n)
{
    for(int i=1;i<=n;i++)a[i-1]=a[i]*i%mod;
    a[n]=0;
}
void inter(ll *a,int n)
{
    for(int i=n;i>=0;i--)a[i+1]=a[i]*inv[i+1]%mod;
    a[0]=0;
}
ll ln_c[N],ln_d[N];
void getln(int n,ll *a,ll *b)
{
    ll *c=ln_c,*d=ln_d;
    getlen(2*n);
    for(int i=0;i<len;i++)c[i]=b[i],d[i]=0;
    getinv(n,d,c),devir(c,n);
    getlen(2*n);
    mul(a,c,d);
    inter(a,n);
    for(int i=n;i<len;i++)a[i]=0;
}
ll exp_c[N];
void getexp(int p,ll *a,ll *b)
{
    if(p==1)return a[0]=1,(void)1;
    getexp((p+1)/2,a,b);
    ll *c=exp_c;
    getlen(2*p);
    for(int i=0;i<len;i++)c[i]=0;
    getln(p,c,a);
    for(int i=0;i<p;i++)c[i]=(b[i]-c[i]+mod)%mod;
    c[0]=(c[0]+1)%mod;
    getlen(2*p);
    mul(a,a,c);
    for(int i=p;i<len;i++)a[i]=0;
}
ll sqrt_c[N],sqrt_d[N];
void getsqrt(int p,ll *a,ll *b)
{
    if(p==1)return a[0]=1,(void)1;
    getsqrt((p+1)/2,a,b);
    ll *c=sqrt_c,*d=sqrt_d;
    getlen(2*p);
    for(int i=0;i<len;i++)c[i]=d[i]=0;
    for(int i=0;i<p;i++)d[i]=b[i];
    getinv(p,c,a);
    getlen(2*p),mul(c,c,d);
    for(int i=0;i<p;i++)a[i]=(a[i]+c[i])*inv[2]%mod;
}
ll a[M],b[M],c[M];
int n;
void work_mul()
{
    int n,m;
    scanf("%d%d",&n,&m);
    for(int i=0;i<=n;i++)scanf("%lld",&a[i]);
    for(int i=0;i<=m;i++)scanf("%lld",&b[i]);
    getlen(n+m);
    mul(c,a,b);
    for(int i=0;i<=n+m;i++)printf("%lld ",c[i]);
}
void read()
{
    scanf("%d",&n);
    for(int i=0;i<n;i++)scanf("%lld",&a[i]);
}
void write(ll *a)
{
    for(int i=0;i<n;i++)printf("%lld ",a[i]);
}
void work_inv()
{
    read();
    getinv(n,b,a);
    write(b);
}
void work_sqrt()
{
    read();
    getsqrt(n,b,a);
    write(b);
}
void work_ln()
{
    read();
    getln(n,b,a);
    write(b);
}
void work_exp()
{
    read();
    getexp(n,b,a);
    write(b);
}
void Init()
{
    inv[1]=1;
    for(int i=2;i<=lim;i++)
        inv[i]=mod-mod/i*inv[mod%i]%mod;
}
int main()
{
    //freopen("data.in","r",stdin);
    Init();
    work_exp();
    return 0;
}

猜你喜欢

转载自www.cnblogs.com/lishuyu2003/p/12123705.html