题解 P4726 【【模板】多项式指数函数】

· · 题解

题意

给定一个多项式A(x)，求一个多项式B(x)，满足

B(x)\equiv e^{A(x)}(\bmod{\:x^n})

题解

默认大家都会多项式求逆和多项式求ln，不会的请左转这里和这里

预知识1：牛顿迭代

先看一道小问题：手算，如何对一个10^9大小的数a开方？

如果你二分，那么预计要二分30次左右，每一次要算一个5位数的平方。

于是牛顿法产生了。

我们本质上是要求f(x)=x^2-a精确到整数的零点。

假设我们已经求得了一个近似值x_0，那么我们只需要过(x_0,f(x_0))这个点作这个函数图像的切线，取切线与x轴的交点作为新的x_0。

丢一张图：

可以看到，这样每次求切线，近似值会迅速逼近精确值，这样可以大大减少工作量。实践证明对201723039（我的考号，这是考数学的时候发生的一个真实的故事）开方初值取到40000只需要迭代四次就可以精确到整数。

假设我们要求一个函数f(x)的零点，初始近似值是x_0。

则切线方程为

y=f'(x_0)(x-x_0)+f(x_0)

（不会的好好看看高中数学选修2-2）

令y=0，得到x=x_0-\dfrac{f(x_0)}{f'(x_0)}。

放到多项式上也是同理的。

假设我们要求一个函数G(x)，满足F(G(x))\equiv 0。

利用上面的式子，我们可以每一次令

G(x)=G_0(x)-\dfrac{F(G_0(x))}{F'(G_0(x))}

这样它可以迅速逼近真实值。

本质上，迭代一次精度就可以翻倍。如果F(G_0(x))\equiv 0(\bmod{\:x^{n\over 2}})，那么F(G(x))\equiv 0(\bmod{\:x^{n}})，证明我也不会，只能感性理解。

于是我们可以每次将n除以2，递归进去做，然后用上面那个牛顿迭代的式子去做。

用牛顿迭代推多项式exp

B(x)\equiv e^{A(x)}(\bmod{\:x^n}) \ln B(x)-A(x)\equiv 0(\bmod{\:x^n})

所以我们只需要使得F(G(x))=\ln G(x)-A(x)\equiv 0。

\Big(F\big(G_0(x)\big)\Big)'=\dfrac{G_0'(x)}{G_0(x)}

（复合函数求导法则不能忘啊）

带入上面那个式子得到：

G(x)=G_0(x)-\dfrac{\ln G_0(x)-A(x)}{\dfrac{G_0'(x)}{G_0(x)}}=\dfrac{G_0(x)(1-\ln G_0(x)+A(x))}{G_0'(x)}

所以，我们只需要每次迭代做一遍多项式求逆，再做一遍多项式求ln，然后再做一遍多项式乘法，就可以得到答案。

时间复杂度为

T(n)=T\left(n\over 2\right)+O(n\log n)=O(n\log n)

没错，怎么牛顿迭代都是O(n\log n)！

下面是代码：

#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
#define ll long long
using namespace std;
const ll MOD=998244353,inv2=499122177;
int limit,l,r[400005];
ll ni[400005];
void calc_ni(int n)
{
    ni[1]=1;
    for(int i=2;i<=n;i++)
      ni[i]=MOD-(MOD/i)*ni[MOD%i]%MOD;
}
ll quick_pow(ll x,ll a)
{
    ll ans=1;
    while(a)
    {
        if(a&1)ans=ans*x%MOD;
        x=x*x%MOD;
        a>>=1;
    }
    return ans;
}
void NTT(ll*A,int type)
{
    for(int i=0;i<limit;i++)
      if(i<r[i])swap(A[i],A[r[i]]);
    for(int mid=1;mid<limit;mid<<=1)
    {
        ll Wn=type==1?quick_pow(3,(MOD-1)/(mid<<1)):quick_pow(3,MOD-1-(MOD-1)/(mid<<1));
        for(int R=mid<<1,j=0;j<limit;j+=R)
        {
            ll w=1;
            for(int k=0;k<mid;k++,w=w*Wn%MOD)
            {
                ll x=A[j+k],y=w*A[j+mid+k]%MOD;
                A[j+k]=(x+y)%MOD;
                A[j+mid+k]=(x-y+MOD)%MOD;
            }
        }
    }
    if(type==-1)
    {
        ll inv=quick_pow(limit,MOD-2);
        for(int i=0;i<limit;i++)A[i]=A[i]*inv%MOD;
    }
}
ll A[400005];
void Inv(ll*a,ll*b,int n)
{
    if(n==1)
    {
        b[0]=quick_pow(a[0],MOD-2);
        return;
    }
    Inv(a,b,(n+1)>>1);
    limit=1,l=0;
    while(limit<(n<<1))limit<<=1,l++;
    for(int i=0;i<limit;i++)
      r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
    for(int i=0;i<limit;i++)A[i]=i<n?a[i]:0;
    NTT(A,1);
    NTT(b,1);
    for(int i=0;i<limit;i++)
      b[i]=b[i]*(2+MOD-A[i]*b[i]%MOD)%MOD;
    NTT(b,-1);
    for(int i=n;i<limit;i++)b[i]=0;
}
ll a2[400005];
void Ln(ll*a,ll*b,int n)
{
    for(int i=0;i<(n<<2);i++)b[i]=0;
    Inv(a,b,n);
    limit=1,l=0;
    while(limit<(n<<1))limit<<=1,l++;
    for(int i=0;i<limit;i++)
      r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
    for(int i=0;i<n-1;i++)a2[i]=a[i+1]*(i+1)%MOD;
    for(int i=n-1;i<limit;i++)a2[i]=0;
    NTT(a2,1);
    NTT(b,1);
    for(int i=0;i<limit;i++)
      b[i]=b[i]*a2[i]%MOD;
    NTT(b,-1);
    for(int i=n-1;i>0;i--)
      b[i]=b[i-1]*ni[i]%MOD;
    for(int i=n;i<limit;i++)b[i]=0;
    b[0]=0;
}
ll lnb[400005];
void Exp(ll*a,ll*b,int n)
{
    if(n==1)
    {
        b[0]=1;
        return;
    }
    Exp(a,b,(n+1)>>1);
    Ln(b,lnb,n);
    limit=1,l=0;
    while(limit<(n<<1))limit<<=1,l++;
    for(int i=0;i<limit;i++)
      r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
    for(int i=0;i<n;i++)lnb[i]=a[i]>=lnb[i]?a[i]-lnb[i]:a[i]-lnb[i]+MOD;
    for(int i=n;i<limit;i++)lnb[i]=b[i]=0;
    lnb[0]++;
    NTT(lnb,1);
    NTT(b,1);
    for(int i=0;i<limit;i++)b[i]=b[i]*lnb[i]%MOD;
    NTT(b,-1);
    for(int i=n;i<limit;i++)b[i]=0;
}
int n;
ll f[400005],g[400005];
int main()
{
    scanf("%d",&n);
    for(int i=0;i<n;i++)scanf("%lld",&f[i]);
    calc_ni(n);
    Exp(f,g,n);
    for(int i=0;i<n;i++)printf("%lld ",g[i]);
    printf("\n");
    return 0;
}

常数巨大，慎用

upd: BUG 修复