多项式学习笔记

mod998244353 · 2020-05-10 18:19:29 · 个人记录

求大佬指错QaQ。

个人推荐的题单：多项式全家桶

顺便搬运下学习路线：

• 前置知识
  • 复数
  • 数论（原根）
  • 生成函数 / 形式幂级数
  • 微积分基础
  • 各种推式子方法
  • 各种卡常技巧
• 基本操作
  • 快速傅里叶变换 / FFT
  • 快速数论变换 / NTT
• 常见技巧
  • 分治 FFT
  • 多项式牛顿迭代（用于推导式子）
  • 任意模数 FFT：三模数 NTT / 拆系数 FFT

• 进阶科技

  • 多项式乘法逆 / 求逆

  前置：牛顿迭代

  • 多项式对数函数 / 求 ln

  前置：求导、积分、求逆

  • 多项式指数函数 / 求 exp

  前置：牛顿迭代、求 ln

  另：求逆 / 求 ln / 求 exp 都可以分治 FFT，并且求 exp 的话分治可能比牛迭快而且好写很多

  • 多项式幂函数 / 求幂

  前置：求 ln、求 exp

  • 多项式开根

  前置：求幂 / 牛顿迭代、求逆（、二次剩余）

  • 多项式除法 / 取模

  前置：求逆

  • 多项式多点求值

  前置：分治 FFT（、取模）

  • 多项式快速插值

  前置：多点求值、拉格朗日插值

FFT

可以用这题练习：P1919 A*B Problem升级版（FFT快速傅里叶）

前置知识：

复数是什么？

可以把复数理解成一个向量或者一个平面直角坐标系上的点。复数有一个实部和一个虚部，正如一个向量（或点）有一个横坐标和一个纵坐标。例如复数5+7i，实部是5，虚部是7，其中i=\sqrt{−1}。可以把它想象成向量(5,7)或点(5,7)。

代码实现：

struct cp {//inline不要去掉，会快很多
    double x,y;
    inline cp (double xx=0,double yy=0) {
        x=xx,y=yy;
    } inline cp operator + (const cp &a) const {
        return cp(x+a.x,y+a.y);
    } inline cp operator - (const cp &a) const {
        return cp(x-a.x,y-a.y);
    } inline cp operator * (const cp &a) const {
        return cp(x*a.x-y*a.y,x*a.y+y*a.x);
    } inline cp operator / (const cp &a) const {
        double t=a.x*a.x+a.y*a.y;
        return cp((x*a.x+y*a.y)/t,(y*a.x-x*a.y)/t);
    } inline cp operator / (const double &a) const {
        return cp(x/a,y/a);
    } inline cp operator ~ () {
        return cp(x,-y);
    }
} omega[2][MAXN];

FFT是什么？有什么用？我可以学会吗？

最后一个问题纯属凑数

这是一种在O(n\log{n})的时间内将多项式转成点值的算法（这里即下文的n为2的整数次幂，n不是时即取满足2^{k-1}<n<2^k的2^k）。

算两个n-1次多项式A(x),B(x)的积时，普通算法是O(n^2)（高精度乘法即x=10的情况）。

我们会有一个想法：两个点值表示的多项式相乘，只需要O(n)，即C_i=A_iB_i。

for(int i=0;i<n;++i)c[i]=a[i]*b[i];

所以我们算多项式乘法，可以考虑先用DFT(离散傅里叶变换)转成点值，相乘后再用IDFT转回多项式。

DFT用到的点的横坐标

这个要用到n个复数，他们不是随机找的，而是把单位圆（圆心为原点，半径长度为1）n等分，取这n个点表示的复数，即分别以这这n个点的横坐标为实部，纵坐标为虚部。我画的图太丑了，自己画下图吧。

从点(1,0)开始，逆时针从0开始将这n个点编号，标到k的点记作\omega^k_n。而且可以求出它们的值：\omega^k_n=\cos\frac{k}{n}2\pi+i\sin\frac{k}{n}2\pi（即它对应点(\cos\frac{k}{n}2\pi,\sin\frac{k}{n}2\pi)），可以看出，\omega^k_n是\omega^1_n的k次方，且\omega_n^n=1，所以\omega^1_n被称为n次单位根。

为什么要以单位根为x代入？

设(y_0,y_1,\cdots,y_{n-1})为A(x)=\sum\limits_{i=0}^{n-1}a_ix^i的DFT（点的横坐标是上述的单位根，省略不写）

再设B(x)=\sum\limits_{i=0}^{n-1}y_ix^i，把上面n个单位根的倒数\omega^0_n,\omega^{-1}_n,\cdots,\omega^{-n+1}_n代入B(x)得到DFT(z_0,z_1,\cdots,z_{n-1})，则有：

结论：把A(x)的DFT换作为另一个多项式B(x)的系数，取单位根的倒数\omega^0_n,\omega^{-1}_n,\cdots,\omega^{-n+1}_n代入B(x)，得到的每个数除以n，得到的是A(x)的各项系数。

但由于DFT仍是O(n^2)，所以有了现在的FFT。

FFT的数学证明：

设A_*(x)=\sum\limits_{i=0}^{\frac{n}{2}-1}a_{2i}x^{i},A_{**}(x)=\sum\limits_{i=0}^{\frac{n}{2}-1}a_{2i+1}x^{i}

则有：A(x)=A_*(x^2)+xA_{**}(x^2)

设k<\frac{n}{2}，把x=\omega^k_n代入：

对于A(\omega^{k+\frac{n}{2}}_n):\\A(\omega^{k+\frac{n}{2}}_n)\\=A_*(\omega^{2k+n}_n)+\omega^{k+\frac{n}{2}}_nA_{**}(\omega^{2k+n}_n)\\=A_*(\omega^{k}_{\frac{n}{2}})-\omega^k_nA_{**}(\omega^{k}_{\frac{n}{2}})

所以，如果我们知道两个多项式A_*(x)和A_{**}(x)分别在(\omega^0_{\frac{n}{2}},\omega^1_{\frac{n}{2}},\cdots,\omega^{\frac{n}{2}-1}_{\frac{n}{2}})的点值表示，就可以O(n)求出A(x)在(\omega^0_n,\omega^1_n,\cdots,\omega^{n-1}_n)处的点值表示了。而A_*(x)和A_{**}(x)都是规模缩小了一半的子问题。分治边界是n=1，此时直接return。

现在你就能写出FFT了

优化FFT

在进行FFT时，原来的位置与最终的位置时有规律的（以n=8为例子）：

一开始：\boxed{0,1,2,3,4,5,6,7}

第一次：\boxed{0,2,4,6},\boxed{1,3,5,7}

第二次：\boxed{0,4},\boxed{2,6},\boxed{1,5},\boxed{3,7}

第三次：\boxed{0},\boxed{4},\boxed{2},\boxed{6},\boxed{1},\boxed{5},\boxed{3},\boxed{7}

规律：一个数a上的数，最后位置在a二进制翻转的数，如6(110)变成了3(001)。

所以可以这样码FFT：先放到最后的位置上，然后不断向上还原，同时求点值表示。（预处理\omega^k_n,\omega_n^{-k}）

#include<bits/stdc++.h>
using namespace std;
const int MAXM=22,MAXN=1<<MAXM;
const double PI=acos(-1.0);
struct cp {
    double x,y;
    inline cp (double xx=0,double yy=0) {
        x=xx,y=yy;
    } inline cp operator + (const cp &a) const {
        return cp(x+a.x,y+a.y);
    } inline cp operator - (const cp &a) const {
        return cp(x-a.x,y-a.y);
    } inline cp operator * (const cp &a) const {
        return cp(x*a.x-y*a.y,x*a.y+y*a.x);
    } inline cp operator / (const cp &a) const {
        double t=a.x*a.x+a.y*a.y;
        return cp((x*a.x+y*a.y)/t,(y*a.x-x*a.y)/t);
    } inline cp operator / (const double &a) const {
        return cp(x/a,y/a);
    } inline cp operator ~ () {
        return cp(x,-y);
    }
} omega[2][MAXN];
const cp I(0,1);
int rev[MAXN];
void FFT(cp*a,int op,int n) {
    for(int i=0,p=n>>1; i<n; ++i) {
        rev[i]=(rev[i>>1]>>1)|(i&1?p:0);
        if(i<rev[i])swap(a[i],a[rev[i]]);
    }
    static cp t;
    for(int m=1,l=2,L=MAXM-1; l<=n; m=l,l<<=1,--L)
        for(int i=0; i<n; i+=l)
            for(int j=0; j<m; ++j) {
                t=omega[op][j<<L]*a[i|j|m];
                a[i|j|m]=a[i|j]-t,a[i|j]=a[i|j]+t;
            }
    if(!op) for(int i=0; i<n; ++i)a[i]=a[i]/n;
}

P1919AC思路

把A(x)放入a的实部，B(x)放入a的虚部，然后计算a^2，因为(x-yi)^2=(x^2-y^2)+2xyi，取虚部最后多除以2即可，就可以只FFT两次了。

P1919完整AC代码：

#include<bits/stdc++.h>
using namespace std;
const int MAXM=22,MAXN=1<<MAXM;
const double PI=acos(-1.0);
struct cp {
    double x,y;
    inline cp (double xx=0,double yy=0) {
        x=xx,y=yy;
    } inline cp operator + (const cp &a) const {
        return cp(x+a.x,y+a.y);
    } inline cp operator - (const cp &a) const {
        return cp(x-a.x,y-a.y);
    } inline cp operator * (const cp &a) const {
        return cp(x*a.x-y*a.y,x*a.y+y*a.x);
    } inline cp operator / (const cp &a) const {
        double t=a.x*a.x+a.y*a.y;
        return cp((x*a.x+y*a.y)/t,(y*a.x-x*a.y)/t);
    } inline cp operator / (const double &a) const {
        return cp(x/a,y/a);
    } inline cp operator ~ () {
        return cp(x,-y);
    }
} omega[2][MAXN];
const cp I(0,1);
int rev[MAXN];
void FFT(cp*a,int op,int n) {
    for(int i=0,p=n>>1; i<n; ++i) {
        rev[i]=(rev[i>>1]>>1)|(i&1?p:0);
        if(i<rev[i])swap(a[i],a[rev[i]]);
    }
    static cp t;
    for(int m=1,l=2,L=MAXM-1; l<=n; m=l,l<<=1,--L)
        for(int i=0; i<n; i+=l)
            for(int j=0; j<m; ++j) {
                t=omega[op][j<<L]*a[i|j|m];
                a[i|j|m]=a[i|j]-t,a[i|j]=a[i|j]+t;
            }
    if(!op) for(int i=0; i<n; ++i)a[i]=a[i]/n;
}
char str1[MAXN],str2[MAXN];
int n,m,res[MAXN];
cp f[MAXN];
int main() {
    for(int i=0; i<MAXN; ++i)
        omega[1][i]=cp(cos(PI*2*i/MAXN),sin(PI*2*i/MAXN)),omega[0][i]=~omega[1][i];
    scanf("%s%s",str1,str2);
    n=strlen(str1),m=strlen(str2);
    for(int i=0; i<n; i++)
        f[i].x=int(str1[n-1-i]^'0');
    for(int i=0; i<m; i++)
        f[i].y=int(str2[m-1-i]^'0');
    for(m+=n,n=1; n<=m; n<<=1);
    FFT(f,1,n);
    for(int i=0; i<n; ++i)f[i]=f[i]*f[i];
    FFT(f,0,n);
    for(int i=0; i<=m; ++i) {
        res[i]+=int(f[i].y/2+0.5);
        if(res[i]>9)res[i+1]+=res[i]/10,res[i]%=10;
    }
    while(!res[m])--m;
    for(int i=m; ~i; --i)
        printf("%d",res[i]);
    return 0;
}

MTT(任意模数NTT)

可以用这题练习：P4245 【模板】任意模数NTT

思路：

给出A(x),B(x)，求A(x)B(x)，系数对mod取模。

先设M=2^{16}

再设\begin{cases}A_1(x)=\sum\limits_{i=0}^{n}{\left\lfloor\frac{a_{i}}{M}\right\rfloor}x\\A_2(x)=\sum\limits_{i=0}^{n}({a_i\operatorname{mod}{M}})x\\B_1(x)=\sum\limits_{i=0}^{n}{\left\lfloor\frac{b_{i}}{M}\right\rfloor}x\\B_2(x)=\sum\limits_{i=0}^{n}({b_i\operatorname{mod}{M}})x\end{cases}

这样就有：

这样确保了不会炸精度，可我们就要做8次FFT了。

优化1：

设两个多项式A(x),B(x)且A(x),B(x)虚部都为0，且需要作DFT。

设\begin{cases}P(x)=A(x)+iB(x)\\Q(x)=A(x)-iB(x)\end{cases}

特判：(P(\omega_{n}^{0}))=\text{conj} (Q(\omega_{n}^{n}))=\text{conj} (Q(\omega_{n}^{0}))

又由于A,B虚部为0，则P,Q的每一项共轭，所以P,Q每一项的点值也共轭，所以P做完FFT时，用O(n)即可求出Q的FFT，然后再用O(n)就可以求出A,B的FFT啦：

void twice_FFT(cp *a,cp *b,int n) {
    for(int i=0; i<n; ++i)a[i]=a[i]+b[i]*I;
    FFT(a,1,n);
    b[0]=~a[0];
    for(int i=1; i<n; ++i)b[i]=~a[n-i];
    for(int i=0; i<n; ++i) {
        cp x=a[i],y=b[i];
        a[i]=(x+y)*0.5,b[i]=(y-x)*0.5*I;
    }
}

用优化1，就可以用两次FFT求出A_1,A_2,B_1,B_2的FFT了。

优化2

乘法运算之后，要对四个多项式A_1B_1,A_2B_1,A_1B_2,A_2B_2作FFT，可以这样：

设\begin{cases}P=A_1B_1\\Q=A_2B_1+A_1B_2+iA_2B_2\end{cases}

由于A_1,A_2,B_1,B_2虚部为0，那么A_1B_1,A_2B_1,A_1B_2,A_2B_2卷起来的虚部一定为0，那么此时P的实部和虚部分别为A_1B_1,0，Q的实部和虚部分别为A_2B_1+A_1B_2,A_2B_2。

最后答案就是P的实部\times2^{32}+Q的实部\times2^{16}+Q的虚部

void MTT(int*a,int*b,int m,int*ans,const int mod) {
    static cp f[MAXN],g[MAXN],p[MAXN],q[MAXN],t,f0,f1,g0,g1;
    int n=1;
    for(; n<m; n<<=1);
    for(int i=0; i<m; ++i) {
        f[i]=cp(a[i]>>16,a[i]&65535);
        g[i]=cp(b[i]>>16,b[i]&65535);
    }
    FFT(f,1,n),FFT(g,1,n);
    for(int i=0; i<n; ++i) {
        t=~f[i?n-i:0],f0=(t-f[i])*0.5*I,f1=(t+f[i])*0.5;
        t=~g[i?n-i:0],g0=(t-g[i])*0.5*I,g1=(t+g[i])*0.5;
        p[i]=f1*g1,q[i]=f1*g0+f0*g1+f0*g0*I;
    }
    FFT(p,0,n),FFT(q,0,n);
    for(int i=0; i<n; ++i)
        ans[i]=(((ll)(p[i].x+0.5)%mod<<32)%mod+((ll)(q[i].x+0.5)%mod<<16)%mod+(ll)(q[i].y+0.5))%mod;
    clr(f,n),clr(g,n),clr(p,n),clr(q,n);
}

然后是P4245的AC代码：

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
const int MAXM=18,MAXN=1<<MAXM;
const double PI=acos(-1.0);
struct cp {
    double x,y;
    cp (double xx=0,double yy=0) {
        x=xx,y=yy;
    } cp operator + (const cp &a) const {
        return cp(x+a.x,y+a.y);
    } cp operator - (const cp &a) const {
        return cp(x-a.x,y-a.y);
    } cp operator * (const cp &a) const {
        return cp(x*a.x-y*a.y,x*a.y+y*a.x);
    } cp operator / (const cp &a) const {
        double t=a.x*a.x+a.y*a.y;
        return cp((x*a.x+y*a.y)/t,(y*a.x-x*a.y)/t);
    } cp operator / (const double &a) const {
        return cp(x/a,y/a);
    } cp operator ~ () {
        return cp(x,-y);
    }
} omega[2][MAXN];
const cp I(0,1);
int rev[MAXN];
void FFT(cp*a,int op,int n) {
    for(int i=0,p=n>>1; i<n; ++i) {
        rev[i]=(rev[i>>1]>>1)|(i&1?p:0);
        if(i<rev[i])swap(a[i],a[rev[i]]);
    }
    static cp t;
    for(int m=1,l=2,L=MAXM-1; l<=n; m=l,l<<=1,--L)
        for(int i=0; i<n; i+=l)
            for(int j=0; j<m; ++j) {
                t=omega[op][j<<L]*a[i|j|m];
                a[i|j|m]=a[i|j]-t,a[i|j]=a[i|j]+t;
            }
    if(!op) for(int i=0; i<n; ++i)a[i]=a[i]/n;
}
inline void clr(int*f,int n) {
    memset(f,0,n<<2);
}
inline void cpy(int*f,int*g,int n) {
    memcpy(f,g,n<<2);
}
inline void clr(cp*f,int n) {
    memset(f,0,n<<4);
}
inline void cpy(cp*f,cp*g,int n) {
    memcpy(f,g,n<<4);
}
void MTT(int*a,int*b,int m,int*ans,const int mod) {
    static cp f[MAXN],g[MAXN],p[MAXN],q[MAXN],t,f0,f1,g0,g1;
    int n=1;
    for(; n<m; n<<=1);
    for(int i=0; i<m; ++i) {
        f[i]=cp(a[i]>>16,a[i]&65535);
        g[i]=cp(b[i]>>16,b[i]&65535);
    }
    FFT(f,1,n),FFT(g,1,n);
    for(int i=0; i<n; ++i) {
        t=~f[i?n-i:0],f0=(t-f[i])*0.5*I,f1=(t+f[i])*0.5;
        t=~g[i?n-i:0],g0=(t-g[i])*0.5*I,g1=(t+g[i])*0.5;
        p[i]=f1*g1,q[i]=f1*g0+f0*g1+f0*g0*I;
    }
    FFT(p,0,n),FFT(q,0,n);
    for(int i=0; i<n; ++i)
        ans[i]=(((ll)(p[i].x+0.5)%mod<<32)%mod+((ll)(q[i].x+0.5)%mod<<16)%mod+(ll)(q[i].y+0.5))%mod;
    clr(f,n),clr(g,n),clr(p,n),clr(q,n);
}
int n,m,mod,a[MAXN],b[MAXN],ans[MAXN];
cp f[MAXN];
int main() {
    for(int i=0; i<MAXN; ++i)
        omega[1][i]=cp(cos(PI*2*i/MAXN),sin(PI*2*i/MAXN)),omega[0][i]=~omega[1][i];
    scanf("%d%d%d",&n,&m,&mod);
    for(int i=0; i<=n; ++i)scanf("%d",a+i);
    for(int i=0; i<=m; ++i)scanf("%d",b+i);
    MTT(a,b,n+m,ans,mod);
    for(int i=0; i<=n+m; i++)printf("%d ",ans[i]);
    return 0;
}

NTT

可以用这题练习：P3803 【模板】多项式乘法（FFT）

思路：

在FFT中，我们需要用到复数，可复数需要用double类型计算，精度太低，那有没有什么东西能够代替复数且解决精度问题呢？

这个东西，叫原根

原根

阶

若\gcd(a,p)=1，且p>1，对于a^n\equiv\operatorname{(mod}p)最小的n，我们称之为a模p的阶，记作\delta_p(a)

原根

定义：

设p是正整数，a是整数，若\delta_p(a)等于\phi(p)，则称a为模p的一个原根。（注意原根的个数是不唯一的，如果模数p有原根，那么它一定有\phi(\phi(p))个原根）

如：\delta_7(3)=\phi(p)=6，有3为7的原根。

一个非常重要的定理：

若P\in\operatorname{prime}，p为P的原根，则p^i\mod P的结果两两不同。

为什么原根能代替单位根进行运算

在FFT中，我们用到了单位根的四条性质，而原根也满足这四条性质，这样我们最终可以得到：

再把FFT中的\omega_n全部替换掉即可。

除个F(x)：G(x)-G_*(x)\equiv0\operatorname{(mod}{x^{\left\lceil\frac{n}{2}\right\rceil}})

平方并展开：G(x)^2-2G(x)G_*(x)+G_*(x)^2\equiv0\operatorname{(mod}{x^n})

乘个F(x)：F(x)G(x)^2-2F(x)G(x)G_*(x)+F(x)G_*(x)^2\equiv0\operatorname{(mod}{x^n})

代入F(x)G(x)\equiv1\operatorname{(mod}{x^n})后化简可得：

所以就可以这样递归啦，边界：当多项式F(x)只有一项，那么显然G_0是F_0的逆元。

时间复杂度：T(n)=T(\frac{n}{2})+O(n\log n)=O(n\log n)

代码（配合前面基本操作食用）：

void poly_inv(int*f,int m) {
    static int g[MAXN],sav[MAXN],n;
    g[0]=fastpow(f[0],mod-2);
    for(n=1; n<m; n<<=1);
    for(int len=2; len<=n; len<<=1) {
        cpy(sav,f,len);
        NTT(sav,1,len<<1),NTT(g,1,len<<1);
        for(int i=0; i<len<<1; ++i) g[i]=(2ll+mod-sav[i]*1ll*g[i]%mod)*g[i]%mod;
        NTT(g,0,len<<1),clr(g+len,len);
    }
    cpy(f,g,m),clr(sav,n<<1),clr(g,n<<1);
}

学了MTT（上面有）的同学可以试试这题：P4239 任意模数多项式乘法逆。

多项式半在线卷积

可以用这题来练习：P4721 【模板】分治 FFT

这题其实就是要我们用强制中序遍历转移的cdq分治。我们考虑每次计算跨越隔板的贡献。（以样例一为例，与结果无关的值用?表示，不与G相加）

一开始G=(1,0,0,0)

隔板：G(\boxed{1|0},0,0)

与F的第一部分(0,3)做卷积，得到(?,3)，与G相加：

隔板：G(\boxed{1,3|0,0})

与F的第二部分(0,3,1,2)做卷积，得到(?,?,10,5)，与G相加：

隔板：G(1,3,\boxed{10|5})

与F的第一部分(0,3)做卷积，得到(?,30)，与G相加：

这就是答案啦！

时间复杂度：T(n)=2T(\frac{n}{2})+O(n\log n)=O(n\log^2n)

代码实现：

int s1[MAXN],s2[MAXN];
void cdq_NTT(int*F,int*G,int l,int r) {
    if (l+1==r)return ;
    int mid=(l+r)>>1,n=r-l,m=n>>1;
    cdq_NTT(F,G,l,mid);
    cpy(s1,G+l,m);
    clr(s1+m,m);
    NTT(s1,1,n);
    cpy(s2,F,n);
    NTT(s2,1,n);
    px(s1,s2,n);
    NTT(s1,0,n);
    for (int i=m; i!=n; ++i) {
        G[l+i]+=s1[i];
        if(G[l+i]>=mod)G[l+i]-=mod;
    }
    cdq_NTT(F,G,mid,r);
}
void cdq_times(int*f,int*g,int m) {
    int n;
    for(g[0]=n=1; n<m; n<<=1);
    cdq_NTT(f,g,0,n);
}

多项式ln

多项式求导和积分

小建议：最好先自行了解简单微积分。

定义多项式的求导：\dfrac{d}{dx}F(x)表示F(x)求导后的结果（简写为F^{\prime}(x)）。

然后我们根据求导四则运算之一(u(x)+v(x))^{\prime}=u^{\prime}(x)+v^{\prime}(x)和单项式求导的一般形式(Cx^a)^{\prime}=Cax^{a-1}可以得到：

然后多项式不定积分（求导的逆运算，符号为\int\!F(x)dx）也就可以求了。

代码实现：

void poly_dao(int*f,int n) {
    for(int i=1; i<n; ++i)
        f[i-1]=(0ll+f[i])*i%mod;
    f[n-1]=0;
}
void poly_jifen(int*f,int n) {
    for(int i=n; i; --i)f[i]=(f[i-1]+0ll)*inv[i]%mod;
    f[0]=0;
}

多项式ln

可以用这题练习：【模板】多项式对数函数（多项式 ln）

用复合函数求导法则（\dfrac{dG(F(x))}{dx}=G^{\prime}(F(x))F^{\prime}(x)）对式子求导，得到\ln^{\prime}(A(x))A^{\prime}(x)=B^{\prime}(x)

再把\ln^{\prime}(A(x))=A(x)^{-1}代入得：

两边再同时求不定积分：

时间复杂度：O(n\log n)(同求逆)

代码：（是不是很短）

void poly_ln(int*f,int n) {
    static int g[MAXN];
    cpy(g,f,n),poly_inv(g,n),poly_dao(f,n);
    times(f,g,n,n),clr(g,n<<1);
    poly_jifen(f,n-1);
}

多项式exp

可以用这题练习：P4726 【模板】多项式指数函数（多项式 exp）

泰勒展开

将F(x)在x_0处展开：F(x)=\sum\limits_{i=0}^{\infty}\frac{F^{(i)}(x_0)}{i!}(x-x_0)^i（F^{(n)}(x)表示F(x)的n阶导数）

多项式牛顿迭代

当G(F(x))\equiv0\pmod{x^{n}},F(x)\equiv\!F_*(x)\pmod{x^{\lceil\frac{n}{2}\rceil}}时，F(x)\equiv\!F_*(x)-\frac{G(F_*(x))}{G^{\prime}(F_*(x))}（F(x)为\pmod{x^n}意义下的解，F_*(x)为\pmod{x^{\lceil\frac{n}{2}\rceil}}意义下的解，例子：下面的\exp）

证明：（不想看就直接看exp吧）

我们考虑把G(F(x))在F_*(x)处泰勒展开：

将G(F(x))\equiv0\pmod{x^n}代入上式后化简可得F(x)\equiv\!F_*(x)-\frac{G(F_*(x))}{G^{\prime}(F_*(x))}\pmod{x^n}

多项式exp

根据e^{A(x)}\equiv\!B(x)\pmod{x^n}可得\ln(B(x))-A(x)\equiv0\pmod{x^n}，我们就设G(B(x))={\ln(B(x))-A(x)}。

将这个式子的B(x)看作主元，A(x)看作常数，求导后可得G^{\prime}(B(x))=B^{-1}(x)

将G(B_*(x))={\ln(B_*(x))-A(x)}和G^{\prime}(B_*(x))=B_*^{-1}(x)代入B(x)\equiv\!B_*(x)-\frac{G(B_*(x))}{G^{\prime}(B_*(x))}\pmod{x^n}可得

也就是B(x)\equiv(1-\ln(B_*(x))+A(x))B_*(x)\pmod{x^n}，靠这个倍增即可。

当A只有一项时，B_0=e^{A_0}=e^0=1。

时间复杂度：T(n)=T(\frac{n}{2})+O(n\log n)=O(n\log n)

代码如下：

void poly_exp(int*f,int m) {
    static int g[MAXN],h[MAXN],n;
    for(assert(!f[0]),n=1; n<m; n<<=1);
    g[0]=1;
    for(int len=2; len<=n; len<<=1) {
        cpy(h,g,n),poly_ln(h,len);
        for(int i=1; i<len; ++i) h[i]=(f[i]-h[i]+mod)%mod;
        h[0]=(1+mod+f[0]-h[0])%mod;
        times(g,h,len,len);
    }
    cpy(f,g,m),clr(g,n<<1),clr(h,n<<1);
}

半在线版多项式exp

提取x^{n+1}项：(n+1)B_{n+1}=\sum\limits_{i=0}^n(i+1)A_{i+1}B_{n-i}

化简：B_n=\frac{1}{n}\sum\limits_{i=1}^niA_iB_{n-i}，用分治NTT即可，时间复杂度O(n\log^2n)，这玩意常数极小，跑得比上面的牛顿迭代还快（\frac{1}{n}在分治边界乘上、i在分治前乘上即可）

void cdq_exp(int*F,int*G,int l,int r) {
    if (l+1==r)
        if(l>0)G[l]=1ll*fastpow(l,mod-2)*G[l]%mod;
        else ;
    else {
        int mid=(l+r)>>1,n=r-l,m=n>>1;
        cdq_exp(F,G,l,mid);
        cpy(s1,G+l,m);
        clr(s1+m,m);
        NTT(s1,1,n);
        cpy(s2,F,n);
        NTT(s2,1,n);
        px(s1,s2,n);
        NTT(s1,0,n);
        for (int i=m; i!=n; ++i) {
            G[l+i]+=s1[i];
            if(G[l+i]>=mod)G[l+i]-=mod;
        }
        cdq_exp(F,G,mid,r);
    }
}
int SSSSS[MAXN];
void poly_cdq_exp(int*F,int m) {
    for(int i=0; i!=m; ++i)SSSSS[i]=F[i]*1ll*i%mod;
    clr(F,m);
    int n=1;
    for(F[0]=1; n<m; n<<=1);
    cdq_exp(SSSSS,F,0,n);
    clr(SSSSS,n);
    clr(F+m,n-m);
}

多项式快速幂及加强版

可以做下这两道题：

• P5245 【模板】多项式快速幂
• P5273 【模板】多项式幂函数 (加强版)

普通版

void poly_pow(int*f,int m,int k) {
    poly_ln(f,m);
    for(int i=0; i<m; ++i)f[i]=(0ll+f[i])*k%mod;
    poly_exp(f,m);
}
void poly_cdq_pow(int*f,int m,int k) {
    poly_ln(f,m);
    for(int i=0; i<m; ++i)f[i]=(0ll+f[i])*k%mod;
    poly_cdq_exp(f,m);
}

加强版

设B_*(x)满足B_*^2(x)\equiv A(x)\pmod{x^{\lceil\frac{n}{2}\rceil}}，则B(x)-B_*(x)\equiv0\pmod{x^{\lceil\frac{n}{2}\rceil}}

平方并展开：B^2(x)+B_*^2(x)\equiv2B(x)B_*(x)\pmod{x^n}

把B^2(x)\equiv A(x)\pmod{x^n}代入后化简可得：

递归边界：n=1时，b_0=\sqrt{a_0}=1

int TTT[MAXN],TTTT[MAXN];
void poly_sqrt(int*f,int m) {
    int n=1;
    for(; n<m; n<<=1);
    TTT[0]=1;
    for(int len=2; len<=n; len<<=1) {
        for(int i=0,p=len>>1; i<p; ++i) {
            TTTT[i]=TTT[i]<<1;
            if(TTTT[i]>=mod)TTTT[i]-=mod;
        }
        poly_inv(TTTT,len);
        NTT(TTT,1,len);
        px(TTT,TTT,len);
        NTT(TTT,0,len);
        for (int i=0; i<len; ++i) {
            TTT[i]+=f[i];
            if(TTT[i]>=mod)TTT[i]-=mod;
        }
        times(TTT,TTTT,len,len);
    }
    cpy(f,TTT,m);
    clr(TTT,n<<1);
    clr(TTTT,n<<1);
}

加强版

就是还要求模意义开根，P5491 【模板】二次剩余。

以下讲二次剩余：

解的数量

严格来讲，非0数n是二次剩余当且仅当方程x^2 \equiv n\pmod{p}，上述方程无解的非0数n称作非二次剩余。

假设有多组解，对于其中不同的两个解x_1,x_2有x_1^2\equiv x_2^2\pmod{p}即(x_1-x_2)(x_1+x_2)\equiv0\pmod{p}

又因为x_1\not=x_2，所以x_1+x_2\equiv0\pmod{p}也就是说两个不相等的解一定是相反数，换言之，该方程只有两个解，且它们互为相反数。而当p为奇素数时模意义的两个相反数不会相等，因为奇偶性不同。

欧拉准则

如何快速判断一个正整数n是否为模p意义下的二次剩余？

以下设q=\frac{p-1}{2}

费马小定理n^{p-1}\equiv1\pmod{p}可以化成n^{2q}\equiv1\pmod{p}

也就是说n^q\equiv1或-1

若n是二次剩余，则n^q\equiv(x^2)^q\equiv x^{2q}\equiv1\pmod{p}

设n^q\equiv1\pmod{p}，将n表示为g^k（g为模p意义下的原根），那么有g^{kq}\equiv1\pmod{p}

又因为g为模p意义下的原根，所以p-1|kq即2|k，所以n为二次剩余且一个模意义开根为g^{\frac{k}{2}}。

因为n^q\mod p不是1就是-1，所以当n^q\equiv-1\pmod{p}时n不为二次剩余。

Cipolla算法

对于方程x^2\equiv n\pmod{p}

先随机找到一个a使a^2-n不是模p意义下的二次剩余。

接下来定义i^2\equiv a^2-n\pmod{p}

但是a^2 - n不是二次剩余，怎么找得到这样一个 i？

类比实数域到复数域的推广，定义这样一个i ，然后可以将所有数表示为a+bi 的形式， 其中a,b都是模p意义下的整数，类似于实部和虚部。

可以得到(a+i)^{p+1}\equiv n\pmod{p}

\color{blue}\text{证明：}

那么(a+i)^{\frac{p+1}{2}}即是一个解，其相反数是另一个解。

设存在(a+bi)^2\equiv n\pmod{p}且b\not=0

展开：a^2+b^2i^2-n\equiv-2abi\pmod{p}

将i^2\equiv a^2-n\pmod{p}代入：

式子的左边“虚部”为0，那么式子右边的虚部也一定为0，也就是说 ab\equiv0\pmod{p} ，又因为b\not=0，所以a=0。

因为b^{-2}和n都为模p意义下的二次剩余，所以nb^{-2}为模p意义下的二次剩余，与i^2不是模p意义下的二次剩余矛盾。

int squareI;
struct cp {
    int real,imag;
    cp(int real=0,int imag=0):real(real),imag(imag) {}
    inline bool operator==(cp y) {
        return real==y.real&&imag==y.imag;
    }
    inline cp operator*(cp y) {
        return cp(((0ll+real)*y.real+(0ll+squareI)*imag%mod*y.imag)%mod,((imag+0ll)*y.real+(0ll+real)*y.imag)%mod);
    }
};
cp power(cp x,int k) {
    cp res=1;
    for(; k; x=x*x,k>>=1)
        if(k&1)
            res=res*x;
    return res;
}
inline int modsqrt(int n) {
    n%=mod;
    if(n==0)return 0;
    srand(time(0));
    int a=rand()%mod;
    for(squareI=((0ll+a)*a%mod-n+mod)%mod; !a||power(squareI,mod>>1)==1; squareI=((0ll+a)*a%mod-n+mod)%mod)
        a=rand()%mod;
    int ans=power(cp(a,1),inv2).real;
    return min(ans,mod-ans);
}
int TTT[MAXN],TTTT[MAXN];
void poly_sqrt(int*f,int m) {
    int n=1;
    for(; n<m; n<<=1);
    TTT[0]=f[0]==1?1:modsqrt(f[0]);
    for(int len=2; len<=n; len<<=1) {
        for(int i=0,p=len>>1; i<p; ++i) {
            TTTT[i]=TTT[i]<<1;
            if(TTTT[i]>=mod)TTTT[i]-=mod;
        }
        poly_inv(TTTT,len);
        NTT(TTT,1,len);
        px(TTT,TTT,len);
        NTT(TTT,0,len);
        for (int i=0; i<len; ++i) {
            TTT[i]+=f[i];
            if(TTT[i]>=mod)TTT[i]-=mod;
        }
        times(TTT,TTTT,len,len);
    }
    cpy(f,TTT,m);
    clr(TTT,n<<1);
    clr(TTTT,n<<1);
}

多项式带余除法

可以用这题练习：P4512 【模板】多项式除法

设F_r(x)=x^nF(x^{-1})，则它表示F(x)系数翻转后的结果，如：

化简题目给的式子：R(x)补位到m-1次。

这样就能用求逆求Q_r(x)，翻转得Q(x)再带回题目的式子求R(x)即可：

void poly_div(int *f,int *g,int n,int m) { //F=F/G,G=F%G
    static int q[MAXN],t[MAXN];
    int L=n-m+1;
    fan(g,m),cpy(q,g,L),fan(g,m);
    fan(f,n),cpy(t,f,L),fan(f,n);
    poly_inv(q,L),times(q,t,L,L);
    fan(q,L);
    savtimes(q,g,L,m,m-1);
    for (int i=0; i<m-1; ++i) {
        f[i]-=q[i];
        if(f[i]<0)f[i]+=mod;
    }
    clr(f+m-1,L),clr(q,n),clr(t,n);
}

多项式三角函数

可以用这题练习：P5264 多项式三角函数

欧拉公式：\cos(x)+i\sin(x)=e^{ix}

也就是\begin{cases}\cos(F(x))=\dfrac{\exp(iF(x))+\exp(-iF(x))}{2}\\\\\sin(F(x))=\dfrac{\exp(iF(x))-\exp(-iF(x))}{2i}\end{cases}

其中i\equiv\sqrt{-1}\equiv\sqrt{\operatorname{mod}-1}

然后\tan(F(x))=\dfrac{\sin(x)}{\cos(x)}（虽然题目没求但写一下也是可以的）

const int I=modsqrt(mod-1),inv2I=fastpow(I<<1,mod-2);
int HHH[MAXN],JJJ[MAXN];
void poly_cos(int*f,int m) {
    for(int i=0; i<m; ++i)f[i]=f[i]*1ll*I%mod;
    poly_exp(f,m);
    cpy(HHH,f,m);
    poly_inv(HHH,m);
    for(int i=0; i<m; ++i)f[i]=(f[i]+HHH[i])*1ll*inv2%mod;
    clr(HHH,m);
}
void poly_sin(int*f,int m) {
    for(int i=0; i<m; ++i)f[i]=f[i]*1ll*I%mod;
    poly_exp(f,m);
    cpy(HHH,f,m);
    poly_inv(HHH,m);
    for(int i=0; i<m; ++i)f[i]=(f[i]+mod-HHH[i])*1ll*inv2I%mod;
    clr(HHH,m);
}
void poly_tan(int*f,int n) {
    cpy(JJJ,f,n);
    poly_sin(f,n);
    poly_cos(JJJ,n);
    poly_inv(JJJ,n);
    times(f,JJJ,n,n);
    clr(JJJ,n);
}

多项式反三角函数

可以用这题练习：P5265 多项式反三角函数

直接摆公式：

代入F(x)：（注意复合函数求导法则）

还有\arccos F(x)=-\arcsin F(x)写一下也没关系的吧（

void poly_asin(int*f,int n) {
    static int g[MAXN];
    cpy(g,f,n);
    int m=1;
    while(m<n+n)m<<=1;
    NTT(g,1,m),px(g,g,m),NTT(g,0,m),clr(g+n,m-n);
    for(int i=0; i<n; ++i)g[i]=g[i]?mod-g[i]:0;
    ++g[0];
    poly_sqrt(g,n);
    poly_inv(g,n);
    poly_dao(f,n);
    times(f,g,n,n);
    poly_jifen(f,n);
    clr(g,n);
}
void poly_acos(int*f,int n) {
    poly_asin(f,n);
    for(int i=0; i<n; ++i)f[i]=f[i]?mod-f[i]:0;
}
void poly_atan(int*f,int n) {
    static int g[MAXN];
    cpy(g,f,n);
    int m=1;
    while(m<n+n)m<<=1;
    NTT(g,1,m),px(g,g,m),NTT(g,0,m),clr(g+n,m-n);
    ++g[0];
    poly_inv(g,n);
    poly_dao(f,n);
    times(f,g,n,n);
    poly_jifen(f,n);
    clr(g,n);
}

拉格朗日插值法

可以用这题练习：P4781 【模板】拉格朗日插值法

作用：

在算法竞赛中，我们常常会碰到一类题目，题目中给出了n+1个点，让我们求由这些点构成的多项式在某一位置的取值。

拉格朗日插值法可以在O(n^2)的复杂度内完美解决这个问题，比高斯消元求出多项式的系数的O(n^3)复杂度优且不存在精度问题。

公式：

如果我们把x_i带入的话，第i项的可以化简成y_i，其余每项的分子中都会有x_i-x_i，这样它们就都被消去了。因此拉格朗日插值法的正确性是可以保证的。

直接套式子即可ACP4781

#include <bits/stdc++.h>
using namespace std;
const int mod=998244353;
int ans,x[2222],y[2222],n,k;
inline int read() {
    int x=0,c=getchar();
    for(; c<=47||c>=58; c=getchar());
    for(; c>=48&&c<=57; c=getchar())x=(x<<3)+(x<<1)+(c&15);
    return x;
}
int fastpow(int a,int k) {
    int base=1;
    for(; k; a=((0ll+a)*a)%mod,k>>=1)
        if(k&1)
            base=((0ll+base)*a)%mod;
    return base%mod;
}
int Lagrange(int n,int xi,int*x,int*y) {
    int ans=0;
    for (int i=1; i<=n; ++i) {
        int s1=1,s2=1;
        for(int j=1; j<=n; ++j)
            if (i!=j) {
                s1=1ll*s1*(xi-x[j])%mod;
                s2=1ll*s2*(x[i]-x[j])%mod;
                s1+=s1>>31&mod;
                s2+=s2>>31&mod;
            }
        ans+=1ll*y[i]*s1%mod*fastpow(s2,mod-2)%mod;
        if(ans>=mod)ans-=mod;
    }
    return ans<0?ans+mod:ans;
}
int main() {
    n=read(),k=read();
    for(int i=1; i<=n; ++i)x[i]=read(),y[i]=read();
    printf("%d\n",Lagrange(n,k,x,y));
    return 0;
}

多项式多点求值

可以用这题练习：P5050 【模板】多项式多点求值

我们设G_i(x)=x-a_i，那么我们就要求F(x)\mod G_i(x)

可是O(n^2\log n)比暴力还慢，那我们就考虑分治：

假设当前分治区间为[l,r]，设H_0(x)=\prod\limits_{i=l}^{mid}G_i(x)，H_1(x)=\prod\limits_{i=mid}^{r}G_i(x)。

而且目前我们已经算出F(x)\mod \prod\limits_{i=l}^{r}G_i(x)，那我们往下递归就行了。

代码：（减小常数要在递归边界写暴力，我用了秦九昭展开）

int gl[MAXM+1][MAXN];
void qfpre(int lev,int l,int r,int *x) {
    if(l==r) {
        gl[lev][l<<1]=mod-x[l];
        gl[lev][l<<1|1]=1;
        return ;
    }
    const int mid=(l+r)>>1;
    qfpre(lev+1,l,mid,x);
    qfpre(lev+1,mid+1,r,x);
    cpy(&gl[lev][l<<1],&gl[lev+1][l<<1],mid-l+2);
    savtimes(&gl[lev][l<<1],&gl[lev+1][(mid+1)<<1],mid-l+2,r-mid+1,r-l+2);
}
void queryf(int lev,int l,int r,int *f,int *x,int *y) {
    static int s1[MAXN],s2[MAXN],b[33]= {1};
    if (r<=l+2048) {
        const int L1=l<<1,L2=L1+31,R=r+l-1,F=f[r+l];
        int xx,j,now;
        for (int i=l; i<=r; ++i) {
            xx=x[i],now=F;
            for(int k=1; k!=33; ++k)b[k]=b[k-1]*1llu*xx%mod;
            for(j=R; j>=L2; j-=32) {
                now=((
                         1ull*now    *b[32]+1ull*f[j   ]*b[31]+1ull*f[j- 1]*b[30]+1ull*f[j- 2]*b[29]+
                         1ull*f[j- 3]*b[28]+1ull*f[j- 4]*b[27]+1ull*f[j- 5]*b[26]+1ull*f[j- 6]*b[25]+
                         1ull*f[j- 7]*b[24]+1ull*f[j- 8]*b[23]+1ull*f[j- 9]*b[22]+1ull*f[j-10]*b[21]+
                         1ull*f[j-11]*b[20]+1ull*f[j-12]*b[19]+1ull*f[j-13]*b[18]+1ull*f[j-14]*b[17]
                     )%mod+
                     1ull*f[j-15]*b[16]+1ull*f[j-16]*b[15]+1ull*f[j-17]*b[14]+1ull*f[j-18]*b[13]+
                     1ull*f[j-19]*b[12]+1ull*f[j-20]*b[11]+1ull*f[j-21]*b[10]+1ull*f[j-22]*b[ 9]+
                     1ull*f[j-23]*b[ 8]+1ull*f[j-24]*b[ 7]+1ull*f[j-25]*b[ 6]+1ull*f[j-26]*b[ 5]+
                     1ull*f[j-27]*b[ 4]+1ull*f[j-28]*b[ 3]+1ull*f[j-29]*b[ 2]+1ull*f[j-30]*b[ 1]+f[j-31]
                    )%mod;
            }
            for(; j>=L1; --j)now=(1ull*now*xx+f[j])%mod;
            y[i]=now;
        }
        return ;
    }
    const int mid=(l+r)>>1,siz=r-l+1;
    cpy(s1,&f[l<<1],siz),poly_div(s1,&gl[lev+1][l<<1],siz,mid-l+2);
    cpy(s2,&f[l<<1],siz),poly_div(s2,&gl[lev+1][(mid+1)<<1],siz,r-mid+1);
    clr(f+(l<<1),(r-l)<<1|1);
    cpy(&f[l<<1],s1,mid-l+1),cpy(&f[(mid+1)<<1],s2,r-mid);
    clr(s1,siz),clr(s2,siz);
    queryf(lev+1,l,mid,f,x,y),queryf(lev+1,mid+1,r,f,x,y);
}
void poly_queryf(int *f,int *x,int *y,int n,int m) {
    qfpre(0,0,m-1,x);
    if(n>m)poly_div(f,gl[0],n,m+1);
    queryf(0,0,m-1,f,x,y);
}

多项式快速插值

可以用这题练习：P5158 【模板】多项式快速插值

以下内容自由元X一律大写：

拉格朗日插值：

1. 计算每个\prod\limits_{i=1,i\not=k}^n(x_k-x_i)

设G(X)=\prod\limits_{i=1}^n(X-x_i)

那么我们就是要求\lim_{X\to x_k}\dfrac{G(x)}{X-x_k}=\prod\limits_{i=1,i\not=k}^n(x_k-x_i)

这个式子就是\dfrac{0}{0}，由洛必达法则得到：

然后用多项式多点求值就可以求这些数了

再设V_k=\dfrac{y_k}{G^{\prime}(x_k)}，则现在有：

可以考虑分治：

设F_{[l,r]}(X)=\sum\limits_{k=l}^r[V_k\prod\limits_{i=l,k\not=i}^r(X-x_i)]，mid=\dfrac{l+r}{2}

边界：F_{k,k}=V_k，而且多点求值的G与分治NTT的G是一个东西，不用算两次。

代码：

void polf(int lev,int l,int r,int*f,int*x,int*v) {
    if (l==r) {
        f[l<<1]=v[l];
        return ;
    }
    const int mid=(l+r)>>1,siz=r-l+1,nlev=lev+1;
    const int L=l<<1;
    polf(nlev,l,mid,f,x,v),polf(nlev,mid+1,r,f,x,v);
    cpy(s1,&f[L],mid-l+1);
    savtimes(s1,&gl[nlev][(mid+1)<<1],mid-l+1,r-mid+1,siz);
    cpy(s2,&f[(mid+1)<<1],r-mid);
    savtimes(s2,&gl[nlev][L],r-mid,mid-l+2,siz);
    clr(f+L,r-l<<1);
    for (int i=0; i<siz; ++i)
        f[i+L]=s1[i]+s2[i]>=mod?s1[i]+s2[i]-mod:s1[i]+s2[i];
    clr(s1,siz),clr(s2,siz);
}
void poly_polf(int m,int*f,int*x,int*y) {
    static int g[MAXN],y2[MAXN];
    qfpre(0,0,m-1,x);
    cpy(g,gl[0],m+1),poly_dao(g,m+1);
    queryf(0,0,m-1,g,x,y2);
    for(int i=0; i<m; i++)y2[i]=y[i]*1ll*getinv(y2[i])%mod;
    polf(0,0,m-1,f,x,y2);
    clr(g,m),clr(y2,m);
}

多项式复合函数

可以用这题练习：P5373 【模板】多项式复合函数

令L=\lceil\sqrt{n}\rceil

预处理G^i(x),G^{iL}(x)，复杂度O(n\sqrt{n}\log n)

由拉格朗日反演得：

预处理(\dfrac{x}{F(x)})^{i},(\dfrac{x}{F(x)})^{iL}，复杂度O(n\sqrt{n}\log n)

我们设f(d,x)=\prod\limits_{i=1}^dx+i，

则n!=\prod\limits_{i=0}^{s-1}f(s,is)\times\prod\limits_{i=s^2+1}^ni（这里最多要乘3s个数）。

这样我们只要求出f(s,0),f(s,s),\cdots,f(s(s-1),s)的值即可。

初始条件有f(1,0)=1。

倍增要有以下两个操作：

1. 知f(d,0),f(d,s),\cdots,f(d,ds)求f(d+1,0),f(d+1,s),\cdots,f(d+1,(d+1)s)
2. 知f(d,0),f(d,s),\cdots,f(d,ds)求f(2d,0),f(2d,s),\cdots,f(2d,2ds)

把d加1

首先有f(2d,x)=f(d,x)f(d,d+x)

那我们就可以考虑求f(d,0),f(d,s),\cdots,f(d,2ds)和f(d,d),f(d,d+s),\cdots,f(2ds+d)

我们设g(x)=f(d,sx)

我们考虑用g(0),g(1),g(d)求g(d+1+0),g(d+1+1),\cdots,g(d+1+d)

用g(0),g(1),\cdots,g(2d)求g(d/s+0),g(d/s+1),\cdots,g(d/s+2d)。

我们就可以写一个函数void calc(int delta,int cur,int*ip,int*op)

来用g(0),g(1),\cdots,g(\operatorname{cur})求g(\Delta+0),g(\Delta+1),\cdots,g(\Delta+\operatorname{cur})

之后乘2的方法如下：

calc(cur+1,cur,val,fv1);//计算g(cur+1)~g(2d+1)
cpy(&val[cur+1],fv1,cur+1),val[cur<<1|1]=0;//g(2d+1)不要
calc(cur*1ll*getinv(n)%mod,cur<<1,val,fv2);//计算g(d/s)~g(d/s+2d)
cur<<=1;//上一行的n实际是s，看完整代码就知道了
for(int i=0; i<=cur; ++i)val[i]=val[i]*1ll*fv2[i]%mod;

为了写那个函数，我们尝试使用拉格朗日插值：

设F(i)=\dfrac{h(i)}{(-1)^{\operatorname{cur}-i}i!(\operatorname{cur}-i)!},G(i)=\dfrac{1}{\Delta+x-i},H(x)=\dfrac{h(\Delta+x-\operatorname{cur})}{\prod\limits_{j=0}^{\operatorname{cur}}(\Delta+x-j)}

就有H(x)=\sum\limits_{i+j=\operatorname{cur},i,j\geq0}F(i)G(j)，用MTT卷下就可以推出H(x+\operatorname{cur})。

知道了H(x+\operatorname{cur})之后我们就可以推出h(x)来了,不难发现H(n+k)和h(n)的比值是一段区间当中数字的乘积，这个系数在可以2-pointer扫一遍在O(\operatorname{cur}+\log mod)的时间内求出。

这样一次是O(s\log s)的。

P5282完整AC代码：（时间复杂度O(\sqrt{n}\log n)）

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
inline int read() {
    int x=0,c=getchar();
    for(; c<=47||c>=58; c=getchar());
    for(; c>=48&&c<=57; c=getchar())x=(x<<3)+(x<<1)+(c&15);
    return x;
}
const int MAXM=18,MAXN=1<<MAXM;
struct cp {
    cp (double xx=0,double yy=0) {
        x=xx,y=yy;
    }
    double x,y;
    cp operator + (cp const &B) const {
        return cp(x+B.x,y+B.y);
    }
    cp operator - (cp const &B) const {
        return cp(x-B.x,y-B.y);
    }
    cp operator * (cp const &B) const {
        return cp(x*B.x-y*B.y,x*B.y+y*B.x);
    }
    cp operator / (cp const &B) const {
        double t=B.x*B.x+B.y*B.y;
        return cp((x*B.x+y*B.y)/t,(y*B.x-x*B.y)/t);
    }
    cp operator / (double const &B) const {
        return cp(x/B,y/B);
    }
    cp operator ~ () {
        return cp(x,-y);
    }
} omega[2][MAXN];
const double PI=acos(-1.0);
const cp I(0,1),II(0,-0.5);
int tf,tr[MAXN];
void tpre(const int limit) {
    if(tf==limit)return ;
    tf=limit;
    for(int i=0,p=limit>>1; i<limit; ++i)
        if(i&1)tr[i]=(tr[i>>1]>>1)|p;
        else tr[i]=tr[i>>1]>>1;
}
void FFT(cp*f,const bool op,const int n) {
    tpre(n);
    for(int i=0; i<n; ++i)
        if(i<tr[i])
            swap(f[i],f[tr[i]]);
    for(int len=1,p=2,L=1; p<=n; len=p,p<<=1,++L) {
        for(int i=0; i<n; i+=p) {
            for(int j=0; j<len; j++) {
                int u=i|j,v=u|len;
                cp tt=omega[op][j<<MAXM-L]*f[v];
                f[v]=f[u]-tt;
                f[u]=f[u]+tt;
            }
        }
    }
    if(!op) {
        for(int i=0; i<n; ++i)
            f[i]=f[i]/n;
    }
}
void clr(int*f,int n) {
    memset(f,0,n<<2);
}
void cpy(int*f,int*g,int n) {
    memcpy(f,g,n<<2);
}
void clr(cp*f,int n) {
    memset(f,0,n<<4);
}
void cpy(cp*f,cp*g,int n) {
    memcpy(f,g,n<<4);
}
int mod;
int inv[MAXN],ifac[MAXN];
int fastpow(int a,int k) {
    int base=1;
    for(; k; a=((0ll+a)*a)%mod,k>>=1)
        if(k&1)
            base=((0ll+base)*a)%mod;
    return base%mod;
}
int getinv(int a) {
    return a<MAXN?inv[a]:fastpow(a,mod-2);
}
void getpre() {
    const int p=min(MAXN,mod);
    inv[0]=inv[1]=ifac[0]=ifac[1]=1;
    for(int i=2; i<p; ++i) {
        inv[i]=inv[mod%i]*(mod+0ull-mod/i)%mod;
        ifac[i]=ifac[i-1]*1ull*inv[i]%mod;
    }
}
void MTT(const int*a,const int*b,const int m,int *ans) {
    static cp f[MAXN],g[MAXN],p[MAXN],q[MAXN];
    int nn=1;
    for(; nn<m; nn<<=1);
    const int len=nn;
    for(int i=0; i<m; ++i) {
        f[i]=cp(a[i]>>16,a[i]&65535);
        g[i]=cp(b[i]>>16,b[i]&65535);
    }
    FFT(f,1,len),FFT(g,1,len);
    cp t=~f[0],f0=(f[0]-t)*II,f1=(f[0]+t)*0.5;
    t=~g[0];
    cp g0=(g[0]-t)*II,g1=(g[0]+t)*0.5;
    p[0]=f1*g1,q[0]=f1*g0+f0*g1+f0*g0*I;
    for(int i=1; i<len; ++i) {
        t=~f[len-i],f0=(f[i]-t)*II,f1=(f[i]+t)*0.5;
        t=~g[len-i],g0=(g[i]-t)*II,g1=(g[i]+t)*0.5;
        p[i]=f1*g1,q[i]=f1*g0+f0*g1+f0*g0*I;
    }
    FFT(p,0,len),FFT(q,0,len);
    for(int i=0; i<len; ++i)
        ans[i]=((((ll)(p[i].x+0.5)%mod<<16)%mod+(ll)(q[i].x+0.5)<<16)%mod+((ll)(q[i].y+0.5))%mod)%mod;
    clr(f,len),clr(g,len),clr(p,len),clr(q,len);
}
namespace GetFac {
    int f[MAXN],g[MAXN],h[MAXN];
    inline void calc(int delta,int cur,int*ip,int*op) {
        const int cur2=cur+cur,cur3=cur2+cur;
        int len=1;
        while(len<=cur3)len<<=1;
        for(int i=0; i<=cur; ++i)f[i]=ip[i]*1ll*ifac[i]%mod*(ll)ifac[cur-i]%mod;
        for(int i=cur-1; i>=0; i-=2)f[i]=mod-f[i];
        for(int i=0; i<=cur2; ++i)g[i]=getinv((0ll+delta+mod+i-cur)%mod);
        clr(&f[cur+1],len-cur-1),clr(&g[cur+cur+1],len-cur2-1);
        MTT(f,g,len,h);
        int p1=delta-cur,p2=delta,xs=1;
        for(int i=p1; i<=p2; ++i)xs=(xs*1ll*i)%mod;
        for(int i=0; i<=cur; ++i) {
            op[i]=h[i+cur]*1ll*xs%mod;
            ++p2;
            xs=xs*1ll*getinv(p1)%mod*p2%mod;
            ++p1;
        }
    }
    int val[MAXN],fv1[MAXN],fv2[MAXN];
    inline void solve(int n) {
        const int invn=getinv(n);
        int hb=0;
        for(int p=n; p; p>>=1)++hb;
        val[0]=1;
        for(int z=hb,cur=0; z>=0; --z) {
            if(cur!=0) {
                calc(cur+1,cur,val,fv1);
                cpy(&val[cur+1],fv1,cur+1),val[cur<<1|1]=0;
                calc(cur*1ll*invn%mod,cur<<1,val,fv2);
                cur<<=1;
                for(int i=0; i<=cur; ++i)val[i]=val[i]*1ll*fv2[i]%mod;
            }
            if((n>>z)&1) {
                for(int i=0; i<=cur; ++i)val[i]=(1ll*n*i+cur+1ll)%mod*val[i]%mod;
                cur|=1,val[cur]=1;
                for(int i=1; i<=cur; ++i)val[cur]=(1ll*n*cur+i)%mod*val[cur]%mod;
            }
        }
    }
    inline int GetF(int n) {
        if(n<=1024) {
            int jc=1;
            for(int i=2; i<=n; ++i)
                jc=jc*1ll*i%mod;
            return jc;
        }
        int bl=sqrt(n),res=1,i=0,id=0;
        for(solve(bl); ; i+=bl,++id) {
            if(i+bl>n) {
                for(++i; i<=n; ++i)res=res*1ll*i%mod;
                return res;
            }
            res=res*1ll*val[id]%mod;
        }
    }
    inline int Get(int n) {
        if(n>mod-(n+1))return n&1?getinv(GetF(mod-1-n)):mod-getinv(GetF(mod-1-n));
        return GetF(n);
    }
}
int n,m,a[MAXN],b[MAXN],ans[MAXN];
int main() {
    for(int i=0; i<MAXN; ++i)
        omega[1][i]=cp(cos(PI*2*i/MAXN),sin(PI*2*i/MAXN)),omega[0][i]=~omega[1][i];
    int t=read();
    while(t--) {
        int n=read();
        mod=read();
        getpre();
        printf("%d\n",GetFac::Get(n));
    }
    return 0;
}

update:

准备写的内容：

1. 下降幂多项式乘法
2. 下降幂多项式转普通多项式
3. 普通多项式转下降幂多项式
4. 斯特林数