【luogu P5205】【模板】多项式开根（牛顿迭代）（多项式）
【模板】多项式开根
题目链接：luogu P5205
题目大意
给你一个多项式 A(x)，要你找到一个多项式 B(x) 使得 B^2(x)=A(x) (mod x^n)。
如果有多个解选零次项系数最小的，然后保证 A(0)=1。
思路
在知道如果弄之前，我们要知道一个叫做牛顿迭代的东西。
那我们就要从泰勒展开开始，于是我们就讲讲求导，积分这些东西。
求导与积分
求导和积分是两个互逆的操作，先不收多项式的吧。
求导就是求图像在某个位置的斜率，积分就是求图像的面积。
然后在多项式就是这样：
求导：\(F'(x)=\sum\limits_{i=0}F[i]ix^{i-1}\)
积分：\(F'(x)=C+\sum\limits_{i=0}F[i]x^{i+1}/i\)
然后这里简单贴一下代码：
void dao(int *f, int m) {
for (int i = 1; i < m; i++)
f[i - 1] = mul(f[i], i);
f[m - 1] = 0;
}
void ji(int *f, int m) {
for (int i = m; i >= 1; i--)
f[i] = mul(f[i - 1], inv[i]);
f[0] = 0;
}
然后多项式的复合就直接是把 \(x\) 带入到 \(G(x)\)，就可以得到：
\(F(G(x))=\sum\limits_{i=0}F[i]G^i(x)\)
多项式牛顿迭代
那在知道牛顿迭代之前我们要知道一个叫做泰勒展开的东西。
我们设 \(F^{(n)}(x)\) 为 \(f(x)\) 的 \(n\) 阶导数，那么就有泰勒展开的式子：
\(F(x)=\sum\limits_{i=0}\dfrac{F^{(i)}(a)}{i!}(x-a)^i\)
（不要问我怎么证，我也不会）
（大概就是你要推这个函数，你就尝试着不断拟合，即“得寸进尺”，让它所有阶导的结果都一样）
（所以你这个 \(i\) 弄到越后拟合的就越像）
然后就有一些东西：
麦克劳林级数：其实就是 \(a=0\) 的泰勒展开：
\(F(x)=\sum\limits_{i=0}\dfrac{F^{(i)}(a)}{i!}x^i\)
这个东西最多的应用就是用卷积的形式来定义一些东西，就比如 \(\exp x\) 即 \(e^x=\sum\limits_{i=0}\dfrac{x^i}{i!}\)。（\(\exp x\) 的导数是它自己）
然后就是多项式牛顿迭代了：
首先给一下式子：如果已知 \(G(F(x))=0\)，而且已经得知 \(G(F_*(x))=0\pmod{x^{\frac{n}{2}}}\)，那么就有 \(F(x)=F_*(x)-\dfrac{G(F_*(x))}{G'(F_*(x))}\)。
一般来讲题目中的 \(G\) 函数都会比较简单，可以手动分析爆算。
然后不难看出可以倍增出来，常数项就特殊考虑一下。
然后这里给给证明：
首先 \(F(x)\ \bmod x^n\) 显然也可以作为 \(\bmod x^{\frac{n}{2}}\) 的一组解，所以 \(F(x)=F_*(x)\pmod{x^{\frac{n}{2}}}\)
然后我们把 \(G(F(x))\) 在 \(F_*(x)\) 的位置泰勒展开：
\(G(F(x))=G(F_*(x))+\dfrac{G'(F_*(x))}{1!}(F(x)-F_*(x))+\dfrac{G''(F_*(x))}{2!}(F(x)-F_*(x))^2+...\)
然后根据上面那个相等的式子可以发现 \(F_*(x)\) 和 \(F(x)\) 的后 \(\dfrac{n}{2}\) 项是相同的，所以 \((F(x)-F_*(x))\) 的最低项次也是 \(x^\frac{n}{2}\)。
然后 \((F(x)-F_*(x))^2\) 的最低项次是 \(x^n\)，然后依次就增加下去，然后你会发现从二次方开始因为你是 \(\pmod{x^n}\)，这些项都被模掉变成 \(0\) 了。
然后就只剩下前面的两项，就是 \(G(F(x))=G(F_*(x))+\dfrac{G'(F_*(x))}{1!}(F(x)-F_*(x))\)
然后你整理一下，就得到了：
\(F(x)=F_*(x)-\dfrac{G(F_*(x))}{G'(F_*(x))}\)
然后这里有一些性质以及特点：
因为 \(G(F_*(x))\) 的最低次项至少是 \(x^{\frac{n}{2}}\)，所以分母 \(G'(F_*(x))\) 的精度只需要 \(x^{\frac{n}{2}}\)。
而且记得搞的时候自由元是 \(F(x)\)。
求逆元
我们从 \(B(x)*A(x)=1\pmod{x^n}\)，然后 \(A(x)\) 已知。
那我们可以设 \(G(x)\)，使得 \(G(B(x))=A(x)B(x)-1\pmod{x^n}\)，那这里 \(A(x)\) 就是系数。
然后我们可以导一下：\(G'(B(x))=A(x)\pmod{x^n}\)
那我们同样设 \(B_*(x)\) 为 \(\pmod{x^{\frac{n}{2}}}\) 的解。
然后就牛顿迭代：
\(B(x)=B_*(x)-\dfrac{G(B_*(x))}{G'(B_*(x))}=B_*(x)-\dfrac{A(x)B_*(x)-1}{A(x)}\)
然后分母 \(\dfrac{1}{A(x)}\) 的精度只用达到 \(x^{\frac{n}{2}}\)，所以我们可以直接用 \(B_*(x)\) 代替它，因为是逆元嘛。
\(B(x)=B_*(x)-(A(x)B_*(x)-1)B_*(x)=2B_*(x)-B_*(x)^2A(x)\)
这个跟我们用平方的方法推出来的是一样的。
求开方
类似的方法，设 \(B(x)^2-A(x)=0\)，\(A(x)\) 已知。
设 \(G(x)\)，\(G(B(x))=B(x)^2-A(x)\)，然后 \(G'(B(x))=2B(x)\)。
上牛叠：\(B(x)=B_*(x)-\dfrac{G(B_*(x))}{G'(B_*(x))}=B_*(x)-\dfrac{B_*(x)^2-A(x)}{2B_*(x)}=\dfrac{A(x)+B_*(x)^2}{2B_*(x)}\)
然后就可以用这个来倍增。
然后一开始的常数项就是要看给的 \(A(x)\) 的常数项，就是它的二次剩余。
然而这道题直接说了常数项式 \(1\)，那就可以直接搞了。
代码
#include<cstdio>
#include<cstring>
#include<algorithm>
#define mo 998244353
#define clr(f, n) memset(f, 0, sizeof(int) * (n))
#define cpy(f, g, n) memcpy(f, g, sizeof(int) * (n))
using namespace std;
int n, a[800001], an[800001];
int G, Gv;
int read() {
int re = 0, zf = 1;
char c = getchar();
while (c < '0' || c > '9') {
if (c == '-') zf = -zf; c = getchar();
}
while (c >= '0' && c <= '9') {
re = (re << 3) + (re << 1) + c - '0';
c = getchar();
}
return re * zf;
}
void writee(int x) {
if (x > 9) writee(x / 10);
putchar(x % 10 + '0');
}
int write(int x) {
if (x < 0) putchar('-'), x = -x;
writee(x);
}
int mul(int x, int y) {
return 1ll * x * y % mo;
}
int jia(int x, int y) {
return (x + y) % mo;
}
int jian(int x, int y) {
return (x - y + mo) % mo;
}
int ksm(int x, int y) {
int re = 1;
while (y) {
if (y & 1) re = mul(re, x);
x = mul(x, x); y >>= 1;
}
return re;
}
void NTT(int *f, int op, int limit) {
for (int i = 0; i < limit; i++)
if (an[i] > i) swap(f[i], f[an[i]]);
for (int mid = 1; mid < limit; mid <<= 1) {
int Wn = ksm(op == 1 ? G : Gv, (mo - 1) / (mid << 1));
for (int R = (mid << 1), j = 0; j < limit; j += R) {
int w = 1;
for (int k = 0; k < mid; k++, w = mul(w, Wn)) {
int x = f[j | k], y = mul(w, f[j | mid | k]);
f[j | k] = jia(x, y); f[j | mid | k] = jian(x, y);
}
}
}
if (op == -1) {
int invl = ksm(limit, mo - 2);
for (int i = 0; i < limit; i++)
f[i] = mul(f[i], invl);
}
}
void px(int *f, int *g, int n) {
for (int i = 0; i < n; i++)
f[i] = mul(f[i], g[i]);
}
void get_an(int limit, int l_size) {
for (int i = 0; i < limit; i++) an[i] = (an[i >> 1] >> 1) | ((i & 1) << (l_size - 1));
}
void cheng(int *f, int *g, int n, int m) {
static int tmp[800001];
cpy(tmp, g, m);
int limit = 1, l_size = 0;
for (; limit < n + m - 1; limit <<= 1, l_size++);
get_an(limit, l_size);
NTT(f, 1, limit); NTT(tmp, 1, limit);
px(f, tmp, limit); NTT(f, -1, limit);
clr(f + n, limit - n); clr(tmp, limit);
}
void invp(int *f, int m) {
static int tmp[800001], r[800001], w[800001];
int n = 1, l_size = 0;
for (; n < m; n <<= 1);
w[0] = ksm(f[0], mo - 2);
for (int len = 2; len <= n; len <<= 1) {
l_size++;
get_an(len, l_size);
for (int i = 0; i < (len >> 1); i++) r[i] = w[i];
cpy(tmp, f, len);
NTT(tmp, 1, len); NTT(r, 1, len); px(r, tmp, len); NTT(r, -1, len);
clr(r, len >> 1);
cpy(tmp, w, len);
NTT(tmp, 1, len); NTT(r, 1, len); px(r, tmp, len); NTT(r, -1, len);
for (int i = len >> 1; i < len; i++)
w[i] = jian(jia(w[i], w[i]), r[i]);
}
cpy(f, w, m); clr(w, n); clr(r, n); clr(tmp, n);
}
void pow2p(int *f, int limit, int l_size) {
get_an(limit, l_size);
NTT(f, 1, limit); px(f, f, limit); NTT(f, -1, limit);
}
void sqrtp(int *f, int m) {
static int b1[800001], b2[800001];
int limit = 1;
for (; limit < m; limit <<= 1);
b1[0] = 1; int l_size = 0;
for (int len = 2; len <= limit; len <<= 1) {
l_size++;
for (int i = 0; i < (len >> 1); i++) b2[i] = jia(b1[i], b1[i]);
invp(b2, len);
pow2p(b1, len, l_size);
for (int i = 0; i < len; i++) b1[i] = jia(b1[i], f[i]);
cheng(b1, b2, len, len);
}
cpy(f, b1, m); clr(b1, limit + limit); clr(b2, limit + limit);
}
int main() {
// freopen("read.txt", "r", stdin);
G = 3; Gv = ksm(G, mo - 2);
n = read();
for (int i = 0; i < n; i++) a[i] = read();
sqrtp(a, n);
for (int i = 0; i < n; i++) write(a[i]), putchar(' ');
return 0;
}
