多项式总结

此处多项式算法，主要针对基于快速傅里叶变换 (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(1-x)=-\sum_{i=1}^\infty \frac{x^i}{i} $$ 存在推导式: $$\ln F=\int \frac{F'}{F}\equiv \int F'F^{-1}$$ 直接计算即可。 #### 代码: ```cpp 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 ，方便调试和复板子。