Home

# 多项式学习笔记

做多项式题就像嗑药，出多项式题就像贩毒。—— 某 FJOI 某知名选手

# 快速傅里叶变换（Fast Fourier Transform）

快速傅里叶变换，指利用计算机快速高效的计算离散傅里叶变换的方法，简称FFT。

# 离散傅里叶变换（Discrete Fourier Transform)

这是一种将多项式的系数表示法变为点值表示法的一种方法，简称 DFT，时间复杂度为优秀的 O(nlog⁡n)O(n \log n)。

我们有下面的多项式：

A(x)=a0+a1x+a2x2+a3x3+...+anxnA(x) = a_0 + a_1 x + a_2 x^2 + a_3 x^3 + ... +a_n x^n

再定义以下两个多项式：

A1(x)=a0+a2x+a4x2+...+anxn2A2(x)=a1+a3x+a5x2+...+an−1xn2−1\begin{aligned} A_1(x) &= a_0 + a_2 x + a_4 x^2 + ... + a_n x^{\frac{n}{2}} \\ A_2(x) &= a_1 + a_3 x + a_5 x^2 + ... + a_{n-1} x^{\frac{n}{2}-1} \end{aligned}

此时，我们试着用 A1A_1 和 A2A_2 表示 AA：

A(x)=A1(x2)+xA2(x2)A(x) = A_1(x^2) + x A_2(x^2)

我们引入单位根这个神奇的东西

满足 xn=1x^n = 1 的 xx 在复数域上的 nn 个解，记为 ωn\omega_n，称为单位根。
ωn1\omega_n^1 称作 nn 次单位根。
ωnk=(ωn1)k\omega_n^k = (\omega_n^1)^k
ωnk=e2πkin=cos⁡2kπn+sin⁡2kπni\omega_n^k = e^{\frac{2 \pi ki}{n}} = \cos\frac{2k\pi}{n} + \sin\frac{2k\pi}{n}i

我们不难发现单位根有许多方便的性质：

ωnk=ωtntk\omega_n^k = \omega_{tn}^{tk}
ωnk=ωnk%n\omega_n^k = \omega_{n}^{k\%n}
ωnk+n2=−ωnk\omega_n^{k+\frac{n}{2}} = -\omega_n^k

(复数域上乘法的几何意义：模长相乘，辐角相加。)

回到正题，直接做系数表示与点值表示的时间复杂度为 O(n)O(n) ，为加速变换过程，我们选择一些特殊的点值 —— ωn0\omega_n^{0}~ωnn−1\omega_n^{n-1} 这 nn 个，同时为了处理的方便， nn 的值取 22 的整次幂。

将单位根代入 A(x)A(x) 中得到：

A(ωnk)=A1(ωn2k)+ωnkA2(ωn2k)=A1(ωn2k)+ωnkA2(ωn2k)A(ωnk+n2)=A1(ωn2k)−ωnkA2(ωn2k)=A1(ωn2k)−ωnkA2(ωn2k)\begin{aligned} A(\omega_n^k) &= A_1(\omega_n^{2k}) + \omega_n^k A_2(\omega_n^{2k}) \\ &= A_1(\omega_{\frac{n}{2}}^{k}) + \omega_n^k A_2(\omega_{\frac{n}{2}}^{k}) \\ A(\omega_n^{k+\frac{n}{2}}) &= A_1(\omega_n^{2k}) - \omega_n^k A_2({\omega}_n^{2k}) \\ &= A_1(\omega_{\frac{n}{2}}^{k}) - \omega_n^k A_2(\omega_{\frac{n}{2}}^{k}) \end{aligned}

据此，我们就可以分治去算啦。

# 离散傅里叶逆变换 (Inverse Discrete Fourier Transform)

有了 DFT，我们还需要一种将点值表示转化为系数表示的算法，直接转化的时间复杂度为 O(n2)O(n^2), 这是我们不能接受的，而 IDFT 则是最佳选择。

我们有一个多项式：

A(x)=∑i=0n−1aixiA(x) = \sum_{i = 0}^{n-1}a_i x^i

将 nn 个单位根代入得到 nn 个点值，ωni\omega_n^i 对应的点值记为 yiy_i。

令 B(x)=∑i=0n−1yixiB(x) = \sum\limits _{i = 0}^{n-1} y_i x^i ，将 ωnk\omega_n^k 的倒数（共轭复数）ωn−k\omega_n^{-k} 代入：

B(ωn−k)=∑i=1n−1yiωn−ik=∑i=0n−1(∑j=0n−1ajωnij)ωn−ik=∑j=0n−1aj∑i=0n−1(ωnj−k)i=nak∴ak=B(ωn−k)n\begin{aligned} B(\omega_n^{-k}) &= \sum_{i = 1} ^{n-1} y_i \omega_n^{-ik} \\ &= \sum_{i = 0}^{n-1} \left( \sum_{j = 0}^{n-1} a_j \omega_n^{ij} \right) \omega_n^{-ik} \\ &= \sum_{j = 0}^{n-1} a_j \sum_{i = 0}^{n-1} \left( \omega_n^{j-k} \right)^i \\ &= n a_k \\ \therefore a_k &= \frac{B(\omega_n^{-k})}{n} \end{aligned}

总结一下就是：

fk=∑i=0n−1ωnigi⟺gk=1n∑i=0n−1ωn−ifif_k = \sum_{i = 0}^{n-1} \omega_n^i g_i \Longleftrightarrow g_k = \frac{1}{n} \sum_{i = 0}^{n-1} \omega_n^{-i} f_i

# 快速数论变换 (Number Theoretic Transform)

在 FFT 中，复数和正余弦函数的使用会造成不可避免的精度误差
” 在离散正交变换的理论中，已经证明在复数域内，具有循环卷积特性的唯一变换是 DFT，所以在复数域中不存在具有循环卷积性质的更简单的离散正交变换。“
"也就提出了因此提出了以数论为基础的具有循环卷积性质的快速数论变换（NTT），用有限域上的单位根来取代复平面上的单位根。"

# 快速数论变换 (Number Theoretic Transform)

原根与单位根有着一些相似的性质，说明原根或许能在一定程度上代替单位根进行多项式卷积。

一般来说，最小正原根并不大，所以可以枚举求出。

gg 为模 pp 意义下的原根的充要条件为：

gp−1pi≡1(modp)g^{\frac{p-1}{p_i}} \equiv 1 \pmod p 恒不成立，且 gp−1≡1(modp)g^{p-1} \equiv 1 \pmod p 。其中，pip_i 为 pp 的质因子

我们记 gg 为模 pp 意义下的原根，记 gn=gp−1ng_n = g^{\frac{p-1}{n}}

我们可以得到如下性质：

gnn=gn⋅p−1n=gp−1gnn2=gn2⋅p−1n=gp−12gtntk=gnkgnn≡1(modp)gnk≡gnk%n(modp)gnn2≡−1(modp)(gnk+n2)2≡gn2k+n≡gn2k(modp)\begin{aligned} g_n^n &= g^{n\cdot\frac{p-1}{n}} = g^{p-1} \\ g_n^{\frac{n}{2}} &= g^{\frac{n}{2} \cdot \frac{p-1}{n}} = g^{\frac{p-1}{2}} \\ g_{tn}^{tk} &= g_n^k \\ g_n^n &\equiv 1 &\pmod p \\ g_n^k &\equiv g_n^{k\%n} &\pmod p \\ g_n^{\frac{n}{2}} &\equiv -1 &\pmod p \\ \left( g_n^{k+\frac{n}{2}} \right)^2 &\equiv g_n^{2k+n} \equiv g_n^{2k} &\pmod{p} \\ \end{aligned}

借以上性质，我们将 gnkg_n^k 和 gnk+n2g_n^{k+\frac{n}{2}} 代入上文多项式 A(x)A(x) 、 A1(x)A_1(x) 、 A2(x)A_2(x) 与 ωnk\omega_n^k 和 ωnk+n2\omega_n^{k+\frac{n}{2}} 本质无异。

# 快速数论逆变换 (Inverse Number Theoretic Transform)

在 IDFT 中，我们将单位根的倒数（共轭复数）代入上文的多项式 B(x)B(x) 中得到转换后的系数，而在 INTT 中代入的则是原根在模 pp 意义下的逆元。

# Coding Time

#define long long long
auto NTT(long *A, int type) -> void
{
    long x, y, g0, gn;
    for(int i = 0; i < limit; i += 1)
        if(i < r[i]) swap(A[i], A[r[i]]);
    for(int midl = 1; midl < limit; midl <<= 1)
    {
        int len = midl<<1;
        g0 = ksm((type?g:invg), (p-1)/len);
        for(int j = 0; j < limit; j += len)
        {
            gn = 1;
            for(int k = 0; k < midl; k += 1, (gn *= g0) %= p)
            {
                x = A[j+k]%p;
                y = gn*A[j+k+midl]%p;
                A[j+k] = (x+y)%p;
                A[j+k+midl] = (x-y+p)%p;
            }
        }
    }
}

# 任意模数的快速数论变换 (任意模数 NTT)

在 NTT 中，模数是有讲究的，这些质数要可以写成a⋅2k+1a\cdot2^k+1 的形式。
比如：
469762049=7∗226+1,(g=3)469762049 = 7*2^{26}+1, (g = 3)
998244353=119∗223+1，(g=3)998244353 = 119*2^{23}+1， (g = 3)
1004535809=479∗221+1,(g=3)1004535809 = 479*2^{21}+1, (g = 3)
. . . . . .
而 2k2^k 则是 NTT 能够运算的最大长度。

对于题目要求的其他模数，比如 109+710^9+7，或者普通的 NTT 就显得心有余而力不足了。
（道高一尺，魔高一丈）
任意模数 NTT（也可以称作三模数 NTT）在众望中出世。

我们知道 （不知道） ： FFT 的运算结果 ≤NP2\leq N P ^2, 所以我们可以找三个大质数相乘大于这个结果，对三个模数分别计算出结果，最后中国剩余定理合并，然后模提莫要求的模数即可。

当然，直接合并这么大的数是会 Bomb 的，我们可以先合并前两个，再使用一些神奇的操作求出答案。

记最终答案为 ansans, 选择的三个质数分别为 p1,p2,p3p_1, p_2, p_3,

记三次 NTT 的结果为：

ans≡a1(modp1)ans≡a2(modp2)ans≡a3(modp3)ans \equiv a_1 \pmod {p_1} \\ ans \equiv a_2 \pmod {p_2} \\ ans \equiv a_3 \pmod {p_3}

先用 CRT 合并得到：

ans≡a4(modp1p2)ans \equiv a_4 \pmod {p_1 p_2}

转化为等式为：

ans=kp1p2+a4ans = k p_1 p_2 + a_4

若求出 kk 便可得到 ansans：

∵ans≡a3(modp3)∴kp1p2+a4≡a3(modp3)∴k≡(a3−a4)p1−1p2−1(modp3)∴ans≡kp1p2+a4(modp1p2p3)\begin{aligned} \because ans &\equiv a_3 &\pmod{ p_3} \\ \therefore k p_1 p_2 + a_4 &\equiv a_3 &\pmod {p_3} \\ \therefore k &\equiv (a_{3}-a_{4}) p_1^{-1} p_2^{-1} &\pmod{p_3} \\ \therefore ans &\equiv k p_1 p_2 + a_4 &\pmod{p_1 p_2 p_3} \end{aligned}

因为 ansans 是小于 p1p2p3p_1 p_2 p_3 的，所以ans=kp1p2+a4ans = k p_1 p_2 + a_4。（别忘记取模）

# 多项式乘法逆

定义多项式 F−1F^{-1} 为 多项式 FF 的乘法逆元，满足

F∗F−1≡1(modxn)F \ast F^{-1} \equiv 1 \pmod {x^n}

若我们已经得知 F∗H≡1(modnn2)F \ast H \equiv 1 \pmod {n^\frac{n}{2}}，需要我们求 F∗G≡1(modxn)F \ast G \equiv 1 \pmod {x^n}。

首先，F∗G≡1(modxn2)F \ast G \equiv 1 \pmod {x^{\frac{n}{2}}} 也是成立的，接下来就是推式子的时间：

∵F∗H≡1(modxn2)F∗G≡1(modxn2)∴F∗(H−G)≡0(modxn2)(H−G)≡0(modxn2)(H−G)2≡0(modxn)F∗(H2−2HG+G2)≡0(modxn)FH2−2HFG+FG2≡0(modxn)∵F∗G≡1(modxn)∴FH2−2H+G≡0(modxn)G≡2H−FH2(modxn)\begin{aligned} \because F \ast H &\equiv 1 \pmod{x^{\frac{n}{2}}} \\ F \ast G &\equiv 1 \pmod{x^{\frac{n}{2}}} \\ \therefore F \ast (H-G) &\equiv 0 \pmod{x^{\frac{n}{2}}} \\ (H-G) &\equiv 0 \pmod{x^{\frac{n}{2}}} \\ (H-G)^2 &\equiv 0 \pmod{x^n} \\ F \ast(H^2 - 2HG + G^2) &\equiv 0 \pmod{x^n} \\ FH^2 - 2 HFG + FG^2 &\equiv 0 \pmod{x^n} \\ \because F \ast G &\equiv 1 \pmod {x^n} \\ \therefore FH^2 - 2H + G &\equiv 0 \pmod{x^n} \\ G &\equiv 2H - FH^2 \pmod {x^n} \end{aligned}

# Coding Time

#define long long long
auto MUL(long *A, long *B, int len) -> void
{
    long invlim;
    limit = 1, l = 0; 
    while(limit <= (len<<1)) limit <<= 1, l += 1;
    for(int i = 1; i < limit; i += 1) 
        r[i] = (r[i>>1]>>1)|((i&1)<<(l-1));
    for(int i = len; i < limit; i += 1) A[i]=0;
    NTT(A, 1);  NTT(B, 1);
    for(int i = 0; i < limit; i += 1)
        B[i] = (2-A[i]*B[i]%p+p)%p*B[i]%p;
    NTT(B, 0); 
    invlim = ksm(limit, p-2);
    for(int i = 0; i < limit; i += 1) (B[i] *= invlim) %= p;
    for(int i = len; i < limit; i += 1) B[i] = 0;
}

# 多项式对数函数 (多项式 ln⁡\ln)

# 前置知识

# 简单导数

常数和基本初等函数导数公式：

(C)′=0(xμ)′=μxμ−1(sin⁡x)′=cos⁡x(cos⁡x)′=−sin⁡x(ax)′=ax  ln⁡a(a>0,a≠1)(ex)′=ex(log⁡ax)′=1xln⁡a  (a>0,a≠1)(ln⁡x)′=1x\begin{aligned} (C)' &= 0\\ (x^\mu)' &= \mu x ^{\mu-1}\\ (\sin x)' &= \cos x\\ (\cos x)' &= -\sin x\\ (a^x)' &= a^x \; \ln a (a > 0, a \not= 1)\\ (e^x)' &= e^x\\ (\log _a x)' &= \frac{1}{x \ln a} \; (a > 0, a \not= 1)\\ (\ln x)' &= \frac{1}{x} \end{aligned}

复合函数的求导法则：
定理 1: 如果 u=g(x)u = g(x) 在点 xx 可导，y=f(u)y = f(u) 在点 u=g(x)u = g(x)，那么复合函数 y=f(g(x))y = f(g(x)) 在点 xx 可导，且其导数为 d⁡yd⁡x=f′(u)⋅g′(x)\frac{\operatorname{d}y}{\operatorname{d}x} = f'(u) \cdot g'(x) 或 d⁡yd⁡x=d⁡yd⁡u⋅d⁡ud⁡x\frac{\operatorname{d}y}{\operatorname{d}x} = \frac{\operatorname{d}y}{\operatorname{d}u} \cdot \frac{\operatorname{d}u}{\operatorname{d}x}。

# 简单积分

定义 1: 如果在区间 I\mathbf{I} 上，可导函数 F(x)F(x) 的导数为 f(x)f(x)，即 ∀x∈I\forall x \in \mathbf{I}，都有 F′(x)=f(x)F'(x) = f(x) 或 d⁡F(x)=f(x)d⁡x\operatorname{d}F(x) = f(x)\operatorname{d}x
那么函数 F(x)F(x) 就称为 f(x)f(x) 或 f(x)d⁡xf(x)\operatorname{d}x 在区间 I\mathbf{I} 上的一个原函数。

原函数存在定理：连续函数一定有原函数

定义 2: 在区间 I\mathbf{I} 上，函数 f(x)f(x) 的带有任意常数项的原函数称为 f(x)f(x)（或 f(x)d⁡xf(x)\operatorname{d}x) 在区间 I\mathbf{I} 上的不定积分，记作 ∫f(x)d⁡x\int f(x)\operatorname{d}x。
其中，记号 ∫\int 称为积分号，f(x)f(x) 称为被积函数， f(x)d⁡xf(x)\operatorname{d}x 称为被积表达式，xx 称为积分变量。

# 多项式对数函数

我们有多项式 F(x)F(x) 和 G(x)G(x)，记 G≡ln⁡F(modxn)G \equiv {\ln F} {\pmod{x^n}}

推式子时间：

G≡ln⁡F(modxn)G′≡(ln⁡F)′(modxn)G′≡(ln⁡′F)∗F′(modxn)G′≡F′F(modxn)\begin{aligned} G &\equiv \ln F &\pmod{x^n} \\ G' &\equiv (\ln F)' &\pmod{x^n} \\ G' &\equiv (\ln' F) \ast F' &\pmod{x^n} \\ G' &\equiv \frac{F'}{F} &\pmod{x^n} \\ \end{aligned}

我们就可以先把 G′G' 搞出来，然后再积回去，就是答案啦。

# Coding Time

#define long long long
auto main() -> signed
{
    scanf("%d", &n);
    for(int i = 0; i < n; i += 1)
        scanf("%d", &a[i]);
    inv(b, n); 
    der(a, n);
    NTT::MUL(a, b, n);
    itg(a, n);
    for(int i = 0; i < n; i += 1)
        printf("%lld ", a[i]);
    return 0;
}
auto inv(long *B, int n) -> void
{
    if(n == 1) return (void)(B[0] = ksm(a[0], p-2));
    inv(B, (n+1)>>1);
    for(int i = 0; i < n; i += 1) c[i] = a[i];
    NTT::INV(c, B, n);
} 
// 导数 derivative 
auto der(long *A, int n) -> void
{
    for(int i = 1; i < n; i += 1)
        A[i-1] = A[i]*(i)%p;
    A[n-1] = A[n] = 0;
}
// 积分 integral
auto itg(long *A, int n) -> void
{
    for(int i = n; i >= 1; i -= 1)
        A[i] = A[i-1]*ksm(i, p-2)%p;
    A[0] = 0;
}

# 多项式指数函数 (多项式 exp⁡\exp)

# 前置知识

# 泰勒展开 (Taylor)

泰勒中值定理 1: 如果函数在 x0x_0 处有 nn 阶导数，那么存在 x0x_0 的一个邻域，对于该邻域内任一 xx，有
f(x)=f(x0)+f′(x0)(x−x0)+f′′(x0)2!(x−x0)2+...+f(n)(x0)n!(x−x0)n+Rn(x)f(x) = f(x_0) +f'(x_0)(x-x_0) + \frac{f''(x_0)}{2!} (x-x_0)^2+...+\frac{f^{(n)}(x_0)}{n!}(x-x_0)^n + R_n(x)
其中，Rn(x)=ο((x−x0)n)R_n(x) = \omicron((x-x_0)^n).

泰勒展开在 x0x_0 处展开的封闭形式如下: ∑n=0+∞f(n)(x0)n!(x−x0)n\sum\limits_{n=0}^{+\infty} \frac{f^{(n)}(x_0)}{n!}(x-x_0)^n

泰勒展开

# 牛顿迭代 (Newton's method)

牛顿迭代用于求函数零点，通过不断地切线逼近所求值，但最终也只是近似值，迭代的次数越多，精确度越高，误差越小。

若我们要求函数 f(x)f(x) 的零点，初始近似值为 x0x_0，切线方程为 y=f′(x0)(x−x0)+f(x0)y = f'(x_0) (x - x_0) + f(x_0)
令 y=0y = 0，得 x=x0−f(x0)f′(x0)x = x_0 - \frac{f(x_0)}{f'(x_0)}.

还需提的一点是，在多项式中，每迭代一次，精度翻倍。

牛顿迭代

# 多项式指数函数

我们有多项式 F(x)F(x) 和 A(x)A(x)， 记多项式的对数函数为 F≡eA(modxn)F \equiv e^A \pmod {x^n}.

对两边分别进行 ln⁡\ln 运算，得到 ln⁡F≡A(modxn)\ln F \equiv A \pmod {x^n}.

移项，得 ln⁡F−A≡0(modxn)\ln {F} - A {\equiv 0 \pmod {x^n}}

记 G(F)=ln⁡F−AG(F) = \ln F - A，有 G(F)≡0(modxn)G(F) \equiv 0 \pmod{x^n}
另外。我们有 G(F0)≡0(modxn2)G(F_0) \equiv 0 \pmod{x^{\frac{n}{2}}}。

这时，我们就可以带入牛顿迭代中求解，推式子时间到:

G(F0)≡0(modxn2)F≡F0−G(F0)G′(F0)(modxn)F≡F0−F0∗G(F0)(modxn)F≡F0−F0(ln⁡F0−A)(modxn)F≡F0(A−ln⁡F0+1)(modxn)\begin{aligned} G(F_0) &\equiv 0 \pmod{x^{\frac{n}{2}}} \\ F &\equiv F_0 - \frac{G(F_0)}{G'(F_0)} &\pmod{x^n} \\ F &\equiv F_0 - F_0 \ast G(F_0) &\pmod{x^n} \\ F &\equiv F_0 - F_0(\ln F_0 - A) &\pmod{x^n} \\ F &\equiv F_0(A - \ln F_0 + 1) &\pmod{x^n} \\ \end{aligned}

又可以递归求解啦

# Coding Time

#define long long long
auto exp(long *B, int n) -> void
{
    if(n == 1) return (void)(B[0] = 1);
    exp(B, (n+1)>>1);
    memset(c, 0, sizeof(c));
    lnz(B, c, n);
    for(int i = 0; i < n; i += 1) 
        c[i] = (a[i]-c[i]+p)%p;
    c[0] += 1;
    NTT::MUL(B, c, n);
}
auto lnz(long *A, long *B, int n) -> void
{
    memset(d, 0, sizeof(d)); inv(A, B, n); 
    memset(d, 0, sizeof(d)); der(A, d, n);
    NTT::MUL(B, d, n);
    itg(B, n);
}
auto inv(long *A, long *B, int n) -> void
{
    if(n == 1) return (void)(B[0] = ksm(A[0], p-2));
    inv(A, B, (n+1)>>1);
    for(int i = 0; i < n; i += 1) d[i] = A[i];
    NTT::INV(d, B, n);
}
auto der(long *A, long *B, int n) -> void
{
    for(int i = 1; i < n; i += 1)
        B[i-1] = A[i]*i%p;
    B[n-1] = B[n] = 0;
}
auto itg(long *A, int n) -> void
{
    for(int i = n; i >= 1; i -= 1)
        A[i] = A[i-1]*ksm(i, p-2)%p;
    A[0] = 0;
}

# 多项式开根

我们有多项式 F(x)F(x) 和 A(x)A(x)，满足 F2(x)≡A(x)(modxn)F^2(x) \equiv A(x) \pmod{x^n}.

要去求多项式 F(x)F(x)，首先我们要有一个多项式 H(x)H(x) 满足 H2(x)≡A(x)(modxn2)H^2(x) \equiv A(x) \pmod {x^{\frac{n}{2}}}.

显然，F2(x)≡A(x)(modxn2)F^2(x) \equiv A(x) \pmod {x^{\frac{n}{2}}} 成立。

记 G(F)=F2−AG(F) = F^2 - A.

又是令人愉悦的推式子时间：

F2(x)≡A(x)(modxn)H2(x)≡A(x)(modxn2)H2−A≡0(modxn2)G(H)≡0(modxn2)F≡H−G(H)G′(H)(modxn)F≡H−H2−A2H(modxn)F≡H2+A2H(modxn)\begin{aligned} F^2(x) &\equiv A(x) &\pmod {x^n} \\ H^2(x) &\equiv A(x) &\pmod {x^{\frac{n}{2}}} \\ H^2 - A &\equiv 0 &\pmod {x^{\frac{n}{2}}} \\ G(H) &\equiv 0 &\pmod{x^{\frac{n}{2}}} \\ F &\equiv H - \frac{G(H)}{G'(H)} &\pmod{x^n} \\ F &\equiv H - \frac{H^2 - A}{2H} &\pmod{x^n} \\ F &\equiv \frac{H^2 + A}{2H} &\pmod{x^n} \end{aligned}

还是递归下去，求个逆元乘起来就好啦。

# Coding Time

#define long long long
auto sqrt(long *A, long *B, int n) -> void
{
    if(n == 1) return (void)(B[0] = 1);
    sqrt(A, B, (n+1)>>1);
    memset(c, 0, sizeof(c));
    inv(B, c, n);
    for(int i = 0; i < n; i += 1) d[i] = A[i];
    NTT::SRT(d, B, c, n);
}
auto inv(long *A, long *B, int n) -> void
{
    if(n == 1) return (void)(B[0] = ksm(a[0], p-2));
    inv(A, B, (n+1)>>1);
    for(int i = 0; i < n; i += 1) d[i] = A[i];
    NTT::INV(d, B, n);
}

还有 NTT::SRT 里的 :

#define long long long
auto SRT(long *A, long *B, long *C, int len) -> void
{
    limit = 1, l = 0;
    while(limit <= (len<<1)) limit <<= 1, l += 1;
    for(int i = len; i < limit; i += 1) d[i] = 0;
    for(int i = 1; i < limit; i += 1) 
        r[i] = (r[i>>1]>>1)|((i&1)<<(l-1));
    NTT(A, 1); NTT(B, 1); NTT(C, 1);
    for(int i = 0; i < limit; i += 1)
        B[i] = (B[i]*B[i]%p+A[i]%p)*C[i]%p*inv2%p;
    NTT(B, 0); liminv = ksm(limit, p-2);
    for(int i = 0; i < limit; i += 1)
        (B[i] *= liminv) %= p;
    for(int i = len; i < limit; i += 1) B[i] = 0;
}

To  be  continued...\Huge{To \; be \; continued...}

牛顿迭代图片摘自此篇博客