多项式总结

2018-05-11 18:37:51


此处多项式算法,主要针对基于快速傅里叶变换 $(FFT)$ 或者快速数论变换 $(NTT)$ 的多项式演算,主要用来解决各种母函数、排列组合、字符串题目。此处仅介绍相关底层算法。文章可能不定期更新。

[TOC]

1.多项式乘法

点我

一切算法的基础。直接做为 $O(n^2)$ ,利用 $FFT$ 或者 $NTT$ 来实现可以做到 $O(n\log n)$。

配合预处理单位根可以提升精度和速度。

代码:


inline void FFT(comp*F,int tl)
{
    rep(i,0,Len)if(i<rev[i])swap(F[rev[i]],F[i]);
    for(register int i=2,ii=1,t=1;i<=Len;i<<=1,ii<<=1,++t)
        for(register int j=0;j<Len;j+=i)rep(k,0,ii)
        {
            comp x=F[j+k+ii]*g[t][k];
            F[j+k+ii]=F[j+k]-x;
            F[j+k]=F[j+k]+x;
        }
    if(tl==-1)
    {
        reverse(F+1,F+Len);
        rep(i,0,Len)F[i]=F[i]/Len;
    }
}

inline void NTT(int*F,int typ)
{
    Rep(i,1,Len-1)if(i<rev[i])swap(F[i],F[rev[i]]);
    for(register int i=2,ii=1,t=1;i<=Len;i<<=1,ii<<=1,++t)
        for(register int j=0;j<Len;j+=i)rep(k,0,ii)
        {
            register int tt=(uint64)F[j+k+ii]*g[t][k]%mod;
            F[j+k+ii]=ad(F[j+k],mod-tt);
            F[j+k]=ad(F[j+k],tt);
        }
    if(typ==-1)
    {
        reverse(F+1,F+Len);
        register uint64 invn=power(Len,mod-2);
        rep(i,0,Len)F[i]=invn*F[i]%mod;
    }
}

int X[MAXN],Y[MAXN];

inline void mul(int *a,int *b,int *c,int lenl,int lenr)
{
    if(lenl<=300/lenr)
    {
        Rep(i,0,lenl+lenr)c[i]=0;
        Rep(i,0,lenl)Rep(j,0,lenr)c[i+j]=(c[i+j]+a[i]*b[j])%mod;
        return;
    }
    for(Len=2;Len<=lenl+lenr;Len<<=1);
    Rep(i,0,lenl)X[i]=a[i];
    Rep(i,0,lenr)Y[i]=b[i];
    rep(i,lenl+1,Len)X[i]=0;
    rep(i,lenr+1,Len)Y[i]=0;
    NTT(X,1),NTT(Y,1);
    rep(i,0,Len)X[i]=(ll)X[i]*Y[i]%mod;
    NTT(X,-1);
    Rep(i,0,lenl+lenr)c[i]=X[i];
}

注意 $FFT$ 的复数类最好自己手写,比系统自带的快多了。

以下几个操作大多基于取模来实现,存在一个特殊的模数,满足其减 $1$ 后存在较多 $2$ 的因数。无说明一般为 $998244353$ 。

2.多项式求逆

点我

也算是基础吧。

推导过程:

设要求 $F$ 关于 $x^n$ 的逆,则

$f_i(x)\equiv F^{-1}(x)\pmod{x^{2^i}}$

假设已知 $f_i$ ,要求 $f_{i+1}$ ,那么

$F(x)f_i(x)\equiv F(x)f_{i+1}(x)\pmod{x^{2^i}}$

$[f_i(x)-f_{i+1}(x)]\equiv 0\pmod{x^{2^i}}$

两边平方,则

$[f_i(x)-f_{i+1}(x)]^2\equiv 0\pmod{x^{2^{i+1}}}$

$f_i^2(x)-2f_i(x)f_{i+1}(x)+f_{i+1}^2(x)\equiv 0\pmod{x^{2^{i+1}}}$

同乘$F$,得

$f_{i+1}(x)-2f_i(x)+F(x)f_i^2(x)\equiv 0\pmod{x^{2^{i+1}}}$

$f_{i+1}(x)\equiv 2f_i(x)-F(x)f_i^2(x)\pmod{x^{2^{i+1}}}$

利用这个公式直接倍增就好了。时间复杂度利用主定理可以证明为 $O(n\log n)$ 。

代码:


inline void Inv(int*F,int*G,int ln)
{
    Iv[0]=power(F[0],mod-2);
    rep(i,0,ln)X[i]=Y[i]=0;
    for(register int Ln=2;Ln<=ln;Ln<<=1)
    {
        rep(i,0,Ln)X[i]=F[i];
        rep(i,0,(Ln>>1))Y[i]=Iv[i];
        Len=Ln,calrev();
        NTT(X,1),NTT(Y,1);
        rep(i,0,Ln)X[i]=(uint64)X[i]*Y[i]%mod;
        NTT(X,-1);
        rep(i,0,(Ln>>1))X[i]=0;
        NTT(X,1);
        rep(i,0,Ln)X[i]=(uint64)X[i]*Y[i]%mod;
        NTT(X,-1);
        rep(i,(Ln>>1),Ln)Iv[i]=mod-X[i];
    }
    rep(i,0,ln)G[i]=Iv[i];
}

3.多项式除法

点我

已知函数$F$和$G$,求$Q$和$R$,使得

$F=QG+R$

保证 $Deg(R)<Deg(G)<Deg(F)$ 。

推导过程:

先考虑一种操作:

$F_R(x)=x^{Deg(F)}F(\frac{1}{x})$

拆开发现该操作的作用就是将$F$的系数反位。

然后开始化式子。

为了简便,设$Deg(F)=n,Deg(G)=m$,那么可知$Deg(Q)=n-m,Deg(R)=m-1$。

$\because F(x)=G(x)Q(x)+R(x)$

$\therefore F(\frac{1}{x})=G(\frac{1}{x})Q(\frac{1}{x})+R(\frac{1}{x})$

$\therefore x^nF(\frac{1}{x})=x^nG(\frac{1}{x})Q(\frac{1}{x})+x^nR(\frac{1}{x})$

$\therefore x^nF(\frac{1}{x})=x^mG(\frac{1}{x})x^{n-m}Q(\frac{1}{x})+x^{n-m+1}*x^{m-1}R(\frac{1}{x})$

$\therefore F_R(x)=G_R(x)Q_R(x)+x^{n-m+1}*R_R(x)$

将这个式子放在模$n-m+1$次方下,可以发现$R_R$的所有系数正好全部被模掉,失去影响。

$\therefore F_R(x)=G_R(x)Q_R(x)\pmod{x^{n-m+1}}$

$\therefore Q_R(x)=F_R(x)G_R^{-1}(x)\pmod{x^{n-m+1}}$

多项式求逆即可。至于$R$,直接$R=F-QG$就可以求出。时间复杂度 $O(n\log n)$ 。

代码:

static int X[MAXN],Y[MAXN];

inline void Div(int *a,int n,int *b,int m)
{
    memcpy(X,a,sizeof X);memcpy(Y,b,sizeof Y);
    reverse(b,b+m+1);reverse(X,X+n+1);
    Rep(i,n-m+1,n)X[i]=0;
    Inv(b,n-m+1);
    while(Len<=(n<<2))Len<<=1;
    calrev(Len);
    mul(X,b);reverse(X,X+n-m+1);
    Rep(i,n-m+1,n)X[i]=0;
    mul(Y,X);
    Rep(i,0,m-1)b[i]=ad(a[i],mod-Y[i]);
    memcpy(a,X,sizeof X);
}

4.多项式开方

点我

开始有点麻烦了。依旧采用倍增来求。

已知$F$,求$G$使得$G^2\equiv F\pmod{x^n}$

假设求出 $T_0^2\equiv F\pmod{x^{t}}$ ,要求 $T^2\equiv F\pmod{x^{2t}}$ 。

$T_0^2-T^2\equiv 0\pmod{x^t}$

$(T_0+T)(T_0-T)\equiv 0\pmod{x^t}$

因为存在平方,所以存在2个解。设 $T-T_0\equiv 0\pmod{x^{t}}$ ,则再平方有

$T'^2-2T_T+T^2\equiv 0\pmod{x^{2t}}$

又因为 $T^2\equiv F\pmod{x^{2t}}$ ,则

$T'^2-2TT_0+F\equiv 0\pmod{x^{2t}}$

$T\equiv \frac{T_0^2+F}{2T_0}\pmod{x^{2t}}$

直接倍增计算即可。时间复杂度依旧 $O(n\log n)$ 。

代码:

inline void Sqrt(int*F,int*G,int ln)
{
    G[0]=1;
    for(register int Ln=2;Ln<=ln;Ln<<=1)
    {
        mul(G,G,ExX,Ln>>1,Ln>>1);
        rep(i,0,Ln>>1)ExX[i]=ad(mod-ExX[i+(Ln>>1)],F[i+(Ln>>1)]);
        rep(i,Ln>>1,Ln)ExX[i]=0;
        Inv(G,ExY,Ln>>1);
        Len=Ln,calrev();
        NTT(ExX,1),NTT(ExY,1);
        rep(i,0,Ln)ExX[i]=(ll)ExX[i]*ExY[i]%mod;
        NTT(ExX,-1);
        rep(i,0,Ln>>1)G[i+(Ln>>1)]=(ll)ExX[i]*inv2%mod;
    }
}

5.多项式积分和求导

这就不要例题了吧...

设$F(x)=\sum\limits_{i=0}^n a_ix^i$,则

定义:

$\int F(x)=\sum\limits_{i=0}^n \frac{a_i}{i+1}x^{i+1}$

$\frac{dF}{dx}=\sum\limits_{i=1}^n ia_ix^{i-1}$

直接$O(n)$计算就可以了。

代码:

inline void Direv(int *a,int n)
{Rep(i,1,n)a[i-1]=(ll)a[i]*i%mod;a[n]=0;}

inline void Inter(int *a,int n)
{Repe(i,n,1)a[i]=(ll)a[i-1]*power(i,mod-2)%mod;a[0]=0;}

6.多项式取 $\ln$

模板题点我

做完模板可以点我

$\ln$ 的定义为 $\exp$ 的逆变换。具体定义请参考下面的 $\exp$ 。

定义式为:$$ \ln(1-x)=-\sum_{i=1}^\infty \frac{x^i}{i} $$ 存在推导式:

$$\ln F=\int \frac{F'}{F}\equiv \int F'F^{-1}$$

直接计算即可。

代码:

inline void Ln(int *a,int n)
{
    memcpy(X,a,sizeof X);
    Inv(X,n);Direv(a,n);
    while(Len<=(n<<2))Len<<=1;
    calrev(Len);
    mul(a,X);
    Inter(a,n);
    Rep(i,n+1,Len)a[i]=0;
}

7.多项式牛顿迭代

题目请见 $WC2019$ 汪乐平的课件。

采用倍增实现。

已知$G$,求多项式$F$使得$G(F(x))= 0$。

设已知$F_0$满足$G(F_0(x))\equiv 0\pmod{x^t}$,要求$G(F(x))\equiv 0\pmod{x^{2t}}$。

存在递推式

$F(x)≡F_0(x)-\frac{G(F_0(x))}{G'(F_0(x))} \pmod {x^{2t}}$

其中,$G'(F(x))=\frac{\mathbb{d}G}{\mathbb{d}F}$,即 $G$ 对 $F$ 的导数。

直接计算即可。

7.1推导过程:

函数 $G$ 对 $F$ 进行泰勒展开,得:

$$G(F)=G(F_0)+G'(F_0)(F-F_0)+G''(F_0)(F-F_0)^2+\cdots$$

因为$F\equiv F_0\pmod{x^t}$,因此$F-F_0\equiv 0\pmod{x^t}$,因此$(F-F_0)^k\equiv 0\pmod{x^{2t}},k\ge 2$。

因此后面几项在$\bmod x^{2t}$意义下全部为$0$,直接忽略。

因此得到关系式

$$G(F)=G(F_0)+G'(F_0)(F-F_0)$$

整理即得上式。

(代码依 $G$ 不同而不同)

7.2.其他式子利用牛顿迭代推导的方式

7.2.1 求逆

考虑对方程$F(x)G(x)-1=0$ 作牛顿迭代,得:

$\delta(F)=GF-1$

$\delta'(F)=G$

可以得到递推式:$F=F_0-\frac{\delta(F_0)}{\delta'(F_0)}$

$=F_0-\frac{G_0F_0-1}{G_0}$

$=F_0-F_0(G_0F_0-1)$

$=2F_0-G_0F_0^2$

7.2.2 开方

考虑对方程$F^2-G=0$作牛顿迭代,得:

$\delta(F)=F^2-G$

$\delta'(F)=2F$

可以得到递推式:

$F=F_0-\frac{\delta(F_0)}{\delta'(F_0)}$$=F_0-\frac{F_0^2-G_0}{2F_0}$$=\frac{F_0^2+G_0}{2F_0}$

7.2.3 $\exp$

见下文。

7.3 牛顿法的通用优化技巧

使用牛顿法,我们通常需要多次求出$F_0-\frac{G(F_0)}{G'(F_0)}$。有时为了方便,我们会将前面的$F_0$和后面部分共同计算。

可以发现,因为$F\equiv F_0\pmod{x^{\frac{n}{2}}}$,所以后面部分一定有长度为$\frac{n}{2}$的部分全部都是$0$。此时,我们可以将后半部分位移至前面,然后再进行操作。通常可以将多项式乘法的长度减少一半,时间复杂度会大大优化。

例如,多项式开方的公式为$F=F_0+\frac{1}{2}F_0^{-1}(G-F_0^2)$。应用以上方法可以知道,$(G-F_0^2)$有一半都是$0$。那么,如果当前层长度为$n$,从$\frac{n}{2}$转移过来,此时我们做多项式乘法只需要做长度为$n$的而不是$2n$的。同理,求逆也只需要处理到$n$。这样可以省下近一半常数。

8.多项式 $\exp$

点我

由泰勒展开可得 $e^x=\sum_{i=0}^\infty \frac{x^i}{i!}$。因此,我们定义 $\sum_{i=0}^\infty \frac{F(x)^i}{i!}=\exp{F}$。介于 $e$ 在膜域下无定义,一般默认$F(x)$ 常数项为 $0$ 。

已知函数 $F$ ,求 $G\equiv \exp{F}$ 。保证 $F$ 常数项为 $0$ 。

可以知道 $\exp(\ln P)=P$ 。

所以 $\ln G-F=0$ 。

直接套用上方多项式牛顿迭代的公式,可得

$G=G_0(1-\ln G_0+F)$

倍增即可。时间复杂度 $O(n\log n)$ 。

该算法常数巨大,经验算大约在 $70\sim 80$ 左右。因此可以粗略认为是 $O(n\log^2 n)$ 的。

最近发现常数过大还没有分治 $FFT$ 跑得快。

代码:

static int ExX[MAXN],ExY[MAXN];

inline void Exp(int *F,int *G,int ln)
{
    G[0]=1;
    for(register int Lx=2;Lx<=ln;Lx<<=1)
    {
        rep(i,Lx>>1,Lx)G[i]=0;
        Ln(G,ExX,Lx);
        rep(i,0,Lx>>1)ExX[i]=ad(F[i+(Lx>>1)],mod-ExX[i+(Lx>>1)]);
        rep(i,0,Lx>>1)ExY[i]=G[i];
        rep(i,Lx>>1,Lx)ExX[i]=0;
        Len=Lx;
        calrev();
        NTT(ExY,1),NTT(ExX,1);
        rep(i,0,Len)ExX[i]=(ll)ExX[i]*ExY[i]%mod;
        NTT(ExX,-1);
        rep(i,0,Len>>1)G[i+(Len>>1)]=ExX[i];
    }
}

9.多项式幂次

例题见10.拉格朗日反演。

直接上式子吧...

当 $F[0]=1$ 时:

$$F^k=\exp(k\ln F)$$

否则,设最低项为 $a_tx^t$ ,则

$$F^k=a_t^kx^{tk}\exp(k\ln \frac{F}{a_tx^t})$$

直接按照式子计算即可。时间复杂度$O(n\log n)$ 。

注意:

直接求 $\ln$ 再 $\exp$ 的方法仅适用于常数项为1的情况(因为多项式 $\exp$ 计算时默认常数项为 $1$ )。因此当常数项不为 $1$ 时需要将多项式的常数项强制变成 $1$ 。多判一下即可(不过大部分时间都不用管)。

代码:

inline void Pow(int*F,int*G,int ln,ll k)
{
    int lst=ln+1;
    Rep(i,0,ln)if(F[i]){lst=i;break;}
    if(ln&&(__int128)lst*k>ln)
    {memset(G,0,sizeof(int)*(ln+1));return;}
    int md=ln-k*lst,bs=F[lst],iv=power(bs,mod-2);
    Rep(i,lst,ln)G[i]=(ll)G[i]*iv%mod;
    Ln(G+lst,G,md);
    k%=mod;
    Rep(i,0,md)G[i]=(ll)G[i]*k%mod;
    Exp(G,G,md);
    bs=power(bs,k);
    Repe(i,md,0)G[i+lst*k]=(ll)G[i]*bs%mod;
    memset(G,0,sizeof(int)*(lst*k));
}

9.5.多项式 $k$ 次方根

在膜意义下,存在一下转化:

$$\sqrt[k]F=F^{\frac{1}{k}}=F^{k^{-1}}$$

直接套用上面的幂次即可。复杂度为 $O(n\log n)$ 。

代码就咕了。

10.拉格朗日反演

点我

已知 $G$ ,求函数 $F$ 使得 $F[G(x)]=x$ 。

经过复杂的推理,存在关系式:

$$[x^n] F(x)=\frac{1}{n}[x^{n-1}](\frac{x}{G(x)})^n$$

其中$[x^n]$表示多项式$x^n$项系数。

直接按照式子模拟即可。时间复杂度 $O(n\log n)$ 。

代码:

inline int Lagrange(int *a,int n)
{
    Rep(i,0,n-1)a[i]=a[i+1];
    Inv(a,n);
    pow(a,n,n);
    return (ll)a[n-1]*power(n,mod-2)%mod;
}

10.5 扩展拉格朗日反演

其实就是另一个式子。

$$[x^n]H[F(x)]=\frac{1}{n}[x^{n-1}]H'(x)(\frac{x}{G(x)})^n$$

11.分治FFT

点我

已知 $F[0]=1$ ,求 $F$ 满足$$F[i]=\sum_{j=1}^i F[i-j]*G[j]$$

其实哪怕是$F[i]=aF[i-b]+c+\sum_{j=1}^i dF[i-j]*G[j]$也是可以解决的。

考虑利用分治。

假设我们求出了$l\to mid$的答案,要求这些点对$mid+1\to r$的影响,那么对右半边点$x$的贡献为:

$w_x=\displaystyle\sum_{i=l}^{mid} f[i]g[x-i]$

这部分可以利用卷积来快速计算。计算完以后,答案直接加到答案数组就可以了。

需要注意的是,如果要求左边点对右边点的影响,首先整个区间以左对该区间的贡献应该先求出。所以分治过程为先分治左边,在求出中间,然后在递归右边。

时间复杂度 $O(n\log^2n)$ 。

有一个简单卡常方法。求一层的时候,我们会做长度为 $md-l+1+r-l+1\approx 1.5(r-l+1)$ 的卷积,但是这没有必要。一个显然的观察是我们只会使用在 $md-l+1\sim r-l+1$ 内的值,而前半截内容显然是没有用的。因此我们可以使用长度为 $r-l+1$ 的卷积。一般情况下可以快 $50\%$ 。

代码(求 $\exp$):

void cdq_FFT(int l,int r)
{
    if(l==r)
    {
        if(!l)F[l]=1;
        else F[l]=(ll)F[l]*iv[l]%mod;
        return;
    }
    int md=(l+r)>>1;cdq_FFT(l,md);
    for(Len=2;Len<=r-l;Len<<=1);
    calrev();
    memset(X,0,sizeof(int)*Len);
    memset(Y,0,sizeof(int)*Len);
    memcpy(X,F+l,sizeof(int)*(md-l+1));
    memcpy(Y,a,sizeof(int)*(r-l));
    NTT(X,1),NTT(Y,1);
    rep(i,0,Len)X[i]=(ll)X[i]*Y[i]%mod;
    NTT(X,-1);
    Rep(i,md+1,r)F[i]=ad(F[i],X[i-l-1]);
    cdq_FFT(md+1,r);
}

11.5.特殊的分治FFT

点我

先假设你会前面的推式子部分。

你现在得到了一个长得很丑的式子,大概长这样:$$f_n=(n-1)f_{n-1}+\sum_{i=1}^{n-2}(i-1)f_if_{n-i}$$ 换句话说,你得到了一个类似于$F=F+F*F$ 的式子。这个式子利用上面的分治 $FFT$ 方法当然是求不出的了。

这时候,你就需要对你的分治 $FFT$ 进行魔改了。魔改结果大概长这样:

void cdqFFT(int *f,int l,int r)
{
    if(l==r)
    {
        f[l]=ad(f[l],(uint64)(l-1)*f[l-1]%mod);
        return;
    }
    int mid=(l+r)>>1;
    cdqFFT(f,l,mid);
    for(Len=r-l+1;Len>(Len&-Len);Len+=(Len&-Len));
    Rep(i,0,Len)A[i]=B[i]=0;
    Rep(i,l,mid)A[i-l]=(uint64)f[i]*(i-1)%mod,B[i-l]=f[i];
    calrev(Len);mul(A,B);
    Rep(i,mid+1,r)if(i>=l*2)f[i]=ad(f[i],A[i-2*l]);
    if(l^2)
    {
        Rep(i,0,Len)A[i]=B[i]=0;
        Rep(i,2,min(l-1,r-l))A[i-2]=f[i];
        Rep(i,l,mid)B[i-l]=f[i];
        mul(A,B);
        Rep(i,mid+1,r)if(i>=l+2)
            f[i]=ad(f[i],(uint64)(i-2)*A[i-l-2]%mod);
    }
    cdqFFT(f,mid+1,r);
}

稍微解释一下。

首先,前面的 $(n-1)f_{n-1}$ 明显是可以在递归最底层在加上的。这个很好处理。

问题就在于后面的 $\sum_{i=1}^{n-2}(i-1)f_if_{n-i}$。这部分很难求。但是你可以分两种情况来进行讨论:

  1. 所计算区间左端点是全局左端点。原来怎么做就怎么做。但是注意不要在卷积中用到还没有计算完毕的项。
  2. 所计算区间不是全局左端点。可以知道会存在一些项之前并没有计算过。我们将这些项合并一起计算。具体来说就是把 $(i-1)f_if_{n-i}$ 和 $(n-i-1)f_{n-i}f_i$ 一起计算。那么这两项加在一起就变成了$(n-2)f_if_{n-i}$ 。直接加上即可。

其余部分和一般分治$FFT$相同。注意处理几个边界情况。

注意,我程序内部有一个 $if(l\otimes 2)$ ,这是判断左端点是否为全局左端点的。因为这道题的全局左端点是 $2$ 。

12.常系数齐次线性递推

点我

主要就是推式子。

推导过程和代码见线性递推

13.多项式求三角函数

经过审稿管理员提示,我发现原来这东西还是有用的。。。

已知$F(x)$,求$p[F(x)]$。其中, $p(x)$ 为某一个三角函数 $(\sin,\cos,\tan,\sec,\csc,\cot)$ 。

首先,可以知道,如果求出了$\sin$ 和$\cos$ ,那么其它的三角函数也可以轻易推出。那么问题就是如何求出这两个三角函数。

根据欧拉定理可知:

$$e^{w_4x}=\cos(x)+w_4\sin(x)$$

其中, $w_4$ 为四次单位根,也就是常用的 $i$ 。

那么同理可得

$$e^{w_4F(x)}=\cos[F(x)]+w_4\sin[F(x)]$$

将$F(x)$ 的系数都乘以 $w_4$ 然后 $\exp$ 就可以了。时间复杂度 $O(n\log n)$ 。

注意以上的内容是 $FFT$ 的求法。如果要用 $NTT$ 的话也不是不行。

$$\cos(x)=\sum_{i=0}^\infty \frac{[i\bmod 2=0](-1)^\frac{i}{2}x^i}{i!}$$

$$\cosh(x)=\frac{e^x+e^{-x}}{2}=\sum_{i=0}^\infty \frac{[i\bmod 2=0]x^i}{i!}$$

可得

$$\cos(x)=\frac{e^{w_4x}+e^{-w_4x}}{2}$$

$$\cos(F)=\frac{1}{2}(e^{w_4F}-e^{-w_4F})$$

因为 $w_4$ 在常用 $NTT$ 模数下一般都有定义,且实系数函数的 $\cos$ 肯定也是实系数函数,因此这个公式计算出来肯定不含复系数。

至于 $\sin$ ,可以利用类似的方法来计算。即:

$$\sin(x)=\sum_{i=0}^\infty \frac{[i\bmod2=1](-1)^{\frac{i-1}{2}}x^i}{i!}$$

$$\sinh(x)=\frac{e^x-e^{-x}}{2}=\sum_{i=0}^\infty \frac{[i\bmod2=1]x^i}{i!}$$

推得:

$$\sin(x)=\frac{w_4^3}{2}(e^{w_4x}-e^{-w_4x})$$

同理得:

$$\sin F=\frac{w_4^3}{2}(e^{w_4F}-e^{-w_4F})$$

代码:

static int Sn[MAXN],Cs[MAXN];

inline void Cos(int*F,int*G,int ln)
{
    Rep(i,0,ln)Cs[i]=(ll)F[i]*w4%mod;
    Exp(Cs,Cs,ln),Inv(Cs,G,ln);
    Rep(i,0,ln)G[i]=((ll)Cs[i]+G[i])*inv2%mod;
}

const int w43div2=(ll)w4*w4%mod*w4%mod*inv2%mod;

inline void Sin(int*F,int*G,int ln)
{
    Rep(i,0,ln)Sn[i]=(ll)F[i]*w4%mod;
    Exp(Sn,Sn,ln),Inv(Sn,G,ln);
    Rep(i,0,ln)G[i]=((ll)Sn[i]+mod-G[i])*w43div2%mod;
}

14.多项式反三角函数

这个东西和上面的没有半毛钱关系。。。

根据积分表可以知道:

$$\arcsin x = \int \frac{1}{\sqrt{1-x^2}}$$

$$\arctan x = \int \frac{1}{1+x^2}$$

那么同理可得:

$$\arcsin F = \int \frac{F'}{\sqrt{1-F^2}}$$

$$\arctan F = \int \frac{F'}{1+F^2}$$

直接计算即可。

代码:

inline void Arcsin(int*F,int*G,int ln)
{
    memcpy(Sn,F,sizeof(int)*(ln+1));
    for(Len=2;Len<=ln*2;Len<<=1);calrev();
    NTT(Sn,1);rep(i,0,Len)Sn[i]=(ll)Sn[i]*Sn[i]%mod;
    NTT(Sn,-1);
    Sn[0]=1;rep(i,ln+1,Len)Sn[i]=0;
    Rep(i,1,ln)Sn[i]=mod-Sn[i];
    Sqrt(Sn,Sn,ln);
    Inv(Sn,Sn,ln);
    Deriv(F,G,ln);
    mul(G,Sn,G,ln,ln);
    Rep(i,ln+1,ln<<1)G[i]=0;
    Inter(G,G,ln);
}

inline void Arctan(int*F,int*G,int ln)
{
    memcpy(Sn,F,sizeof(int)*(ln+1));
    for(Len=2;Len<=ln*2;Len<<=1);calrev();
    NTT(Sn,1);rep(i,0,Len)Sn[i]=(ll)Sn[i]*Sn[i]%mod;
    NTT(Sn,-1),Sn[0]=1;
    rep(i,ln+1,Len)Sn[i]=0;Inv(Sn,Sn,ln);
    Deriv(F,G,ln),mul(G,Sn,G,ln,ln);
    Rep(i,ln+1,ln<<1)G[i]=0;
    Inter(G,G,ln);
}

15.多项式多点求值

点我

已知函数 $F$ 和数组 $\{a_i\}$ ,求 $\{F(a_i)\}$ 。

直接计算是 $O(nm)$ 的。此时,我们可以考虑降幂。

构造一个函数 $D_{l,r}=\prod_{i=l}^r (x-a_i)$ 。

可以发现,$\forall k\in [l,r],D_{l,r}(a_k)=0 $。

因此,$\forall k\in [l,r],F(a_k)=F(a_k)\bmod D_{l,r}(a_k) $。

于是,设 $F_{l,r}(x)$ 表示 $a_l$ 到 $a_r$ 在原函数上的点构成的多项式,则:

$$F_{l,mid}(x)=F_{l,r}(x)\bmod D_{l,mid}(x)$$

直接分治即可。时间复杂度为 $O(n\log^2 n)$ 。注意 $D$ 的计算可以直接分治乘法每次乘起来,不影响复杂度。

代码:

static int solv[18][2][MAXN];

inline void calc(int*a,int l,int r,int lev,int dir)
{
    if(l==r)
    {
        solv[lev][dir][0]=mod-a[l];
        solv[lev][dir][1]=1;return;
    }
    int md=(l+r)>>1;
    calc(a,l,md,lev+1,0),calc(a,md+1,r,lev+1,1);
    mul(solv[lev+1][0],solv[lev+1][1],solv[lev][dir],md-l+1,r-md);
}

static int pol[18][MAXN],Q[MAXN];

void getnum(int*a,int*ans,int l,int r,int lev,int n)
{
    if(n>r-l)
    {
        calc(a,l,r,0,0);
        Div(pol[lev-(lev>0)],solv[0][0],Q,pol[lev],n,r-l+1);
        n=r-l;
    }
    else if(lev){Rep(i,0,n)pol[lev][i]=pol[lev-1][i];}
    if(l==r)
    {
        ans[l]=pol[lev][0];
        return;
    }
    int md=(l+r)>>1;
    getnum(a,ans,l,md,lev+1,n);
    getnum(a,ans,md+1,r,lev+1,n);
}

16.多项式快速插值

点我

给定二维平面上点集$\{(x_i,y_i)\}$(保证$x_i$互不相同),求函数$F(x)$,满足所有点均在该函数上。

核心内容是优化拉格朗日插值。

具体过程见多项式快速插值

Ex1.MTT

点我还有

上面的所有函数基于模数为一个 $NTT$ 模数(即 $a*2^k+1$ ,其中$k$比较大)。如果不是的话,可以将多项式乘法换成以下的函数:

inline void mul(int*F,int*G)
{
    rep(i,0,Len)A[i]=comp(F[i]>>15,F[i]&32767),
        C[i]=comp(G[i]>>15,G[i]&32767);
    DBLFFT(A,B,1),DBLFFT(C,D,1);
    rep(i,0,Len)
    {
        register comp t=A[i];
        A[i]=B[i]*C[i]+A[i]*D[i];
        B[i]=B[i]*D[i];
        C[i]=C[i]*t;
    }
    DBLFFT(A,B,-1),FFT(C,-1);
    rep(i,0,Len)
    {
        register int lz=(((ll)floor(C[i].re+0.5))%mod*73741817
             +((ll)floor(A[i].re+0.5))%mod*32768+(ll)floor(B[i].re+0.5))%mod;
        F[i]=i<=k?lz:0;
    }
}

其中 $DBLFFT$ 是共轭优化过的 $FFT$ 。下面一章有代码。

简单解释其实就是将系数拆成 $2^{15}a+b$ ,然后利用等式 $(2^{15}A+B)(2^{15}C+D)=2^{30}AC+2^{15}(BC+AD)+BD$ 直接计算。原理为 $double$ 的精度可以承受住 $10^{14}$ ,因此能够支撑 $2^{15}*2^{15}*10^5=107374182400000$ 以内的数字运算。

具体还是看 $2016$ 毛爷爷写的国家队论文吧。

顺便一提,配合共轭优化后可以做到 $4$ 次 $DFT$ ,基本和一般的乘法没什么区别了。

Ex2.BlueStein算法

这个算例题吧

已知多项式 $F(x)$ ,求 $\forall i\in[0,k),A_i=F(w_k^i)$ 。其中 $w_k$ 为 $k$ 次单位根。

被这个东西送到了退役边缘。。。只怪我自己没有认真看论文。。。

因为 $w_k^t=w_k^{t\bmod k}$ ,因此可以认为 $F$ 的最高次幂为 $k-1$ 。

设 $F$ 的 $x^t$ 次方系数为 $f_t$ 。

那么,将式子写出,得:

$$F(w_k^t)=\sum_{i=0}^{k-1}f_iw_k^{ti}$$

$$=\sum_{i=0}^{k-1}f_iw_k^{-(k-t)i}$$

设这个东西为 $ans_{k-t}$ ,则 $F(w_k^t)=ans_{k-t}$ 。

运用套路,得:

$$ti=\binom{t+i}{2}-\binom{t}{2}-\binom{i}{2}$$

因此有:

$$ans_t=\sum_{i=0}^{k-1} f_iw_k^{\binom{i-t}{2}-\binom{-t}{2}-\binom{i}{2}}$$

$$=w_k^{-\binom{-t}{2}}\sum_{i=0}^{k-1} f_iw_k^{-\binom{i}{2}}*w_k^{\binom{-(t-i)}{2}}$$

因为我们每一项都要卷到 $k-1$ ,因此考虑将系数直接后移 $k-1$ 位,可以得到新的式子:

$$ans_t=w_k^{-\binom{-t}{2}}\sum_{i=0}^{k-1+t} f_iw_k^{-\binom{i}{2}}*w_k^{\binom{k-1-(k-1+t-i)}{2}}$$

我们设 $A_i=f_iw_k^{-\binom{i}{2}}$ , $B_i=w_k^{\binom{k-1-(k-1+t-i)}{2}}$ ,则

$$ans_t=w_k^{-\binom{-t}{2}}[k-1+t](A*B)$$

时间复杂度 $O(n\log n)$ 。

这个过程被称为 $CZT$ 。 $ICZT$ 只需要将 $F_1$ 到 $F_{k-1}$ 的部分翻转,然后整体除掉 $k$ 即可。

算出这个东西后,我们就可以将他当做以 $w_k$ 为公比的插值。因此如果将多个式子相卷,得到的结果将会是这些式子在膜 $x^k$ 意义下的循环卷积。有的时候会有用。

代码:

inline void Blue_Stein(int*F,int*G,int ln)
{
    rep(i,0,ln)A[i]=(ll)F[i]*power(ome,ln-(ll)i*(i-1)/2%ln,mod)%mod;
    Rep(i,0,2*ln-2)B[i]=power(ome,ln+(ll)(ln-1-i)*(ln-2-i)/2%ln,mod);
    mul(A,B,A,ln-1,2*ln-2);
    reverse(A+1,A+k);
    rep(i,0,ln)G[i]=(ll)A[i+k-1]*power(ome,ln-(ll)(-i)*(-i-1)/2%ln,mod)%mod;
}

注意使用的时候大部分情况下模数都不会很友好,因此结合上面的 $MTT$ 或者直接上 $FFT$ 是一般情况。

Ex3.常数优化

根据经验可以发现,多项式方面的常数差异非常大,有的时候别人跑过了但是你没跑过就凉了。因此,常数优化是非常重要的。限于篇幅就只推荐这篇这篇了。后面那篇的详细内容和不详细证明可以参考毛爷爷的国家队论文。

此处放出共轭优化的代码。该代码可以同时对 $F$ 和 $G$ 进行 $DFT$ 或者 $IDFT$ 操作。


static comp Z[MAXN];

inline void DBLFFT(comp*F,comp*G,int tl)
{
    if(tl==1)
    {
        memcpy(Z,F,sizeof(comp)*Len);
        FFT(Z,1);
        memcpy(F,Z,sizeof(comp)*Len);
        reverse(F+1,F+Len);
        rep(i,0,Len)F[i].im=-F[i].im;
        rep(i,0,Len)G[i]=(Z[i]-F[i])/2,F[i]=(Z[i]+F[i])/2
            ,swap(G[i].re,G[i].im),G[i].im=-G[i].im;
    }
    else
    {
        rep(i,0,Len)F[i]=comp(F[i].re-G[i].im,F[i].im+G[i].re);
        FFT(F,-1);
        rep(i,0,Len)G[i]=comp(F[i].im,0),F[i]=comp(F[i].re,0);
    }
}

Ex4.快速子集变换

点我

和 $FFT$ 的形式差不多。代码和内容在这里

Ex4.5快速子集运算

限于篇幅在另一篇文章作说明。

总集:帕秋莉的超级多项式

已知 $n-1$ 次多项式 $F$ 和 $k$ ,求

$$\left[\left(1+\ln(1+\frac{1}{\exp(\int\frac{1}{\sqrt{F}})})\right)^k\right]'$$

这道题汇集了几乎所有的多项式基础操作,有兴趣的可以尝试做一下,练练手。

代码:

#include<bits/stdc++.h>
#define Rep(i,a,b) for(register int i=(a);i<=(b);++i)
#define Repe(i,a,b) for(register int i=(a);i>=(b);--i)
#define rep(i,a,b) for(register int i=(a);i<(b);++i)
#define pb push_back
#define mx(a,b) (a>b?a:b)
#define mn(a,b) (a<b?a:b)
typedef unsigned long long uint64;
typedef unsigned int uint32;
typedef long long ll;
using namespace std;

namespace IO
{
    const uint32 Buffsize=1<<15,Output=1<<24;
    static char Ch[Buffsize],*S=Ch,*T=Ch;
    inline char getc()
    {
        return((S==T)&&(T=(S=Ch)+fread(Ch,1,Buffsize,stdin),S==T)?0:*S++);
    }
    static char Out[Output],*nowps=Out;

    inline void flush(){fwrite(Out,1,nowps-Out,stdout);nowps=Out;}

    template<typename T>inline void read(T&x)
    {
        x=0;static char ch;T f=1;
        for(ch=getc();!isdigit(ch);ch=getc())if(ch=='-')f=-1;
        for(;isdigit(ch);ch=getc())x=x*10+(ch^48);
        x*=f;
    }

    template<typename T>inline void write(T x,char ch='\n')
    {
        if(!x)*nowps++='0';
        if(x<0)*nowps++='-',x=-x;
        static uint32 sta[111],tp;
        for(tp=0;x;x/=10)sta[++tp]=x%10;
        for(;tp;*nowps++=sta[tp--]^48);
        *nowps++=ch;
    }
}
using namespace IO;

void file(void)
{
    FILE*DSD=freopen("polynomial.in","r",stdin);
    FILE*CSC=freopen("polynomial.out","w",stdout);
}

const int MAXN=1<<21;

namespace poly
{
    const int mod=998244353,gen=3;

    static int g[23],iv[MAXN];

    inline int power(int u,int v)
    {
        register int sm=1;
        for(;v;v>>=1,u=(uint64)u*u%mod)if(v&1)
            sm=(uint64)sm*u%mod;
        return sm;
    }

    inline void predone()
    {
        Rep(i,1,21)g[i]=power(gen,(mod-1)>>i);
        iv[1]=1;
        Rep(i,2,1e5)iv[i]=mod-(uint64)mod/i*iv[mod%i]%mod;
    }

    static int Len,rev[MAXN];

    inline void calrev()
    {
        int II=log(Len)/log(2)-1;
        Rep(i,1,Len-1)rev[i]=rev[i>>1]>>1|(i&1)<<II;
    }

    inline int ad(int u,int v){return(u+=v)>=mod?u-mod:u;}

    inline void NTT(int*F,int typ)
    {
        Rep(i,1,Len-1)if(i<rev[i])swap(*(F+i),*(F+*(rev+i)));
        for(register int i=2,ii=1,t=1;i<=Len;i<<=1,ii<<=1,++t)
        {
            register int wn=*(g+t);
            for(register int j=0;j<Len;j+=i)
            {
                register int w=1;
                rep(k,0,ii)
                {
                    register int tt=(uint64)*(F+j+k+ii)*w%mod;
                    *(F+j+k+ii)=ad(*(F+j+k),mod-tt);
                    *(F+j+k)=ad(*(F+j+k),tt);
                    w=(uint64)w*wn%mod;
                }
            }
        }
        if(typ==-1)
        {
            reverse(F+1,F+Len);
            register uint64 invn=power(Len,mod-2);
            rep(i,0,Len)*(F+i)=invn**(F+i)%mod;
        }
    }

    static int X[MAXN],Y[MAXN],Iv[MAXN];

    inline void mul(int *a,int *b,int *c,int lenl,int lenr)
    {
        if((ll)lenl*lenr<=300)
        {
            Rep(i,0,lenl+lenr)X[i]=0;
            Rep(i,0,lenl)Rep(j,0,lenr)
                X[i+j]=(X[i+j]+(ll)a[i]*b[j])%mod;
            Rep(i,0,lenl+lenr)c[i]=X[i];
            return;
        }
        for(Len=2;Len<=lenl+lenr;Len<<=1);
        calrev();
        Rep(i,0,lenl)X[i]=a[i];
        Rep(i,0,lenr)Y[i]=b[i];
        rep(i,lenl+1,Len)X[i]=0;
        rep(i,lenr+1,Len)Y[i]=0;
        NTT(X,1),NTT(Y,1);
        rep(i,0,Len)X[i]=(ll)X[i]*Y[i]%mod;
        NTT(X,-1);
        Rep(i,0,lenl+lenr)c[i]=X[i];
        rep(i,lenl+lenr+1,Len)c[i]=0;
    }

    inline void Inv(int*F,int*G,int ln)
    {
        Iv[0]=power(F[0],mod-2);
        for(register int Ln=2;Ln>>1<=ln;Ln<<=1)
        {
            rep(i,ln+1,Ln)F[i]=0;
            rep(i,0,Ln)X[i]=F[i],Y[i]=0;
            rep(i,0,(Ln>>1))Y[i]=Iv[i];
            Len=Ln,calrev();
            NTT(X,1),NTT(Y,1);
            rep(i,0,Ln)X[i]=(uint64)X[i]*Y[i]%mod;
            NTT(X,-1);
            rep(i,0,(Ln>>1))X[i]=0;
            NTT(X,1);
            rep(i,0,Ln)X[i]=(uint64)X[i]*Y[i]%mod;
            NTT(X,-1);
            rep(i,(Ln>>1),Ln)Iv[i]=mod-X[i];
        }
        Rep(i,0,ln)G[i]=Iv[i];
    }

    static int ExX[MAXN],ExY[MAXN],Op[MAXN];

    inline void Div(int*F,int*G,int*Q,int*R,int lenf,int leng)
    {
        Rep(i,0,lenf)ExX[i]=F[lenf-i];
        Rep(i,0,leng)ExY[i]=G[leng-i];
        Rep(i,leng+1,lenf-leng)ExY[i]=0;
        Inv(ExY,ExY,lenf-leng);
        Rep(i,lenf-leng+1,lenf)ExX[i]=0;
        Rep(i,lenf-leng+1,leng)ExY[i]=0;
        mul(ExX,ExY,ExY,lenf-leng,lenf-leng);
        Rep(i,lenf-leng+1,(lenf-leng)<<1)ExY[i]=0;
        Rep(i,0,(lenf-leng)>>1)swap(ExY[i],ExY[lenf-leng-i]);
        mul(ExY,G,ExX,lenf-leng,leng);
        assert(F[leng]==ExX[leng]);
        Rep(i,0,leng-1)R[i]=ad(F[i],mod-ExX[i]);
        Rep(i,0,lenf-leng)Q[i]=ExY[i];
    }

    const int inv2=(mod+1)>>1;

    inline void Sqrt(int*F,int*G,int ln)
    {
        Op[0]=sqrt(F[0]);
        for(register int Ln=2;Ln>>1<=ln;Ln<<=1)
        {
            mul(Op,Op,ExX,Ln>>1,Ln>>1);
            rep(i,0,Ln>>1)ExX[i]=ad(mod-ExX[i+(Ln>>1)],F[i+(Ln>>1)]);
            rep(i,Ln>>1,Ln)Op[i]=ExX[i]=0;
            Inv(Op,ExY,Ln>>1);
            Len=Ln,calrev();
            NTT(ExX,1),NTT(ExY,1);
            rep(i,0,Ln)ExX[i]=(ll)ExX[i]*ExY[i]%mod;
            NTT(ExX,-1);
            rep(i,0,Ln>>1)Op[i+(Ln>>1)]=(ll)ExX[i]*inv2%mod;
        }
        Rep(i,0,ln)G[i]=Op[i];
    }

    inline void Deriv(int*F,int*G,int ln)
    {Rep(i,1,ln)G[i-1]=(uint64)F[i]*i%mod;G[ln]=0;}

    inline void Inter(int*F,int*G,int ln)
    {Repe(i,ln,1)G[i]=(uint64)F[i-1]*iv[i]%mod;G[0]=0;}

    static int LnX[MAXN];

    inline void Ln(int*F,int*G,int ln)
    {
        Deriv(F,LnX,ln),Inv(F,G,ln);
        for(Len=2;Len<=ln<<1;Len<<=1);
        rep(i,ln+1,Len)LnX[i]=G[i]=0;
        calrev();
        NTT(LnX,1),NTT(G,1);
        rep(i,0,Len)G[i]=(uint64)G[i]*LnX[i]%mod;
        NTT(G,-1);
        Inter(G,G,ln);
    }

    inline void Exp(int *F,int *G,int ln)
    {
        Op[0]=1;
        for(register int Lx=2;Lx>>1<=ln;Lx<<=1)
        {
            rep(i,Lx>>1,Lx)Op[i]=0;
            Ln(Op,ExX,Lx);
            rep(i,0,Lx>>1)ExX[i]=ad(F[i+(Lx>>1)],mod-ExX[i+(Lx>>1)]);
            rep(i,0,Lx>>1)ExY[i]=Op[i];
            rep(i,Lx>>1,Lx)ExX[i]=ExY[i]=0;
            Len=Lx;
            calrev();
            NTT(ExY,1),NTT(ExX,1);
            rep(i,0,Len)ExX[i]=(ll)ExX[i]*ExY[i]%mod;
            NTT(ExX,-1);
            rep(i,0,Lx>>1)Op[i+(Len>>1)]=ExX[i];
        }
        Rep(i,0,ln)G[i]=Op[i];
    }
}
using poly::power;
using poly::Len;
using poly::calrev;
using poly::NTT;
using poly::mod;
using poly::predone;
using poly::Inv;
using poly::Inter;
using poly::Deriv;
using poly::Ln;
using poly::Exp;
using poly::Sqrt;
using poly::Div;
using poly::ad;

static int n,k,F[MAXN];

int main()
{
    file();
    predone();
    read(n),read(k);
    k%=mod;
    Rep(i,0,n-1)read(F[i]);
    Sqrt(F,F,n-1);
    Inv(F,F,n-1);
    Inter(F,F,n-1);
    Exp(F,F,n-1);
    Inv(F,F,n-1);
    ++F[0];
    Ln(F,F,n-1);
    ++F[0];
    Ln(F,F,n-1);
    Rep(i,0,n-1)F[i]=(ll)F[i]*k%mod;
    Exp(F,F,n-1);
    Deriv(F,F,n-1);
    Rep(i,0,n-1)write(F[i],' ');
    flush();
    return 0;
}

推荐这种东西写成一个 $namespace$ ,方便调试和复板子