目录
@0 - 参考资料@
@0.5 - 多项式平方根@
已知一个多项式 \(A(x)\)，求一个多项式 \(F(x)\)，使得 \(F^2(x) = A(x)\)。
我们仿照多项式求逆的方法，加上一个模数变成 \(F^2(x) = A(x) \mod x^n\)。
然后考虑倍增，假如我们已知：
\[F_0^2(x) \equiv A(x) \mod x^{\lceil\frac{n}{2}\rceil}\]
套路移项，得到：
\[F_0^2(x)-A(x) \equiv 0 \mod x^{\lceil\frac{n}{2}\rceil}\]
套路平方，得到：
\[F_0^4(x)-2F_0^2(x)A(x)+A^2(x) \equiv 0 \mod x^{n}\]
你看右边空荡荡的，心里肯定有点儿不爽对吧：
\[F_0^4(x)+2F_0^2(x)A(x)+A^2(x) \equiv 4F_0^2(x)A(x) \mod x^{n}\]
你看右边已经有个 \(A(x)\) 了，我们就再把 \(A(x)\) 单独孤立出来吧。
\[\dfrac{F_0^4(x)+2F_0^2(x)A(x)+A^2(x)}{4F_0^2(x)} \equiv A(x) \mod x^{n}\]
你看左边已经是一个平方式了哦~
我们令 \(F(x) = \dfrac{F_0^2(x)+A(x)}{2F_0(x)}\)，则 \(F^2(x)\equiv A(x) \mod x^{n}\)。
然后就可以一路倍增上去了。边界情况下是一个模意义下的平方剩余。
这个多项式开方和多项式求逆之间有许多共通的思想，但有些地方是很具有构造性的。
但如果题目要我们求立方根呢？如果题目要我们解高次方程呢？我们这样推导肯定是不合适的。
我们是否可以找到一个通用的方法，对于这两个，甚至更多的问题都适用呢？
@1 - 牛顿迭代法@
@数学上的定义@
本节跟我们要讲的没有任何关系。
对于某一个函数 \(f(x)\) 与函数上某一个点 \((x_0, f(x_0))\)，我们可以用过该点的函数切线解析式去近似地代替整个函数。
根据导函数的定义，过 \((x_0, f(x_0))\) 的函数切线解析式为 \(g(x) = f(x_0)+f'(x_0)*(x-x_0)\)。则我们可以用 \(g(x)\) 去近似替代 \(f(x)\)。
假如 \(g(x)\) 的零点为 \(g(z)=0\)（\(z = x_0 - \frac{f(x_0)}{f'(x_0)}\)），因为 \(g(x)\) 是 \(f(x)\) 的近似，所以也许 \(f(z) = 0\)。
但是这个概率非常的低。于是我们再选择 \((z,f(z))\) 进行迭代，求出 \((z,f(z))\) 对应的 \(g(x)\)，重复十几二十百遍，应该就非常非常接近零点了。
然而……这个方法并不具有普适性，有些时候迭代会出现循环。但是对于多项式的相关操作，它的确是普适的。
根据泰勒展开：
\[f(x)=f(x_0)+f'(x_0)(x-x_0)+\dfrac{f''(x_0)(x-x_0)^2}{2!}+\dots\]
我们只取前两项（因为它们是线性的），也可以得到上面的结果。
@对于多项式的定义@
总结多项式逆元与多项式开方，我们可以得到通用的一些东西。
对于一个自变量和因变量都是多项式的函数 \(G(F(x))\)，我们要找到多项式 \(F_0(x)\)，使得：
\[G(F_0(x))=0\]
当 \(G(F(x)) = A(x)*F(x) - 1\) 时相当于逆元，当 \(G(F(x)) = F^2(x)-A(x)\) 时相当于开方。
怎么弄呢？我们套路取模，转换为 \(G(F(x))\equiv 0\mod x^n\)。
再套路倍增，假如我们已知：
\[G(F_0(x))\equiv 0\mod x^{\lceil\frac{n}{2}\rceil}\]
将 \(G(F(x))\) 在 \(F_0(x)\) 处泰勒展开，得到：
\[G(F(x))\equiv G(F_0(x))+G'(F_0(x))(F(x)-F_0(x))+\dfrac{G''(F_0(x))(F(x)-F_0(x))^2}{2!}+\dots\mod x^n\]
因为在模意义下 \(G(F(x)) \equiv 0\mod x^{\lceil\frac{n}{2}\rceil}\) 只有唯一解，所以：
\[F(x)\equiv F_0(x)\mod x^{\lceil\frac{n}{2}\rceil}\]
然后套路移项 + 套路平方：
\[(F(x)-F_0(x))^2\equiv 0\mod x^n\]
因为泰勒展开第三项及其以后的项都含 \((F(x)-F_0(x))^2\)，所以我们可以直接忽略。
又因为根据定义 \(G(F(x))\equiv 0\mod x^n\)，故：
\[0\equiv G(F_0(x))+G'(F_0(x))(F(x)-F_0(x))\mod x^n\]
变一下形：
\[F(x)\equiv F_0(x)-\dfrac{G(F_0(x))}{G'(F_0(x))}\mod x^n\]
好的我们就可以倍增了。
@2 - 牛顿迭代的应用@
@重新推导 - 多项式逆元@
对于 \(G(F(x)) = A(x)*F(x) - 1\) ，求 \(G(F(x))\equiv 0\mod x^n\) 的 \(F(x)\)。
根据迭代式：
\[F(x)\equiv F_0(x)-\dfrac{G(F_0(x))}{G'(F_0(x))}\equiv F_0(x)-\dfrac{A(x)*F_0(x) - 1}{A(x)}\mod x^n\]
根据严密的推导，我们可以得到：\(\dfrac{A(x)*F_0(x) - 1}{A(x)}\equiv F_0(x)(A(x)*F_0(x) - 1)\mod x^n\)
所以：
\[F(x)\equiv 2F_0(x)-A(x)*F_0^2(x)\mod x^n\]
时间复杂度 \(O(n\log n)\)。
@重新推导 - 多项式平方根@
对于 \(G(F(x)) = F^2(x)-A(x)\) ，求 \(G(F(x))\equiv 0\mod x^n\) 的 \(F(x)\)。
根据迭代式：
\[F(x)\equiv F_0(x)-\dfrac{G(F_0(x))}{G'(F_0(x))}\equiv F_0(x)-\dfrac{F_0^2(x) - A(x)}{2F_0(x)}\mod x^n\]
稍微恒等变换一下，就可以得到我们上边那个公式。
时间复杂度 \(O(n\log n)\)。
@多项式对数函数@
这个其实和牛顿迭代没啥关系。。。
已知 \(F(x) \equiv \ln A(x)\mod x^n\) ，求 \(F(x)\)。
怎么去理解呢？其实它就是个泰勒展开：
\[\ln(1-P(x))=-P(x)-\frac{P(x)^2}{2}-\frac{P(x)^3}{3}-\dots\]
直接求不好做的，我们考虑两边求个导数：
\[F'(x) \equiv \frac{A'(x)}{A(x)}\mod x^n\]
这样我们就可以顺利求出 \(F'(x)\)，再对 \(F'(x)\) 积分一下就可以了。
但是不定积分的常数项？一般来说题目会保证 \(A(x)\) 的常数项为 1，根据泰勒展开可得 \(F(x)\) 的常数项为 0。
时间复杂度 \(O(n\log n)\)。
@多项式指数函数@
已知 \(F(x) \equiv e^{A(x)}\mod x^n\) ，求 \(F(x)\)。
如果我们求导是没啥用的，毕竟 \(e^x\) 的导数还是它自己。
我们考虑两边求对数：
\[\ln F(x) \equiv A(x)\mod x^n\]
然后记 \(G(F(x)) = \ln F(x)-A(x)\)，于是又可以愉快地进行牛顿迭代啦。
\[F(x)\equiv F_0(x)-\dfrac{G(F_0(x))}{G'(F_0(x))}\equiv F_0(x)-(\ln F_0- A(x))*F_0(x)\mod x^n\]
时间复杂度 \(O(n\log n)\)。
@多项式幂函数@
\[F^k(x)=(e^{\ln F(x)})^k=e^{k\ln F(x)}\]
时间复杂度 \(O(n\log n)\)。
@多项式三角函数@
\[e^{i*F(x)}=\cos F(x) + i*\sin F(x)\]
把运算改成支持复数的运算。
时间复杂度 \(O(n\log n)\)。
@3 - 一些参考代码@
@对数函数@
本代码为 luoguP4725 的 AC 代码。
#include<cstdio>
#include<algorithm>
using namespace std;
const int G = 3;
const int MAXN = 400000;
const int MOD = 998244353;
int pow_mod(int b, int p) {
int ret = 1;
while( p ) {
if( p & 1 ) ret = 1LL*ret*b%MOD;
b = 1LL*b*b%MOD;
p >>= 1;
}
return ret;
}
int inv[MAXN + 5];
void init() {
inv[1] = 1;
for(int i=2;i<=MAXN;i++)
inv[i] = 1LL*(MOD - MOD/i)*inv[MOD%i]%MOD;
}
struct Polynomial{
void poly_copy(int *A, int *B, int n) {
for(int i=0;i<n;i++)
A[i] = B[i];
}
void poly_clear(int *A, int l, int r) {
for(int i=l;i<r;i++)
A[i] = 0;
}
void ntt(int *A, int n, int type) {
for(int i=0,j=0;i<n;i++) {
if( i < j ) swap(A[i], A[j]);
for(int l=(n>>1);(j^=l)<l;l>>=1);
}
for(int s=2;s<=n;s<<=1) {
int t = (s>>1);
int u = (type == 1) ? pow_mod(G, (MOD-1)/s) : pow_mod(G, (MOD-1)-(MOD-1)/s);
for(int i=0;i<n;i+=s) {
for(int j=0,p=1;j<t;j++,p=1LL*p*u%MOD) {
int x = A[i+j], y = 1LL*A[i+j+t]*p%MOD;
A[i+j] = (x + y)%MOD, A[i+j+t] = (x + MOD - y)%MOD;
}
}
}
if( type == -1 ) {
for(int i=0;i<n;i++)
A[i] = 1LL*A[i]*inv[n]%MOD;
}
}
int tmp1[MAXN + 5];
void poly_inv(int *A, int *B, int n) {
if( n == 1 ) {
B[0] = pow_mod(A[0], MOD-2);
return ;
}
int len; for(len = 1;len < (n<<1);len <<= 1);
poly_inv(A, B, (n + 1) >> 1);
poly_copy(tmp1, A, n); poly_clear(tmp1, n, len);
ntt(tmp1, len, 1); ntt(B, len, 1);
for(int i=0;i<len;i++) B[i] = 1LL*B[i]*(2 + MOD - 1LL*tmp1[i]*B[i]%MOD)%MOD;
ntt(B, len, -1); poly_clear(B, n, len);
poly_clear(tmp1, 0, len);
}
int tmp2[MAXN + 5], tmp3[MAXN + 5];
void poly_mul(int *A, int *B, int *C, int n, int m) {
int len; for(len = 1;len < (n+m-1);len <<= 1);
poly_copy(tmp2, A, n), poly_copy(tmp3, B, m);
ntt(tmp2, len, 1), ntt(tmp3, len, 1);
for(int i=0;i<len;i++) C[i] = 1LL*tmp2[i]*tmp3[i]%MOD;
ntt(C, len, -1);
poly_clear(C, n+m-1, len);
poly_clear(tmp2, 0, len); poly_clear(tmp3, 0, len);
}
void poly_int(int *A, int *B, int n) {
for(int i=n-1;i>=0;i--)
A[i+1] = 1LL*B[i]*inv[i+1]%MOD;
A[0] = 0;
}
void poly_deri(int *A, int *B, int n) {
for(int i=1;i<n;i++)
A[i-1] = 1LL*B[i]*i%MOD;
}
int tmp4[MAXN + 5];
void poly_ln(int *A, int *B, int n) {
poly_inv(A, B, n); poly_deri(tmp4, A, n);
poly_mul(tmp4, B, B, n, n);
poly_int(B, B, n-1);
}
}oper;
int f[MAXN + 5], g[MAXN + 5];
int main() {
init(); int n;
scanf("%d", &n);
for(int i=0;i<n;i++)
scanf("%d", &f[i]);
oper.poly_ln(f, g, n);
for(int i=0;i<n;i++)
printf("%d ", g[i]);
puts("");
}@指数函数@
本代码为 luoguP4726 的 AC 代码。
#include<cstdio>
#include<algorithm>
using namespace std;
const int G = 3;
const int MAXN = 400000;
const int MOD = 998244353;
int pow_mod(int b, int p) {
int ret = 1;
while( p ) {
if( p & 1 ) ret = 1LL*ret*b%MOD;
b = 1LL*b*b%MOD;
p >>= 1;
}
return ret;
}
int inv[MAXN + 5];
void init() {
inv[1] = 1;
for(int i=2;i<=MAXN;i++)
inv[i] = 1LL*(MOD - MOD/i)*inv[MOD%i]%MOD;
}
struct Polynomial{
void poly_copy(int *A, int *B, int n) {
for(int i=0;i<n;i++)
A[i] = B[i];
}
void poly_clear(int *A, int l, int r) {
for(int i=l;i<r;i++)
A[i] = 0;
}
void ntt(int *A, int n, int type) {
for(int i=0,j=0;i<n;i++) {
if( i < j ) swap(A[i], A[j]);
for(int l=(n>>1);(j^=l)<l;l>>=1);
}
for(int s=2;s<=n;s<<=1) {
int t = (s>>1);
int u = (type == 1) ? pow_mod(G, (MOD-1)/s) : pow_mod(G, (MOD-1)-(MOD-1)/s);
for(int i=0;i<n;i+=s) {
for(int j=0,p=1;j<t;j++,p=1LL*p*u%MOD) {
int x = A[i+j], y = 1LL*A[i+j+t]*p%MOD;
A[i+j] = (x + y)%MOD, A[i+j+t] = (x + MOD - y)%MOD;
}
}
}
if( type == -1 ) {
for(int i=0;i<n;i++)
A[i] = 1LL*A[i]*inv[n]%MOD;
}
}
int tmp1[MAXN + 5];
void poly_inv(int *A, int *B, int n) {
if( n == 1 ) {
B[0] = pow_mod(A[0], MOD-2);
return ;
}
int len; for(len = 1;len < (n<<1);len <<= 1);
poly_inv(A, B, (n + 1) >> 1);
poly_copy(tmp1, A, n);
ntt(tmp1, len, 1); ntt(B, len, 1);
for(int i=0;i<len;i++) B[i] = 1LL*B[i]*(2 + MOD - 1LL*tmp1[i]*B[i]%MOD)%MOD;
ntt(B, len, -1); poly_clear(B, n, len);
poly_clear(tmp1, 0, len);
}
int tmp2[MAXN + 5], tmp3[MAXN + 5];
void poly_mul(int *A, int *B, int *C, int n, int m) {
int len; for(len = 1;len < (n+m-1);len <<= 1);
poly_copy(tmp2, A, n), poly_copy(tmp3, B, m);
ntt(tmp2, len, 1), ntt(tmp3, len, 1);
for(int i=0;i<len;i++) C[i] = 1LL*tmp2[i]*tmp3[i]%MOD;
ntt(C, len, -1); poly_clear(C, n+m-1, len);
poly_clear(tmp2, 0, len); poly_clear(tmp3, 0, len);
}
void poly_int(int *A, int *B, int n) {
for(int i=n-1;i>=0;i--)
A[i+1] = 1LL*B[i]*inv[i+1]%MOD;
A[0] = 0;
}
void poly_deri(int *A, int *B, int n) {
for(int i=1;i<n;i++)
A[i-1] = 1LL*B[i]*i%MOD;
}
int tmp4[MAXN + 5];
void poly_ln(int *A, int *B, int n) {
poly_inv(A, B, n); poly_deri(tmp4, A, n);
poly_mul(tmp4, B, B, n, n);
poly_clear(B, n, 4*n); poly_int(B, B, n-1);
poly_clear(tmp4, 0, n);
}
int tmp5[MAXN + 5], tmp6[MAXN + 5];
void poly_exp(int *A, int *B, int n) {
if( n == 1 ) {
B[0] = 1;
return ;
}
int len; for(len = 1;len < (n<<1);len <<= 1);
poly_exp(A, B, (n + 1) >> 1);
poly_copy(tmp5, A, n); poly_ln(B, tmp6, n);
ntt(tmp5, len, 1), ntt(tmp6, len, 1), ntt(B, len, 1);
for(int i=0;i<len;i++) B[i] = 1LL*B[i]*(1 + MOD - tmp6[i] + tmp5[i])%MOD;
ntt(B, len, -1); poly_clear(B, n, len);
poly_clear(tmp5, 0, len), poly_clear(tmp6, 0, len);
}
}oper;
int f[MAXN + 5], g[MAXN + 5];
int main() {
init(); int n;
scanf("%d", &n);
for(int i=0;i<n;i++)
scanf("%d", &f[i]);
oper.poly_exp(f, g, n);
for(int i=0;i<n;i++)
printf("%d ", g[i]);
puts("");
}@4 - 算法应用@
例题：@COGS - 2189@ 帕秋莉的超级多项式
已知 \(F(x)\) 是个 n 次多项式，\(k\) 是一个整数，且已知：
\[G(x)\equiv(1+ \ln (1+\frac{1}{exp(\int \frac{1}{\sqrt{F(x)}})}))^k\mod x^n\]
求 \(G'(x)\)，先取模再求导。
这道题……我暂时没有找到哪里可以交，但是找网上的程序对拍可以拍上几千组……应该是 AC 了的。
solution：按题意模拟 qwq。
主程序里面要把多项式的长度扩展成 2 的整数次幂才能 AC，但是依照我的理解是不用扩展的，不知道为什么。
#include<cmath>
#include<cstdio>
#include<algorithm>
using namespace std;
const int G = 3;
const int MAXN = 1000000;
const int MOD = 998244353;
int pow_mod(int b, int p) {
int ret = 1;
while( p ) {
if( p & 1 ) ret = 1LL*ret*b%MOD;
b = 1LL*b*b%MOD;
p >>= 1;
}
return ret;
}
int inv[MAXN + 5];
void init() {
inv[1] = 1;
for(int i=2;i<=MAXN;i++)
inv[i] = 1LL*(MOD - MOD/i)*inv[MOD%i]%MOD;
}
struct Polynomial{
void poly_copy(int *A, int *B, int n) {
for(int i=0;i<n;i++)
A[i] = B[i];
}
void poly_clear(int *A, int l, int r) {
for(int i=l;i<r;i++)
A[i] = 0;
}
void ntt(int *A, int n, int type) {
for(int i=0,j=0;i<n;i++) {
if( i < j ) swap(A[i], A[j]);
for(int l=(n>>1);(j^=l)<l;l>>=1);
}
for(int s=2;s<=n;s<<=1) {
int t = (s>>1);
int u = (type == 1) ? pow_mod(G, (MOD-1)/s) : pow_mod(G, (MOD-1)-(MOD-1)/s);
for(int i=0;i<n;i+=s) {
for(int j=0,p=1;j<t;j++,p=1LL*p*u%MOD) {
int x = A[i+j], y = 1LL*A[i+j+t]*p%MOD;
A[i+j] = (x + y)%MOD, A[i+j+t] = (x + MOD - y)%MOD;
}
}
}
if( type == -1 ) {
for(int i=0;i<n;i++)
A[i] = 1LL*A[i]*inv[n]%MOD;
}
}
int tmp1[MAXN + 5];
void poly_inv(int *A, int *B, int n) {
if( n == 1 ) {
B[0] = pow_mod(A[0], MOD-2);
return ;
}
int len; for(len = 1;len < (n<<1);len <<= 1);
poly_inv(A, B, (n + 1) >> 1);
poly_copy(tmp1, A, n);
ntt(tmp1, len, 1); ntt(B, len, 1);
for(int i=0;i<len;i++) B[i] = 1LL*B[i]*(2 + MOD - 1LL*tmp1[i]*B[i]%MOD)%MOD;
ntt(B, len, -1); poly_clear(B, n, len);
poly_clear(tmp1, 0, len);
}
int tmp2[MAXN + 5], tmp3[MAXN + 5];
void poly_mul(int *A, int *B, int *C, int n, int m) {
int len; for(len = 1;len < (n+m-1);len <<= 1);
poly_copy(tmp2, A, n), poly_copy(tmp3, B, m);
ntt(tmp2, len, 1), ntt(tmp3, len, 1);
for(int i=0;i<len;i++) C[i] = 1LL*tmp2[i]*tmp3[i]%MOD;
ntt(C, len, -1); poly_clear(C, n+m-1, len);
poly_clear(tmp2, 0, len); poly_clear(tmp3, 0, len);
}
void poly_int(int *A, int *B, int n) {
for(int i=n-1;i>=0;i--)
A[i+1] = 1LL*B[i]*inv[i+1]%MOD;
A[0] = 0;
}
void poly_deri(int *A, int *B, int n) {
for(int i=1;i<n;i++)
A[i-1] = 1LL*B[i]*i%MOD;
}
int tmp4[MAXN + 5];
void poly_ln(int *A, int *B, int n) {
poly_inv(A, B, n); poly_deri(tmp4, A, n);
poly_mul(tmp4, B, B, n, n);
poly_clear(B, n, 4*n); poly_int(B, B, n-1);
poly_clear(tmp4, 0, n);
}
int tmp5[MAXN + 5], tmp6[MAXN + 5];
void poly_exp(int *A, int *B, int n) {
if( n == 1 ) {
B[0] = 1;
return ;
}
int len; for(len = 1;len < (n<<1);len <<= 1);
poly_exp(A, B, (n + 1) >> 1);
poly_copy(tmp5, A, n); poly_ln(B, tmp6, n);
ntt(tmp5, len, 1), ntt(tmp6, len, 1), ntt(B, len, 1);
for(int i=0;i<len;i++) B[i] = 1LL*B[i]*(1 + MOD - tmp6[i] + tmp5[i])%MOD;
ntt(B, len, -1); poly_clear(B, n, len);
poly_clear(tmp5, 0, len), poly_clear(tmp6, 0, len);
}
int tmp7[MAXN + 5];
void poly_pow(int *A, int *B, int n, int k) {
poly_ln(A, tmp7, n);
for(int i=0;i<n;i++) tmp7[i] = 1LL*k*tmp7[i]%MOD;
poly_exp(tmp7, B, n); poly_clear(tmp7, 0, n);
}
int tmp8[MAXN + 5], tmp9[MAXN + 5];
void poly_sqrt(int *A, int *B, int n) {
if( n == 1 ) {
B[0] = sqrt(A[0]);
return ;
}
int len; for(len = 1;len < (n<<1);len <<= 1);
poly_sqrt(A, B, (n + 1) >> 1);
poly_copy(tmp8, A, n); poly_inv(B, tmp9, n);
ntt(tmp8, len, 1), ntt(tmp9, len, 1), ntt(B, len, 1);
for(int i=0;i<len;i++) B[i] = 1LL*inv[2]*tmp9[i]%MOD*(tmp8[i] + 1LL*B[i]*B[i]%MOD)%MOD;
ntt(B, len, -1); poly_clear(B, n, len);
poly_clear(tmp8, 0, len), poly_clear(tmp9, 0, len);
}
}oper;
int f[MAXN + 5], g[MAXN + 5];
int main() {
init(); int n, k;
scanf("%d%d", &n, &k);
for(int i=0;i<n;i++)
scanf("%d", &f[i]);
int m; for(m = 1;m <= n;m <<= 1);
oper.poly_sqrt(f, g, m); oper.poly_copy(f, g, m); oper.poly_clear(g, 0, m);
oper.poly_inv(f, g, m); oper.poly_copy(f, g, m); oper.poly_clear(g, 0, m);
oper.poly_int(g, f, m); oper.poly_copy(f, g, m); oper.poly_clear(g, 0, m);
oper.poly_exp(f, g, m); oper.poly_copy(f, g, m); oper.poly_clear(g, 0, m);
oper.poly_inv(f, g, m); oper.poly_copy(f, g, m); oper.poly_clear(g, 0, m);
f[0]++;
oper.poly_ln(f, g, m); oper.poly_copy(f, g, m); oper.poly_clear(g, 0, m);
f[0]++;
oper.poly_pow(f, g, m, k); oper.poly_copy(f, g, m); oper.poly_clear(g, 0, m);
oper.poly_deri(g, f, n); oper.poly_copy(f, g, n); oper.poly_clear(g, 0, n);
for(int i=0;i<n;i++)
printf("%d", f[i]), (i == n-1) ? printf("\n") : printf(" ");
}
1727
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
