lswzk

· · 个人记录

除法、开根、exp、多点插值、快速求值、拉格朗日插值的板子还没有整理，重心拉格朗日插值法还不会，正在补坑中。

\color{red}{\text{约定：}}

$2.w$代表单位根，$w_m$表示$m$次单位根。 $3.A$代表一个数列。 $4.g$表示原根。 $5.F^{(z)}(x)$表示$F(x)$的$z$次导数。 ### $\color{red}{\text{多项式系列之零 底层知识：}}

多项式的表示：

多项式可以通过系数数列A表示，a_i是x_i的系数。

多项式可以通过点值表示，对于一个n次多项式，取n种不同的x取值带入F(x)，得到n个值，在取相同这n个数的意义下，可以唯一的表示这个多项式。

多项式乘法：

定义F(x)\oplus G(x)=\sum_{i=0}^n\sum_{j=0}^i f_ix^i\times g_{j-i}x^{j-i}，在系数表示之下相乘复杂度\Theta(n^2)，在点值表示之下F(x)\oplus G(x)=A_f\times A_g=\sum_{i=1}^n a_{fi}\times a_{gi}，复杂度\Theta(n)。

复数：

复数一般情况下可以表示成a+bi的形式，a,b是实数，i=\sqrt{-1}。

复数的幅角：平面直角坐标系上点(a,b)所在的任意角。

复数的模长：\sqrt{a^2+b^2}

两个复数相乘：(a+bi)\times(c+di)=ac+adi+bci-bd=(ac-bd)+(ad+bc)i，复数相乘之后，模长等于原来两个复数的模长的乘积，幅角的角度等于原来两个幅角的和。

复数可以加减乘除，可以和实数一样的带入F(x)。

单位根：

在单位圆上从w_m^0=(1,0)开始平均取m个点，从0开始编号，分别是w_m^0,w_m^1,w_m^2,w_m^3\cdots w_m^{m-1}。

画图观察可得：

$2.w_m^{-k}=(cos(\frac {-k} n 2\pi),sin(\frac {-k} n 2\pi))$所代表的复数 $3.w_m^m\;=w_m^0=(1,0) 4.w_m^m\;=-w_m^{\frac m 2} 5.w_{2m}^{2k}\;=w_m^k 6.w_m^{k+\frac m 2}\!=-w_m^k

DFT&IDFT：

在科学的数学函数意义上DFT是讲一个函数转化成三角函数的加减乘除的形式，三角函数的系数是原函数系数与点值之间的变换规律。IDFT是DFT的逆变换。

原根：

（NTT用）

$2.g$存在的条件：$p=2,4,q^a,2q^a$，$q$是奇素数。 $3.$如何求$g$：把$\varphi(p)$进行质因数分解$\varphi(p)=\prod p_i^{a_i}$，如果对于任意的$p_i$，总有$g^{\frac {\varphi(p)} {p_i}}\neq 1(mod~p)$，暴力枚举即可。 #### **CRT合并：** （NTT用） 求解${\begin{cases}x\equiv a_1 (mod~p_1)\\x\equiv a_2 (mod~p_2)\end{cases}}

由x\equiv a_1 (mod~p_1)，得

\begin{aligned}x-p_1y&=a_1\\x&=a_1+p_1y\end{aligned}

带入二式，得

\begin{aligned}a_1+p_1y&\equiv a_2(mod~p_2)\\p_1y&\equiv a_2-a_1(mod~p_2)\end{aligned}

若\gcd(p1,p2)==1，用逆元直接除便可；否则通过exgcd可求得y，若无解则方程组无解。

最后x=p_1y+a_1(mod~p_1p_2)。

泰勒展开：

（牛顿迭代用）

简单来说，F(x)在x_0处展开就是

F(x)=\sum_{i=0}^{\infty}\frac{F^{(i)}(x_0)}{i!}(x-x_0)^i \color{red}{\text{咕}}\color{green}{\text{咕}}\color{red}{\text{咕}}\color{green}{\text{咕}}\color{red}{\text{咕}}\color{green}{\text{咕}}\color{red}{\text{咕}}\color{green}{\text{咕}}\color{red}{\text{咕}}\color{green}{\text{咕}}\color{red}{\text{咕}}\color{green}{\text{咕}}\color{red}{\text{咕}}\color{green}{\text{咕}}\color{red}{\text{咕}}\color{green}{\text{咕}}

牛顿迭代：

（开根、exp用）

牛顿迭代是用来求解零点问题的一种方法。

已知G(z) (这是一个函数，如G(z)=z^2-H(x)那么我们就是对H(x)开根) 求G(F(x))\equiv 0(mod\;x^n)的F(x)。我们尝试递归求解。

当$n==1$时，$G(F(x))\equiv0(mod\;x)$直接求解即可，与$G(x)$的具体形式有关。 $2.$**转移**： 假设我们已经得知$G(J(x))\equiv0(mod\;x^{\lceil\frac n 2\rceil})$： 我们运用泰勒展开把$G(F(x))$在$J(x)$处展开，得 $$G(F(x))\equiv\sum_i\frac{G^{(i)}(J(x))}{i!}(F(x)-J(x))^i\qquad(mod\;x^n)$$ 由于$F(x)$和$J(x)$的前$\lceil\frac n 2\rceil$相同，所以$F(x)-J(x)$的最低非零项大于$\lceil\frac n 2\rceil$，所以$(F(x)-J(x))^2$的最低非零项大于$2\lceil\frac n 2\rceil$取模后为$0$。 $$\begin{aligned}G(F(x))&\equiv G(J(x))+G^{(1)}(J(x))(F(x)-J(x))&(mod\;x^n)\\F(x)&\equiv\frac{G(F(x))-G(J(x))}{G^{(1)}(J(x))}+J(x)&(mod\;x^n)\end{aligned}$$ 由于$G(F(x))\equiv 0\;(mod\;x^n): F(x)\equiv-\frac{G(J(x))}{G^{(1)}(J(x))}+J(x)\quad(mod\;x^n)

拉格朗日插值法：

（多项式快速插值用）

给定点集$X=\{(x_0,y_0),(x_1,y_1)\cdots\}$，求一个函数$F(x)$，使得$\forall x_i,F(x_i)=y_i 构造函数$L_i(x)=\prod_{i\ne j}\frac{x-x_j}{x_i-x_j}$。 当$x=x_i$时，$L_i(x)=1$。 当$x\in X,x\ne x_i$时，$L_i(x)=0

这样，F(x)=\sum_{i}L_i(x)\times y_i

复杂度\Theta(n^2)

在朴素拉格朗日插值法中，如果我们新插入一个节点，修改的复杂度是$\Theta(n^2)$的，因为每个函数都要修改，而函数的修改要先减去它再修改再加回来，加减的复杂度是$\Theta(n)$的。 我们推一下式子： $$\begin{aligned}F(x)&=\sum_{i}L_i(x)\times y_i\\&=\sum_{i}y_i\prod_{i\ne j}\frac{x-x_j}{x_i-x_j}\\&=\sum_{i}\frac{\prod_{j}x-x_j}{x-x_i}\prod_{i\ne j}\frac{y_i}{x_i-x_j}\\\end{aligned}$$ 令$w_i=\prod_{i\ne j}\frac{1}{x_i-x_j},l(x)=\prod_jx-x_j \begin{aligned}F(x)&=\sum_i\frac{l(x)w_i}{x-x_i}y_i\\&=l(x)\sum_i\frac{w_i}{x-x_i}y_i\end{aligned}

修改时，我们只需要

\color{red}{\text{多项式全集之一 FFT:}}

什么是FFT：

FFT是利用DFT的特殊性质，把w带入x从而\Theta(nlogn)求一个系数多项式的点值表示，所以叫FDFT。

w的具体应用：

设$F(x)$的系数是$A$，在$w_m$的DFT下点值是$B$，$G(x)$的系数是$B$，在$w_m^{-1}$的DFT下点值是$C$。 $$\begin{aligned}c_k&=\sum_{i=0}^{m-1}b_iw_m^{-ki}\\&=\sum_{i=0}^{m-1}(\sum_{j=0}^{m-1}a_jw_m^{ji})w_m^{-ki}\\&=\sum_{j=0}^{m-1}a_j\sum_{i=0}^{m-1}w_m^{(j-k)i}\end{aligned}$$ 当$j-k=0$时$\sum_{i=0}^{m-1}w_m^{(j-k)i}=m$，否则根据等比数列求和公式得 $$\frac{w_m^{(j-k)m}-1}{w_m^{j-k}-1}=\frac{w_m^{m(j-k)}-1}{w_m^{j-k}-1}=\frac{1-1}{w_m^{j-k}-1}=0$$ 由此可得：$c_k=m\times a_k$，$a_k=\frac{c_k}{m}

综上所述，对于点值取的w相反数做DFT再除以m可得到系数。

直接将$w$带入多项式做DFT需要复杂度$\Theta(n^2)$，我们利用$w$的性质优化： 把$F(x)$按照奇偶分裂，$F(x)=(a_0x^0+a_2x^2+a_4x^4+a_6x^6\cdots)+(a_1x^1+a_3x^3+a_5x^5+a_7x^7\cdots)

我们令

\begin{aligned}F_0(x)&=a_0x^0+a_2x^1+a_4x^2+a_6x^3\cdots+a_{m-2}x^{\frac m 2}\\F_1(x)&=a_1x^0+a_3x^1+a_5x^2+a_6x^3\cdots+a_{m-1}x^{\frac m 2-1}\end{aligned}

我们可以发现F(x)=F_0(x^2)+xF_1(x^2)

现在我们把w_m带入，令k<\frac m 2：

\begin{aligned}F(w_m^k)&=F_0(w_m^{2k})+w_m^kF_1(w_m^{2k})\\&=F_0(w_{\frac m 2}^{k})+w_m^kF_1(w_{\frac m 2}^{k})\\F(w_m^{k+\frac m 2})&=F_0(w_m^{2k+m})+w_m^{k+\frac m 2}F_1(w_m^{2k+m})\\&=F_0(w_{\frac m 2}^{k} )-w_m^{k}F_1(w_{\frac m 2}^{k})\end{aligned}

我们知道取w_{\frac m 2}时，A_0,A_1的取值，就可以算出A(w_m)，而A_0,A_1的长度都为A的一半，所以可以递归计算。

非递归优化FFT：

画图可知，递归版FFT最底层结束状态第$i$个位置的项是$i$二进制翻转后的结果。我们可以$\Theta(n)$的得到最底层的结果，然后向上模拟回溯合并即可。 $2.$蝴蝶变换： 由上述式子： $$\begin{aligned}F(w_m^k)&=F_0(w_{\frac m 2}^{k})+w_m^kF_1(w_{\frac m 2}^{k})\\F(w_m^{k+\frac m 2})&=F_0(w_{\frac m 2}^{k} )-w_m^{k}F_1(w_{\frac m 2}^{k})\end{aligned}$$ 可得在迭代时$w_m^k$与$w_m^{k+\frac m 2}$都只与$w_{\frac m 2}^k,w_{\frac m 2}^{k+\frac m 2}$有关，所以我们可以用临时变量记录下一层的两个信息向上迭代。 #### **共轭复数优化FFT：** $1.$优化原理： （在DFT时） 令 $$\begin{aligned}D(x)=A(x)+iB(x)\\E(x)=A(x)-iB(x)\end{aligned}$$ 那么 $$\begin{aligned}E_D(x)&=conj(D_D(m-x))\\A_D(x)&=\frac{D_D(x)+E_D(x)} 2 \\&=\frac{D_D(x)+conj(D_D(m-x))}2\\B_D(x)&=\frac{D_D(x)-E_D(x)} {2i}\\&=-i\frac{D_D(x)-conj(D_D(m-x))} {2}\end{aligned}$$ $2.$证明： $step~1$： $$\begin{aligned}D_D(w_m^k)&=A_D(w_m^k)+iB_D(w_m^k)\\D_D(w_m^k)&=\sum_{j=0}^{m-1}a_jw_m^{jk}+ib_jw_m^{jk}\\D_D(w_m^k)&=\sum_{j=0}^{m-1}(a_j+ib_j)w_m^{jk}\end{aligned}$$ $step~2$：为方便起见，我们用$X$代替$\frac{2\pi j k} {m} \begin{aligned}E_D(w_m^k)&=A_D(w_m^k)-iB_D(w_m^k)\\&=\sum_{j=0}^{m-1}a_jw_m^{jk}-ib_jw_m^{jk}\\&=\sum_{j=0}^{m-1}(a_j-ib_j)w_m^{jk}\\&=\sum_{j=0}^{m-1}(a_j-ib_j)(cosX+isinX)\\&=\sum_{j=0}^{m-1}(a_jcosX+b_jsinX)+i(a_jsinX-b_jcosX)\\&=conj\big(\sum_{j=0}^{m-1}(a_jcosX+b_jsinX)-i(a_jsinX-b_jcosX)\big)\\&=conj\big(\sum_{j=0}^{m-1}(a_jcosX+b_jsinX)-i(a_jsinX-b_jcosX)\big)\\&=conj\big(\sum_{j=0}^{m-1}(a_jcos(-X)-b_jsin(-X))+i(a_jsin(-X)+b_jcos(-X))\big)\\&=conj\big(\sum_{j=0}^{m-1}(a_j+ib_j)(cos(-X)+isin(-X)))\big)\\&=conj\big(\sum_{j=0}^{m-1}(a_j+ib_j)w_m^{-jk}\big)\\&=conj\big(\sum_{j=0}^{m-1}(a_j+ib_j)w_m^{jm-jk}\big)=conj(D_D(w_m^{m-k}))\\\end{aligned}

而在IDFT时，我们需要

\begin{aligned}D_D(w_m^{-k})&=\sum_{j=0}^{m-1}(a_j+ib_j)w_m^{-jk}\\E_D(w_m^{-k})&=conj\big(\sum_{j=0}^{m-1}(a_j+ib_j)w_m^{jk}\big)=conj(D_D(w_m^{k}))\end{aligned}

数论优化FFT(NTT)：

$1.(g^{\frac {p-1} m})^k$和$w_m^k$都互不相同$(k\in[0,m-1])$。 $2.(g^{\frac {p-1} m})^m=g^{p-1}=1$，$w_m^m=1$。 $3.(g^{\frac {p-1} m})^{\frac m 2}=g^{\frac{p-1} 2}=\sqrt{g^{p-1}}$，由于原根的互不相同，$(g^{\frac {p-1} m})^{\frac m 2}=-1=-g^{p-1}=-(g^{\frac {p-1} m})^m$，$w_m^m=-w_m^{\frac m 2} 4.(g^{\frac {p-1} {2m}})^{2k}=(g^{\frac {p-1} {m}})^{k}$，$w_{2m}^{2k}=w_m^k 5.(g^{\frac {p-1} {m}})^{k+\frac m 2}=(g^{\frac {p-1} {m}})^{k}\times (g^{\frac {p-1} {m}})^{\frac m 2}=-(g^{\frac {p-1} {m}})^{k}$，$w_m^{k+\frac m 2}=-w_m^k

因为有这些共性，所以g^{\frac {p-1} m}可以代替w_m。

喜闻乐见的模板：

FFT模板（共轭优化）

namespace FFT{
    const double pi = acos(-1);
    struct cp{
        double x, y;
        cp() {x = y = 0;}
        cp(double X,double Y) {x = X; y = Y; }
        cp conj() {return (cp) {x, -y};}
    }a[3000005], b[3000005], c[3000005], I(0, 1), d[3000005];

    cp operator+ (const cp &a, const cp &b) {return (cp){a.x + b.x, a.y + b.y}; }
    cp operator- (const cp &a, const cp &b) {return (cp){a.x - b.x, a.y - b.y}; }
    cp operator* (const cp &a, const cp &b) {return (cp){a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x};}
    cp operator* (const cp &a, double b) {return (cp){a.x * b, a.y * b};}
    cp operator/ (const cp &a, double b) {return (cp){a.x / b, a.y / b};}
    struct p_l_e{
        int wz[3000005];
        void FFT(cp *a, int N, int op) {
            for(int i = 0; i < N; i++)
                if (i<wz[i]) swap(a[i],a[wz[i]]);
            for(int le = 2; le <= N; le <<= 1) {
                int mid = le >> 1;
                for(int i = 0; i < N; i += le) {
                    cp x, y, w = (cp) {1, 0};
                    cp wn = (cp){cos(op * pi / mid), sin(op * pi / mid)};
                    for(int j = 0 ; j < mid; j++) {
                        x = a[i + j]; y = a[i + j + mid] * w;
                        a[i + j] = x + y;
                        a[i + j + mid] = x - y;
                        w = w * wn;
                    }
                }
            }
        }
        void D_FFT(cp *a, cp *b, int N, int op){
            for(int i = 0; i < N; i++)  d[i] = a[i] + I * b[i];
            FFT(d, N, op);
            d[N] = d[0];
            for(int i = 0; i < N; i++){
                a[i] = (d[i] + d[N - i].conj()) / 2;
                b[i] = I * (-1) * (d[i] - d[N - i].conj()) / 2;
            }
            d[N] = cp(0, 0);
        }
        void mult(cp *a, cp *b, cp *c, int M){
            int N = 1, len = 0;
            while(N < M) N <<= 1, len++;
            for(int i = 0; i < N; i++)
                wz[i] = (wz[i >> 1] >> 1) | ((i & 1) << (len - 1));
            D_FFT(a, b, N, 1);
            for(int i = 0; i < N; i++)  c[i] = a[i] * b[i];
            FFT(c, N, -1);
            for(int i = 0; i < N; i++)  c[i].x = c[i].x / N;
        }
    }PLE;
    int n, m;
    void main() {
        scanf("%d%d", &n, &m); n++; m++;
        for(int i = 0; i < n; i++) scanf("%lf", &a[i].x);
        for(int i = 0; i < m; i++) scanf("%lf", &b[i].x);
        PLE.mult(a, b, c, n + m - 1); 
        for(int i = 0; i < n + m - 1; i++)  printf("%d ", (int)round(c[i].x));
        return;
    }
}

NTT模板：

namespace NTT{
    typedef long long LL;
    const int mod = 998244353;
    const int g = 3;
    LL a[3000005], b[3000005], c[3000005];
    int n, m;
    LL qpow(LL a, LL b){
        LL ans = 1;
        while(b){
            if(b & 1)   ans = ans * a % mod;
            a = a * a % mod;
            b >>= 1;
        }
        return ans;
    }
    struct p_l_e{
        int wz[3000005];
        void NTT(LL *a, int N, int op) {
            for(int i = 0; i < N; i++) 
                if(i < wz[i]) swap(a[i], a[wz[i]]);
            for(int le = 2; le <= N; le <<= 1) {
                int mid = le >> 1;
                LL wn = qpow(g, (mod - 1) / le);
                if(op == -1) wn = qpow(wn, mod - 2); 
                for(int i = 0; i < N; i += le) {
                    int w = 1, x, y;
                    for(int j = 0; j < mid; j++) {
                        x = a[i + j]; 
                        y = a[i + j + mid] * w % mod; 
                        a[i + j] = (x + y) % mod;
                        a[i + j + mid] = (x - y + mod) % mod;
                        w = w * wn % mod;
                    }
                }
            }
        }
        void mult(LL *a, LL *b, LL *c, int M) {
            int N = 1, len = 0;
            while(N < M)    N <<= 1, len++;
            for(int i = 0; i < N; i++)  
                wz[i] = (wz[i >> 1] >> 1) | ((i & 1) << (len - 1));
            NTT(a, N, 1); NTT(b, N, 1);
            for(int i = 0; i < N; i++)  c[i] = a[i] * b[i] % mod;
            NTT(c, N, -1);
            LL t = qpow(N, mod - 2);
            for(int i = 0; i < N; i++)  c[i] = c[i] * t % mod;
        }
    }PLE;
    void main() {
        scanf("%d%d", &n, &m); n++; m++;
        for(int i = 0; i < n; i++)  scanf("%lld", &a[i]);
        for(int i = 0; i < m; i++)  scanf("%lld", &b[i]);
        PLE.mult(a, b, c, n + m - 1);
        for(int i = 0; i < n + m - 1; i++)  printf("%lld ", c[i]);
    }
}

\color{red}{\text{多项式全集之二 任长任模的FFT:}}

三模NTT现任意模数FFT（MTT）：

$2.$MTT的使用原理：我们对初始多项式取模，那么如果在不取模卷积情况下，答案$x$不会超过$N\times p^2$。我们取三个NTT模数$p_1,p_2,p_3$，分别做多项式乘法，得到$x$分别$mod~p_1,p_2,p_3$的答案，通过CRT合并可以得到$x~mod~p_1p_2p_3$的答案，如果$x<p_1p_2p_3$那么就可以得到准确的答案，再对$p$取模即可。 $3.$CRT合并的小优化： $step~0:$初始式子 $${\begin{cases}x\equiv c_1(mod~p_1)\\x\equiv c_2(mod~p_2)\\x\equiv c_3(mod~p_3)\end{cases}}$$ $step~1:$把一式二式合并(LL范围内)。 $${\begin{cases}x\equiv a(mod~p_1p_2)\\x\equiv c_3(mod~p_3)\end{cases}}$$ $step~2:$再次合并(不需要$long~double$ 快速乘)。 $4.$常用NTT模数： 以下模数的共同$g=3189 p=r\times 2^k+1$|$k$|$g

:-:|:-:|:-:

104857601$|$22$|$3 167772161$|$25$|$3 469762049$|$26$|$3 950009857$|$21$|$7 998244353$|$23$|$3 1004535809$|$21$|$3 2013265921$|$27$|$31 2281701377$|$27$|$3 3221225473$|$30$|$5

拆系数FFT(CFFT)实现任模FFT：

$2.$可按合并DFT的方法优化DFT次数。 #### **$bluestein$算法实现任意长长度FFT：** 可兼容NTT 当$m$不是$2$的幂次的时候，我们从式子入手： 令$X_i=a_iw_m^{\frac {i^2} 2},Y_i=w_m^{\frac{-i^2}2} \begin{aligned}A_k & = \sum_{j=0}^{m-1}a_jw_m^{jk}\\ & = \sum_{j=0}^{m-1}a_jw_m^{\frac{j^2+k^2-{(k-j)}^2}{2}}\\&=w_m^{\frac {k^2} 2}\sum_{j=0}^{m-1}a_jw_m^{\frac{j^2} 2}w_m^{\frac{-{(k-j)}^2}{2}}\\&=w_m^{\frac {k^2} 2}\sum_{j=0}^{m-1}X_jY_{k-j}\end{aligned}

喜闻乐见的模板：

三模NTT模板（注意：不可以MTT回来，因为系数会取模）

namespace MTT{
    typedef long long LL;
    int n, m;
    LL p, mod;
    const LL p1 = 998244353;
    const LL p2 = 1004535809;
    const LL p3 = 104857601;
    const int g = 3189;
    LL a[300005], b[300005], c[300005], cpa[300005], cpb[300005];
    LL c3[300005], c1[300005], c2[300005];
    LL qpow(LL a, LL b, LL mod) {
        LL ans = 1;
        while(b) {
            if(b & 1)   ans = ans * a % mod;
            a = a * a % mod;
            b >>= 1;
        }
        return ans;
    }
    const LL inv12 = qpow(p1, p2 - 2, p2);
    const LL inv123 = qpow(p1 * p2 % p3, p3 - 2, p3);
    struct p_l_e{
        int wz[300005];
        void MTT(LL *a, int N, int op) {
            for(int i = 0; i < N; i++)
                if(i < wz[i]) swap(a[i], a[wz[i]]);
            for(int le = 2; le <= N; le <<= 1) {
                int mid = le >> 1;
                LL wn = qpow(g, (mod - 1) / le, mod);
                if(op == -1) wn = qpow(wn, mod - 2, mod);
                for(int i = 0; i < N ;i += le) {
                    LL w = 1, x, y;
                    for(int j = 0; j < mid; j++) {
                        x = a[i + j];
                        y = a[i + j + mid] * w % mod;
                        a[i + j] = (x + y) % mod;
                        a[i + j + mid] = (x - y + mod) % mod;
                        w = w * wn % mod;
                    }
                }
            }
        }
        void mult(LL *a, LL *b, LL *c, int M) {
            int N = 1, len = 0;
            while(N < M) N <<= 1, len++;
            for(int i = 0; i < N; i++)
                wz[i] = (wz[i >> 1] >> 1) | ((i & 1) << (len - 1));
            MTT(a, N, 1); MTT(b, N, 1);
            for(int i = 0; i < N; i++)  c[i] = a[i] * b[i] % mod;
            MTT(c, N, -1);
            LL t = qpow(N, mod - 2, mod);
            for(int i = 0; i < N; i++)  c[i] = c[i] * t % mod;
        }

    }PLE;
    LL CRT(LL c1, LL c2, LL c3) {
        LL x = (c1 + p1 * ((c2 - c1 + p2) % p2 * inv12 % p2));
        LL y = (x % p + p1 * p2 % p * ((c3 - x % p3 + p3) % p3 * inv123 % p3) % p) % p;
        return y;
    }
    void merge(LL *c1, LL *c2, LL *c3, LL *c, int N) {
        for(int i = 0; i < N; i++)
            c[i] = CRT(c1[i], c2[i], c3[i]);
        return;
    }
    void main() {
        scanf("%d%d%lld", &n, &m, &p); n++; m++;
        for(int i = 0; i < n; i++)  scanf("%lld", &a[i]);
        for(int i = 0; i < m; i++)  scanf("%lld", &b[i]);
        mod = p1; memcpy(cpa, a, sizeof(a)); memcpy(cpb, b, sizeof(b)); PLE.mult(cpa, cpb, c1, n + m - 1);
        mod = p2; memcpy(cpa, a, sizeof(a)); memcpy(cpb, b, sizeof(b)); PLE.mult(cpa, cpb, c2, n + m - 1);
        mod = p3; memcpy(cpa, a, sizeof(a)); memcpy(cpb, b, sizeof(b)); PLE.mult(cpa, cpb, c3, n + m - 1);
        merge(c1, c2, c3, c, n + m - 1);
        for(int i = 0; i < n + m - 1; i++)  printf("%lld ", (c[i] % p + p) % p);
        return;
    }
}

拆系数FFT模板（注意：相同系数的两项可以合并一起IDFT。采用共轭优化法，只进行四次DFT）

namespace CFFT{
    typedef long long LL;
    int n, m, p ,sqrp; 
    int a[300005], b[300005];
    const long double pi = acos(-1);
    struct cp{
        long double x, y;
        cp() {x = y = 0;}
        cp(long double X,long double Y) {x = X; y = Y; }
        cp conj() {return (cp) {x, -y};}
    }ka[300005], kb[300005], ta[300005], tb[300005], kk[300005], kt[300005], tt[300005], c[300005], I(0, 1), d[300005];

    cp operator+ (const cp &a, const cp &b) {return (cp){a.x + b.x, a.y + b.y}; }
    cp operator- (const cp &a, const cp &b) {return (cp){a.x - b.x, a.y - b.y}; }
    cp operator* (const cp &a, const cp &b) {return (cp){a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x};}
    cp operator* (const cp &a, long double b) {return (cp){a.x * b, a.y * b};}
    cp operator/ (const cp &a, long double b) {return (cp){a.x / b, a.y / b};}  
    struct p_l_e{
        int wz[300005];
        void FFT(cp *a, int N, int op){
            for(int i = 0; i < N; i++)
                if(i < wz[i])   swap(a[i], a[wz[i]]);
            for(int le = 2; le <= N; le <<= 1){
                int mid = le >> 1;
                cp x, y, w, wn = (cp){cos(op * 2 * pi / le), sin(op * 2 * pi / le)};
                for(int i = 0; i < N; i += le){
                    w = (cp){1, 0};
                    for(int j = 0; j < mid; j++){
                        x = a[i + j];
                        y = a[i + j + mid] * w;
                        a[i + j] = x + y;
                        a[i + j + mid] = x - y;
                        w = w * wn;
                    }
                }
            } 
        }
        void D_FFT(cp *a, cp *b, int N, int op){
            for(int i = 0; i < N; i++)  d[i] = a[i] + I * b[i];
            FFT(d, N, op);
            d[N] = d[0];
            if(op == 1){
                for(int i = 0; i < N; i++){
                    a[i] = (d[i] + d[N - i].conj()) / 2;
                    b[i] = I * (-1) * (d[i] - d[N - i].conj()) / 2;
                }
            } else {
                for(int i = 0; i < N; i++){
                    a[i] = cp(d[i].x, 0);
                    b[i] = cp(d[i].y, 0);
                }
            }
            d[N] = cp(0, 0);
        }
        void mult(int *a, int *b, int M){
            int N = 1, len = 0;
            while(N < M) N <<= 1, len++;
            for(int i = 0; i < N; i++)
                wz[i] = (wz[i >> 1] >> 1) | ((i & 1) << (len - 1));
            for(int i = 0; i < N; i++){
                ka[i].x = a[i] >> 15;
                kb[i].x = b[i] >> 15;
                ta[i].x = a[i] & 32767;
                tb[i].x = b[i] & 32767;
            }
            D_FFT(ta, ka, N, 1); D_FFT(tb, kb, N, 1);
            for(int i = 0; i < N; i++){
                kk[i] = ka[i] * kb[i];
                kt[i] = ka[i] * tb[i] + ta[i] * kb[i];
                tt[i] = ta[i] * tb[i];
            }
            D_FFT(tt, kk, N, -1); FFT(kt, N, -1);
            for(int i = 0; i < N; i++){
                tt[i] = tt[i] / N;
                kt[i] = kt[i] / N;
                kk[i] = kk[i] / N;
            }
        }
    }PLE;
    void main() {
        scanf("%d%d%d", &n, &m, &p); n++; m++;
        for(int i = 0; i < n; i++)  scanf("%d", &a[i]),a[i] = a[i] % p;
        for(int i = 0; i < m; i++)  scanf("%d", &b[i]),b[i] = b[i] % p;
        PLE.mult(a, b, n + m - 1);
        for(int i = 0; i < n + m - 1; i++)
            printf("%lld ",(((((LL)round(kk[i].x)) % p) << 30) + ((((LL)round(kt[i].x)) % p) << 15) + ((LL)round(tt[i].x)) % p) % p);
    }
}
```cpp struct polynie { CP getw(int m, int k, int op) { return CP(cos(2 * pi * k / m), op * sin(2 * pi * k / m)); } int wz[MAXN]; CP A[MAXN], B[MAXN], C[MAXN]; void FFT(CP *a, int N, int op) { rop(i, 0, N) if(i < wz[i]) swap(a[i], a[wz[i]]); for(int l = 2; l <= N; l <<= 1) { int mid = l >> 1; CP x, y, w, wn = CP(cos(pi / mid), sin(op * pi / mid)); for(int i = 0; i < N; i += l) { w = CP(1, 0); rop(j, 0, mid) { x = a[i + j]; y = w * a[i + j + mid]; a[i + j] = x + y; a[i + j + mid] = x - y; w = w * wn; } } } } void mult(CP *a, CP *b, CP *c, int M) { int N = 1, len = 0; while(N < M) N <<= 1, len++; rop(i, 0, N) wz[i] = (wz[i >> 1] >> 1) | ((i & 1) << (len - 1)); FFT(a, N, 1); FFT(b, N, 1); rop(i, 0, N) c[i] = a[i] * b[i]; FFT(c, N, -1); rop(i, 0, N) c[i].x = c[i].x / N, c[i].y = c[i].y / N; } void blue_stein(CP *a, int M, int op) { int M2 = M << 1; memset(A, 0, sizeof(A)); memset(B, 0, sizeof(B)); memset(C, 0, sizeof(C)); rop(i, 0, M) A[i] = a[i] * getw(M2, 1ll * i * i % M2, op); rop(i, 0, M2) B[i] = getw(M2, 1ll * (i - M) * (i - M) % M2, -op); mult(A, B, C, M2 + M - 1); rop(i, 0, M) a[i] = C[i + M] * getw(M2, 1ll * i * i % M2, op); if(op == -1) rop(i, 0, M) a[i].x = a[i].x / M, a[i].y = a[i].y / M; } }PLE; ``` ### $\color{red}{\text{多项式全集之三 多项式求逆与除法:}}

多项式求逆：

已知$A(x)$，且$F(x)A(x)\equiv 1 (mod~x^n)$，求$F(x) $$\begin{aligned}B(x)&\equiv F(x)^{-1}&(mod\,x^{\lceil \frac n 2 \rceil})\\\text{由于}G(x)&\equiv F(x)^{-1}& (mod\,x^n)\\\text{所以}G(x)&\equiv F(x)^{-1}& (mod\,x^{\lceil \frac n 2 \rceil})\\G(x)& \equiv B(x)&(mod\,x^{\lceil \frac n 2 \rceil})\\G(x)-B(x)& \equiv 0&(mod\,x^{\lceil \frac n 2 \rceil})\\\end{aligned}$$ 两边平方，得： 由于$[G(x)-B(x)]^2$的第$k<n$项为 $$\sum_{i=0}^k[g_ix^i-b_ix^i][g_{k-i}x^{k-i}-b_{k-i}x^{k-i}]$$ $i,k-i$一定有一项$<\frac n 2$，所以 $$\begin{aligned}[G(x)-B(x)]^2&\equiv 0 &(mod~x^n)\\G^2(x)+B^2(x)-2G(x)B(x)&\equiv 0&(mod~x^n)\end{aligned}$$ 两边同乘$A(x)$，得： $$\begin{aligned}A(x)G^2(x)+A(x)B^2(x)-2A(x)G(x)B(x)&\equiv 0& (mod~x^n)\\G(x)+A(x)B^2(x)-2B(x)&\equiv 0& (mod~x^n)\\G(x)&\equiv 2B(x)-A(x)B^2(x)&(mod~x^n)\\G(x)&\equiv B(x)[2-A(x)B(x)]&(mod~x^n)\\\end{aligned}$$ #### **多项式除法：** $1.$问题描述： 已知一个$n$次多项式$A(x)$，一个$m$次多项式$B(x)$，且$A(x)=B(x)C(x)+D(x)$，求$n-m$次多项式$C(x)$，$<m$次多项式$D(x)$。 $2.$推导过程： 由$A(x) = B(x)C(x)+D(x)$得： $$\begin{aligned}A(\frac 1 x)&=B(\frac 1 x)C(\frac 1 x)+D(\frac 1 x)\\x^nA(\frac 1 x)&=x^nB(\frac 1 x)C(\frac 1 x)+x^nD(\frac 1 x)\\x^nA(\frac 1 x)&=x^mB(\frac 1 x)x^{n-m}C(\frac 1 x)+x^{m-1}x^{n-m+1}D(\frac 1 x)\\A_r(x)&=B_r(x)C_r(x)+x^{n-m+1}D(x)\\A_r(x)&=B_r(x)C_r(x)&(mod \;x^{n-m+1})\\C_r(x)&=A_r(x)B_r^{-1}(x)&(mod \;x^{n-m+1})\\\end{aligned}$$ 求逆可得$B_r(x)$，再反转得$B(x)$，然后乘$C(x)$去减$A(x)$得$D(x)$. #### **喜闻乐见的模板：** ```cpp namespace INV{ typedef long long LL; int n, a[300005], b[300005]; const int mod = 998244353; const int g = 3189; int qpow(int a, int b){ int ans = 1; while(b){ if(b & 1) ans = 1ll * ans * a % mod; a = 1ll * a * a % mod; b >>= 1; } return ans; } struct p_l_e{ int wz[300005], i_c[300005]; void NTT(int *a, int N, int op){ for(int i = 0; i < N; i++) if(i < wz[i]) swap(a[i], a[wz[i]]); for(int le = 2; le <= N; le <<= 1){ int mid = le >> 1, wn = qpow(g, (mod - 1) / le); if(op == -1) wn = qpow(wn, mod - 2); for(int i = 0; i < N; i += le){ LL w = 1; int x, y; for(int j = 0; j < mid; j++){ x = a[i + j]; y = w * a[i + j + mid] % mod; a[i + j] = (x + y) % mod; a[i + j + mid] = (x - y + mod) % mod; w = w * wn % mod; } } } } int init(int M){ int N = 1, len = 0; while(N < M) N <<= 1, len++; for(int i = 0; i < N; i++) wz[i] = (wz[i >> 1] >> 1) | ((i & 1) << (len - 1)); return N; } void INV(int *a, int *b, int deg){ if(deg == 1){b[0] = qpow(a[0], mod - 2); return;} INV(a, b, (deg + 1) >> 1); int N = init(deg + deg - 1); for(int i = 0; i < deg; i++) i_c[i] = a[i]; for(int i = deg; i < N; i++) i_c[i] = 0; NTT(b, N, 1);NTT(i_c, N, 1); for(int i = 0; i < N; i++) b[i] = 1ll * b[i] * (2 - 1ll * b[i] * i_c[i] % mod + mod) % mod; NTT(b, N, -1); int t = qpow(N, mod - 2); for(int i = 0; i < N; i++) b[i] = 1ll * b[i] * t % mod; for(int i = deg; i < N; i++) b[i] = 0; } }PLE; void main(){ scanf("%d", &n); for(int i = 0; i < n; i++) scanf("%d", &a[i]); PLE.INV(a, b, n); for(int i = 0; i < n; i++) printf("%d ",b[i]); } } ``` ### $\color{red}{\text{多项式全集之四 多项式ln与exp:}}

多项式ln：

设 $$G(x) =\ln F(x)$$ 两边求导得 $$G'(x)=\frac {F'(x)} {F(x)}$$ 积分回去即可。 $2.$应用： $$e^{F(x)}=\sum_{k \ge 0}\frac{F^k(x)}{k!}=G(x)$$ $$F(x) = \ln G(x)$$ 这个的组合意义是：无序组合。 设$F(x)$，$f_i$表示一些东西，那么这些东西有序组合的方案数为 $$F^0(x) + F^1(x)+F^2(x)+\cdots=\frac 1 {1 - F(x)}$$ 而无序组成的方案数为： $$\frac{F^0(x)}{0!} + \frac {F^1(x)}{1!}+\frac{F^2(x)}{2!}+\cdots=e^{F(x)}$$ 如果无序组合方案数好求，那么求$\ln$就能得到$F(x)$。 [例题](https://www.luogu.org/problemnew/show/P4841) #### **多项式$\exp$：** 我们就是求$F(x)=e^{A(x)}$，变形一下得$\ln(F(x))-A(x)=0$，就是求函数$G(z)=\ln(z)-A(x)$的零点，带入牛顿迭代的公式得： $$\begin{aligned}F(x)&\equiv F_0(x)+\frac{G(F(x))-G(F_0(x))}{G^{(1)}(F_0(x))}&(mod\;x^n)\\\end{aligned}$$ 由于$G(F(x))\equiv 0 (mod\;x^n) \begin{aligned}F(x)\equiv F_0(x)-\frac{G(F_0(x))}{G^{(1)}(F_0(x))}&&(mod\;x^n)\\\end{aligned}

把F(x)看成自变量x,A(x)看成常数来求导

\begin{aligned}F(x)&\equiv F_0(x)-\frac{\ln F_0(x)-A(x)}{\frac{1}{F_0(x)}}&(mod\;x^n)\\F(x)&\equiv F_0(x)(1-{\ln F_0(x)+A(x)})&(mod\;x^n)\end{aligned}

对于终止情况：

\ln F(x)-A(x)=0(mod\;x)

因为都只剩下了常数项，

F(x)=e^{A(x)}

喜闻乐见的代码：

多项式\ln:

namespace PLE_ln{
  struct polyme {
      int li[SZ], wz[SZ];
      void NTT(int *a, int N, int op) {
          rop(i, 0, N) if(i < wz[i]) swap(a[i], a[wz[i]]);
          for(int l = 2; l <= N; l <<= 1) {
              int mid = l >> 1;
              int x, y, w, wn = qpow(g, (mod - 1) / l);
              if(op) wn = qpow(wn, mod - 2);
              for(int i = 0; i < N; i += l) {
                  w = 1;
                  for(int j = 0; j < mid; ++j) {
                      x = a[i + j]; y = 1ll * w * a[i + j + mid] % mod;
                      a[i + j] = (x + y) % mod;
                      a[i + j + mid] = (x - y + mod) % mod;
                      w = 1ll * w * wn % mod;
                  }
              }
          }
      }
      void qd(int *a, int *b, int n) {
          rop(i, 0, n) b[i] = 1ll * a[i + 1] * (i + 1) % mod;
      }
      void jf(int *a, int *b, int n) {
          rop(i, 1, n) b[i] = 1ll * a[i - 1] * qpow(i, mod - 2) % mod;
      }
      void mult(int *a, int *b, int *c, int M) {
          int N = 1, len = 0;
          while(N < M) N <<= 1, len ++;
          rop(i, 0, N) wz[i] = (wz[i >> 1] >> 1) | ((i & 1) << (len - 1));
          NTT(a, N, 0); NTT(b, N, 0);
          rop(i, 0, N) c[i] = 1ll * a[i] * b[i] % mod;
          NTT(c, N, 1);
          int t = qpow(N, mod - 2);
          rop(i, 0, N) c[i] = 1ll * c[i] * t % mod;
      }
      void inv(int *a, int *b, int deg) {
          if(deg == 1) {b[0] = qpow(a[0], mod - 2) % mod; return;}
          inv(a, b, (deg + 1) >> 1);
          rop(i, 0, deg) li[i] = a[i];
          int N = 1, len = 0;
          while(N < deg + deg - 1) N <<= 1, len ++;
          rop(i, 0, N) wz[i] = (wz[i >> 1] >> 1) | ((i & 1) << (len - 1));
          rop(i, deg, N) li[i] = 0;
          NTT(li, N, 0); NTT(b, N, 0);
          rop(i, 0, N) b[i] = 1ll * b[i] * (2 - 1ll * li[i] * b[i] % mod + mod) % mod;
          NTT(b, N, 1);
          int t = qpow(N, mod - 2);
          for(int i = 0; i < N; i++) b[i] = 1ll * b[i] * t % mod;
          rop(i, deg, N) b[i] = 0;
      }
  }PLE;
  int a[SZ], da[SZ], ia[SZ], dla[SZ], la[SZ], n;

  void main() {

      scanf("%d", &n);
      rop(i, 0, n) scanf("%d", &a[i]);
      PLE.qd(a, da, n);
      PLE.inv(a, ia, n);
      PLE.mult(ia, da, dla, n + n - 1);
      PLE.jf(dla, la, n);
      rop(i, 0, n) printf("%d ", la[i]);
  }
}

\color{red}{\text{多项式全集之五 多项式快速幂与开根:}}

多项式快速幂方法一：

直接套用快速幂模板，重载*运算符，O(n\log n\log W)。

多项式快速幂方法二：

#### **多项式开根：** 已知$F^2(x)\equiv A(x)\;\;\;(mod\;x^n)$，求$F(x)$。 移项得：$F^2(x)-A(x)\equiv 0\;\;\;(mod\;x^n)

问题转化成求函数G(z)=z^2-A(x)\equiv 0\;\;\;(mod\;x^n)的零点，带入牛顿迭代得：

\begin{aligned}F(x)&=\frac{G(F(x))-G(J(x))}{G^{(1)}(J(x))} + J(x)&(mod\;x^n)\\F(x)&=\frac{A(x)-J^2(x)}{2J(x)} + \frac{2J^2(x)}{2J(x)}&(mod\;x^n)\\F(x)&=\frac{A(x)+J^2(x)}{2J(x)}&(mod\;x^n)\end{aligned}

多项式求逆即可。

结束状态：F(x)=\sqrt{A_0}二次剩余一下。

\color{red}{\text{多项式全集之六 多项式多点求值与快速插值：}}

多项式多点求值 ：

给定多项式F(x)，X=\{{x_0, x_1, x_2\cdots x_m}\}，求F(x_0),F(x_1)\cdots。

设：

\begin{aligned}X_0&=\{x_0,x_1,x_2,\cdots x_{\lfloor\frac n 2\rfloor}\}\\X_1&=\{x_{\lfloor\frac n 2\rfloor+1},x_{\lfloor\frac n 2\rfloor+2},x_{\lfloor\frac n 2\rfloor+3},\cdots x_n\}\\P_0(x)&=\prod_{i=0}^{\lfloor\frac n 2\rfloor}(x-x_i)\\P_1(x)&=\prod_{i=\lfloor\frac n 2\rfloor+1}^{n}(x-x_i)\end{aligned}

显然当x\in X_0时，P_0(x)=0，P_1同理。

我们设

A_0(x)=F(x)\% P_0(x)

那么

F(x)=P_0(x)D(x)+A_0(x)

当x\in X_0时，F(x)=A_0(x)，我们求出A_0(x)的值就可以得到F(x)的值了，而且A_0(x)长度小于P_0(x)长度(\lfloor\frac n 2\rfloor)。P_1,A_1同理，这样我们把P预处理一下，就可以递归解决了。

多项式快速插值 ：

\color{red}{\text{多项式全集之八 本文提到知识的部分扩展：}}

麦克劳林级数：

泰勒公式中，取x_0=0得到的式子叫麦克劳林级数：

\sum_{i=0}^{\infty}\frac{F^{(i)}(0)}{i!}x^i