多项式的ln、exp、快速幂和开根学习小记

• 不妨又学习了一下多项式的求ln、exp、快速幂和开根操作。

• 这些操作比之前的求逆更上了一层台阶，应用同样很广。

• 多项式求逆等知识在我的博客里有讲：多项式的求逆、取模和多点求值学习小记

多项式ln

• 给出一个次数界为 $n$ 的多项式 $F(x)$ ，要求多项式 $G(x)$ 满足：$$G(x)\equiv lnF(x)(mod{x}^n)$$

• 我们对其两边求导（右边相当于是复合函数求导）：$$G^{\prime }(x)=\frac{F^{\prime }(x)}{F(x)}$$

• 再两边积分，即：$$G(x)=\int \frac{F^{\prime }(x)}{F(x)}dx$$

• 于是应用一下多项式求逆即可在 $O(nlogn)$ 的时间内求出一个多项式的 ln 了。

• 其中求导和积分都是 $O(n)$ ，具体来说就是：$$(x^a{)}^{\prime }=a{x}^{a-1}$$$$\int x^{a}dx=\frac{1}{a+1}{x}^{a+1}$$

• 逐项操作即可求导、积分。

多项式exp

• 给出一个次数界为 $n$ 的多项式 $F(x)$ ，要求多项式 $G(x)$ 满足：$$G(x)\equiv e^{F(x)}(mod{x}^n)$$

• 这个需要用到 牛顿迭代法——一种在实数域和复数域上近似求解方程的方法。

• 牛顿迭代法的一般式如下（要求的根为 $x$）：$$x^{\prime }=x-\frac{f(x)}{f^{\prime }(x)}$$

• 于是我们先将定义式移项：$$G(x)-e^{F(x)}=0$$

• 两边取对数：$$lnG(x)-F(x)=0$$

• 把 $G(x)$ 当做要求的方程的根，则根据牛顿迭代法可得：$$G_{t+1}(x)=G_t(x)-\frac{ln{G}_t(x)-F(x)}{\frac{1}{G_t(x)}}=G_t(x)(1+F(x)-ln{G}_t(x))$$

• 这样每次迭代都使精度（次数界）翻倍，递归处理即可，每次递归过程中要进行一次多项式求ln 。

• 这样的时间复杂度是 $T(n)=T(\frac{n}{2})+O(nlogn)=O(nlogn)$ 的，不过常数较大。

多项式快速幂

• 给出一个次数界为 $n$ 的多项式 $F(x)$ ，要求多项式 $G(x)$ 满足：$$G(x)\equiv F(x{)}^k(mod{x}^n)$$

• 如果像普通一样每次直接 NTT 的话，复杂度会很不理想，特别是当 $k$ 大到一定程度的时候。

• 我们可以两边取对数，则有：$$lnG(x)=lnF(x{)}^k=klnF(x)$$

• 右边多项式系数直接乘 $k$ ，再两边exp还原即可。

• 这样做的话 $k$ 的值可以很大，我们只需要知道它模意义下的值即可。

• 这样多项式快速幂需要用到求ln、exp，但复杂度仍是 $O(nlogn)$ 的。

• 附模板题和我的代码： 洛谷 P5245 【模板】多项式快速幂

• 值得注意的是这里的 $k$ 达到了 $10^{10^5}$ 数量级，我们用高精度除单精度求其模意义下的值即可。

Code

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cctype>
using namespace std;
typedef long long LL;
const int N=1e5+5,L=N<<2,G=3,mo=998244353;
int a[L],b[L],c[L];
int f[L],rev[L],wn[L],inv[L];
char s[N];
inline int read()
{
    int X=0,w=0; char ch=0;
    while(!isdigit(ch)) w|=ch=='-',ch=getchar();
    while(isdigit(ch)) X=(X<<3)+(X<<1)+(ch^48),ch=getchar();
    return w?-X:X;
}
inline int ksm(int x,int y)
{
    int s=1;
    while(y)
    {
        if(y&1) s=(LL)s*x%mo;
        x=(LL)x*x%mo;
        y>>=1;
    }
    return s;
}
inline void NTT(int *y,int len,int ff)
{
    for(int i=0;i<len;i++)
        if(i<rev[i]) swap(y[i],y[rev[i]]);
    for(int h=2,d=len>>1;h<=len;h<<=1,d>>=1)
        for(int i=0,k=h>>1;i<len;i+=h)
            for(int j=0,cnt=0;j<k;j++,cnt+=d)
            {
                int u=y[i+j],t=(LL)y[i+j+k]*wn[cnt]%mo;
                y[i+j]=u+t>=mo?u+t-mo:u+t;
                y[i+j+k]=u-t<0?u-t+mo:u-t;
            }
    if(ff==-1)
    {
        reverse(y+1,y+len);
        int invl=inv[len];
        for(int i=0;i<len;i++) y[i]=(LL)y[i]*invl%mo;
    }
}
void getinv(int len,int num)
{
    if(len==1)
    {
        b[0]=ksm(c[0],mo-2);
        return;
    }
    getinv(len>>1,num-1);
    for(int i=0;i<len>>1;i++) f[i]=c[i];
    for(int i=len>>1;i<len;i++) f[i]=0;
    for(int i=0;i<len;i++) rev[i]=rev[i>>1]>>1|(i&1)<<num-1;
    int w0=ksm(G,(mo-1)/len);
    for(int i=wn[0]=1;i<=len;i++) wn[i]=(LL)wn[i-1]*w0%mo;
    NTT(f,len,1),NTT(b,len,1);
    for(int i=0;i<len;i++) b[i]=(2-(LL)b[i]*f[i]%mo+mo)*b[i]%mo;
    NTT(b,len,-1);
    for(int i=len>>1;i<len;i++) b[i]=0;
}
inline void getln(int len,int num)
{
    getinv(len,num);
    for(int i=0;i<len-1;i++) f[i]=(LL)c[i+1]*(i+1)%mo;
    f[len-1]=0;
    for(int i=len>>1;i<len;i++) f[i]=b[i]=0;
    NTT(f,len,1),NTT(b,len,1);
    for(int i=0;i<len;i++) f[i]=(LL)f[i]*b[i]%mo;
    NTT(f,len,-1);
    for(int i=len-1;i>0;i--) f[i]=(LL)f[i-1]*inv[i]%mo;
    f[0]=0;
}
void getexp(int len,int num)
{
    if(len==1)
    {
        c[0]=1;
        return;
    }
    getexp(len>>1,num-1);
    for(int i=0;i<len;i++) f[i]=b[i]=0;
    getln(len,num);
    f[0]=(a[0]+1-f[0]+mo)%mo;
    for(int i=1;i<len>>1;i++) f[i]=(a[i]-f[i]+mo)%mo;
    for(int i=len>>1;i<len;i++) f[i]=0;
    NTT(c,len,1),NTT(f,len,1);
    for(int i=0;i<len;i++) c[i]=(LL)c[i]*f[i]%mo;
    NTT(c,len,-1);
    for(int i=len>>1;i<len;i++) c[i]=0;
}
int main()
{
    int n=read(),k=0;
    scanf("%s",s+1);
    int m=strlen(s+1);
    for(int i=1;i<=m;i++) k=((LL)k*10+s[i]-'0')%mo;
    for(int i=0;i<n;i++) c[i]=read();
    inv[0]=inv[1]=1;
    for(int i=2;i<L;i++) inv[i]=(LL)(mo-mo/i)*inv[mo%i]%mo;
    int len=1,num=0;
    while(len<n<<1) len<<=1,num++;
    getln(len,num);
    for(int i=0;i<n;i++) a[i]=(LL)f[i]*k%mo;
    for(int i=0;i<len;i++) f[i]=b[i]=c[i]=0;
    getexp(len,num);
    for(int i=0;i<n;i++) printf("%d ",c[i]);
    return 0;
}

多项式开根

• 给出一个次数界为 $n$ 多项式 $A(x)$ ，要求一个次数界为 $n$ 的多项式 $B(x)$ 满足：$$B^2(x)\equiv A(x)(mod{x}^n)(\ast )$$

• 这个在做生成函数题时用得很多，放在这里是因为这可以用 ln、exp 来处理（需要保证 $A_0=1$）。

• 两边取对数，得：$$lnB(x)=\frac{1}{2}lnA(x)$$

• 再两边 exp 还原即可，本质与多项式快速幂是一样的，时间复杂度 $O(nlogn)$ 。

• 当然开根也有其独特的解法，跟求逆一样是倍增来求的。

• 我们已求出 $G(x)$ 满足：$$G^2(x)\equiv A(x)(mod{x}^{\lfloor \frac{n}{2}\rfloor })$$

• 用 (*) 式减上式，得：$$B^2(x)-G^2(x)\equiv 0(mod{x}^{\lfloor \frac{n}{2}\rfloor })$$

• 平方差公式：$$(B(x)+G(x))(B(x)-G(x))\equiv 0(mod{x}^{\lfloor \frac{n}{2}\rfloor })$$

• 这里的根有两个，这里考虑其中一个：$$B(x)-G(x)\equiv 0(mod{x}^{\lfloor \frac{n}{2}\rfloor })$$

• 两边平方（次数界扩大了两倍）：$$(B(x)-G(x){)}^2\equiv 0(mod{x}^n)$$

• 展开来：$$B^2(x)-2B(x)G(x)+G^2(x)\equiv 0(mod{x}^n)$$

• 消去 $B^2(x)$ ，即：$$A(x)-2B(x)G(x)+G^2(x)\equiv 0(mod{x}^n)$$

• 移项即得：$$B(x)\equiv \frac{A(x)+G^2(x)}{2G(x)}(mod{x}^n)$$

• 或者这种形式：$$B(x)\equiv \frac{1}{2}(\frac{A(x)}{G(x)}+G(x))(mod{x}^n)$$

• 于是我们多项式求逆并计算即可，时间复杂度 $O(nlogn)$ 。

• 注意递归底层是 $B_0=\sqrt{A_0}$ ，若是 0、1 就最好了，否则我们得算一下二次剩余（当然必须得有，不然没法开根）！

• 如果题目使得这个值是固定的，那么就可以直接枚举算二次剩余是啥。

• 如果是会变的，也许我们要用用什么 $Cipolla$ 算法……这里就不做讨论了。。

• 有兴趣的同学可以看看 a_crazy_czy大神 的 二次剩余Cipolla算法学习小记 。

• 模板题：洛谷 P5205 【模板】多项式开根

Code

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cctype>
using namespace std;
typedef long long LL;
const int N=1e5+5;
int G=3,mo=998244353,inv2=mo+1>>1;
int a[N<<2],b[N<<2],c[N<<2];
int f[N<<2],rev[N<<2],wn[N<<2];
inline int read()
{
    int X=0,w=0; char ch=0;
    while(!isdigit(ch)) w|=ch=='-',ch=getchar();
    while(isdigit(ch)) X=(X<<3)+(X<<1)+(ch^48),ch=getchar();
    return w?-X:X;
}
void write(int x)
{
    if(x>9) write(x/10);
    putchar(x%10+'0');
}
inline int ksm(int x,int y)
{
    int s=1;
    while(y)
    {
        if(y&1) s=(LL)s*x%mo;
        x=(LL)x*x%mo;
        y>>=1;
    }
    return s;
}
inline void NTT(int *y,int len,int ff)
{
    for(int i=0;i<len;i++)
        if(i<rev[i]) swap(y[i],y[rev[i]]);
    for(int h=2,d=len>>1;h<=len;h<<=1,d>>=1)
        for(int i=0,k=h>>1;i<len;i+=h)
            for(int j=0,cnt=0;j<k;j++,cnt+=d)
            {
                int u=y[i+j],t=(LL)wn[cnt]*y[i+j+k]%mo;
                y[i+j]=u+t>=mo?u+t-mo:u+t;
                y[i+j+k]=u-t<0?u-t+mo:u-t;
            }
    if(ff==-1)
    {
        for(int i=len>>1;i;i--) swap(y[i],y[len-i]);
        int inv=ksm(len,mo-2);
        for(int i=0;i<len;i++) y[i]=(LL)y[i]*inv%mo;
    }
}
void solve1(int len,int num)
{
    if(len==1)
    {
        c[0]=ksm(b[0],mo-2);
        return;
    }
    solve1(len>>1,num-1);
    for(int i=0;i<len;i++) rev[i]=rev[i>>1]>>1|(i&1)<<num-1;
    int w0=ksm(G,(mo-1)/len);
    for(int i=wn[0]=1;i<=len;i++) wn[i]=(LL)wn[i-1]*w0%mo;
    for(int i=0;i<len>>1;i++) f[i]=b[i];
    for(int i=len>>1;i<len;i++) f[i]=0;
    NTT(f,len,1),NTT(c,len,1);
    for(int i=0;i<len;i++) f[i]=(2-(LL)f[i]*c[i]%mo+mo)*c[i]%mo;
    NTT(f,len,-1);
    for(int i=0;i<len>>1;i++) c[i]=f[i];
    for(int i=len>>1;i<len;i++) c[i]=0;
}
void solve(int len,int num)
{
    if(len==1)
    {
        b[0]=1;
        return;
    }
    solve(len>>1,num-1);
    solve1(len,num);
    for(int i=0;i<len>>1;i++) f[i]=a[i];
    for(int i=len>>1;i<len;i++) f[i]=0;
    NTT(f,len,1),NTT(c,len,1);
    for(int i=0;i<len;i++) f[i]=(LL)f[i]*c[i]%mo;
    NTT(f,len,-1);
    for(int i=0;i<len>>1;i++) b[i]=(LL)(f[i]+b[i])*inv2%mo;
    for(int i=len>>1;i<len;i++) b[i]=0;
    for(int i=0;i<len;i++) c[i]=0;
}
int main()
{
    int n=read();
    for(int i=0;i<n;i++) a[i]=read();
    int mx=1,ll=0;
    while(mx<n+n) mx<<=1,ll++;
    solve(mx,ll);
    for(int i=0;i<n;i++) write(b[i]),putchar(' ');
    return 0;
}