【数学】简单的多项式技巧汇总
下面对一些多项式常见操作进行总结
前置芝士
快速数论变换NTT
约定NTT前对于一定长度的范围处理和rev数组初始化函数为\(getrev()\)。
inline void getrev(int len)
{
tt = 1,tw = 0;
while(tt <= len) tt <<= 1,tw++;
rev[0] = 0;
for(int i = 1;i < tt;i++) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (tw - 1));
}
多项式求逆
给定\(n - 1\)次多项式\(F(x)\)，求\(G(x)\)使得\(F(x) * G(x) \equiv 1 \ (mod \ x^n)\)。系数对\(998244353\)取模。
考虑倍增法，假设已经有一个函数\(H(x),F(x) * H(x) \equiv 1\ (mod\ x^{\lceil \frac n2\rceil})\)。我们又知道，\(F(x) * G(x) \equiv 1\ (mod \ x^{\lceil \frac n2 \rceil})\)。
作差得：
\[F(x) * (G(x) - H(x)) \equiv 0\ (mod\ x^{\lceil \frac n2 \rceil}) \]假定F不为0，则有：
\[G(x) - H(x) \equiv 0\ (mod\ x^{\lceil \frac n2 \rceil}) \]平方得：
\[G^2(x) - 2G(x)H(x) + H^2(x) \equiv 0\ (mod\ x^n) \]同时乘上F(x)：
\[G(x) - 2F(x)H(x) + H^2(x)F(x) \equiv 0\ (mod\ x^n) \]\[G(x) \equiv 2F(x)H(x) + F(x)H^2(x)\ (mod\ x^n) \]递归进行计算即可，边界值为\(G(0) = F(0)^{MOD - 2}\)。
（代码中的\(x\)，\(y\)代表将\(x\)的逆元求到\(y\)里面，后面同理）
inline void getinv(int *x,int *y,int len)
{
if(len == 1)
{
y[0] = ksm(x[0],MOD - 2);
return;
}
getinv(x,y,(len + 1) >> 1);
getrev((len << 1) - 1);
for(int i = 0;i < len;i++) c[i] = x[i];
fill(c + len,c + tt,0);
NTT(c,1);NTT(y,1);
for(int i = 0;i < tt;i++) y[i] = (2ll * y[i] % MOD - 1ll * y[i] * y[i] % MOD * c[i] % MOD + MOD) % MOD;
NTT(y,-1);
fill(y + len,y + tt,0);
}
多项式除法/多项式取模
给定\(n\)次多项式\(F(x)\)和\(m\)次多项式\(G(x)\)，求多项式\(Q(x)\)和\(R(x)\)，使：
-
\(Q(x)\)为\(n - m\)次，\(R(x)\)小于\(m\)次。
-
\(F(x) = G(x) * Q(x) + R(x)\)
考虑\(A^R(x)\)为将A系数反转过后取的多项式，有：
\[A^R(x) = A(\frac 1x)x^n \]用这个转化式子：
\[F(x) = Q(x) * G(x) + R(x) \]\[F(\frac 1x) = Q(\frac 1x) * G(\frac 1x) + R(\frac 1x) \]\[x^nF(\frac 1x) = x^{n - m}Q(\frac 1x) * x^mG(\frac 1x) + x^{n - m + 1} * x^{m - 1}R(\frac 1x) \]\[F^R(x) = Q^R(x) * G^R(x) + x^{n - m + 1}R^R(x) \]\[Q^R(x)\equiv F^R(x) * \frac 1{G^R(x)}\ (mod\ x^{n - m + 1}) \]由于\(Q^R(x)\)的次数是\(n - m\)，所以这样刚好能够求出\(Q(x)\)，然后
\[R(x) = F(x) - Q(x) * R(x) \]分别计算即可。
inline void rvs(int *x,int len)
{
for(int i = 0;len - i - 1 > i;i++) swap(x[i],x[len - i - 1]);
}
inline void div(int *x,int len1,int *y,int len2,int *quo,int *rem)
{
for(int i = 0;i < len2;i++) GR[len2 - 1 - i] = y[i];
if(len1 - len2 + 1 <= len2) fill(GR + len1 - len2 + 1,GR + len2,0);
rvs(x,len1);
getinv(GR,c,len1 - len2 + 1);
getrev((len1 + len2) << 1);
for(int i = 0;i < len1;i++) quo[i] = x[i];
fill(quo + len1,quo + tt,0);
NTT(c,1);NTT(quo,1);
for(int i = 0;i < tt;i++) quo[i] = 1ll * quo[i] * c[i] % MOD;
NTT(quo,-1);NTT(c,-1);
fill(quo + len1 - len2 + 1,quo + tt,0);
rvs(quo,len1 - len2 + 1);
rvs(x,len1);
NTT(quo,1);NTT(y,1);
for(int i = 0;i < tt;i++) y[i] = 1ll * y[i] * quo[i] % MOD;
NTT(quo,-1);NTT(y,-1);
fill(quo + len1 - len2 + 1,quo + tt,0);
for(int i = 0;i < tt;i++) rem[i] = (x[i] - y[i] + MOD) % MOD;
fill(rem + len2 - 1,rem + tt,0);
}
多项式ln
事实上多项式求\(ln\)值基本是不会用到的（哪个闲人会让你求多项式的ln值），但是在例如快速幂这些算法中会用到，所以有一定价值。
给定\(F(x)\)，求\(lnF(x)\ mod\ x^n\)，保证\(f_0 = 1\)。
有定理：多项式\(F(x)\)有对数多项式当且仅当\(f_0 = 1\)（证明请读者自行了解）。
设\(T(x) = lnF(x)\)，我们发现多项式很好求导，而\((ln\ x)' = \frac 1x\)，很好计算，考虑求导：
\[T'(x) = F'(x) * \frac 1{F(x)} \]（复合函数求导）
发现这个只需要求逆就好了，容易计算，再将函数做不定积分就可以求出答案。
求导方法：
设一个多项式为\(\sum_{i = 0}^{n - 1}f_ix^i\)，则其导函数为\(\sum_{i = 0}^{n - 2}f_{i + 1}(i + 1)x^i\)。
不定积分反之。
inline void diff(int *x,int len)
{
for(int i = 0;i < len - 1;i++) x[i] = 1ll * x[i + 1] * (i + 1) % MOD;
x[len - 1] = 0;
}
inline void intg(int *x,int len)//len是多项式现在的长度
{
for(int i = len;i >= 1;i--) x[i] = 1ll * x[i - 1] * inv[i] % MOD;
x[0] = 0;
}
inline void getln(int *x,int *y,int len)
{
fill(y,y + len * 3,0);
fill(b,b + len * 3,0);
for(int i = 0;i < len;i++) y[i] = x[i];
getinv(y,b,len);
diff(y,len);
getrev(len << 1);
NTT(y,1);NTT(b,1);
for(int i = 0;i < tt;i++) y[i] = 1ll * y[i] * b[i] % MOD;
NTT(y,-1);
intg(y,len - 1);
fill(y + len,y + tt,0);
}
多项式exp
首先要从一点说起——牛顿迭代法。
牛顿迭代常用于求一个函数的零点，假定我们有一个曲线：
我们任取一个点，做函数的切线，取其与\(x\)轴的交点，然后再以这个点的横坐标为\(x\)，继续以上过程，就可以快速接近一个零点。
假设切点为\(x_0\)，与\(x\)轴交点为\(x_1\)(暂时不知道)，所以切线方程为：
\[y = f'(x_0)(x - x_0) + f(x_0) \]由于与\(x\)轴相交时，\(y = 0\)，解得
\[x = x_0 - \frac{f(x_0)}{f'(x_0)} \]对于多项式来说也一样：
试求\(G(x)\)，使\(F(G(x)) \equiv 0\)，设上次迭代的函数为\(H(x)\),则：
\[G(x) = H(x) - \frac{F(H(x))}{F'(H(x))} \]有结论是，这样迭代每次的精度是翻倍的，也就是说假设\(F(G(x)) \equiv 0\ mod\ x^n\)，则\(F(H(x))\equiv 0\ mod\ x^{\lceil{\frac n2}\rceil}\)。
考虑将这个代入exp，我们发现：
\[G(x) \equiv e^{F(x)}\ mod \ x^n \]则：
\[ln\ G(x) \ - \ F(x)\equiv 0\ mod \ x^n \]假设\(T(G(x)) = ln\ G(x) \ - \ F(x)\)，则\(T(G(x))\equiv 0\ mod\ x^n\)。
对这个函数求导：
\[T'(G(x)) = \frac 1{G(x)} \]因为：
\[(f - g)' = f' - g' \]\[(ln\ x)' = \frac 1x \]\[c'(常数) = 0 \]在这里，由于已知，我们将\(F(x)\)看作一个“常数”
代入：
\[G(x) \equiv H(x) - \frac{ln\ H(x)\ -\ F(x)} {\frac 1{H(x)}}\ mod\ x^n \]\[G(x) \equiv H(x)(1 - ln\ H(x) - F(x))\ mod\ x^n \]递归计算即可。
越复杂使用的辅助数组越多，一定要合理分配辅助数组。
inline void getexp(int *x,int *y,int len)
{
if(len == 1)
{
y[0] = 1;
return;
}
getexp(x,y,(len + 1) >> 1);
fill(e + 0,e + len,0);
getln(y,e,len);
for(int i = 0;i < len;i++) e[i] = (x[i] - e[i] + MOD) % MOD;
e[0]++;
fill(e + len,e + tt,0);
getrev(len << 1);
NTT(e,1);NTT(y,1);
for(int i = 0;i < tt;i++) y[i] = 1ll * e[i] * y[i] % MOD;
NTT(y,-1);
fill(y + len,y + tt,0);
}
多项式快速幂
比较简单：
\(F^k(x) = e^{k\ lnF(x)}\)
套\(ln\)和\(exp\)即可。
有结论：次数\(k\)可以对模数取模而没有影响。
如果\(a_0 \neq 1\)，就没有办法求\(ln\)了。这时我们直接将整个多项式乘上\(a_0\)的逆元，快速幂后再乘上原来\(a_0\)值的\(k\ mod\ \phi(MOD)\)倍即可。
如果\(a_0 = 0\)，就没有办法求逆元了，这时我们将多项式整体左移至第一位不为0为止，快速幂完再右移即可，注意特判次数乘位移数取模前是否大于\(n\)，如果是则全部设为0。
inline void polyksm(int *x,int *y,int pts,int pts2,int len)
{
int pos = 0,org;
while(!x[pos] && pos < len) ++pos;
if(1ll * pts * pos > len) return;
for(int i = 0;i < len - pos;i++) x[i] = x[i + pos];
fill(x + len - pos,x + len,0);
len -= pos;
org = x[0];
for(int i = 0;i < len;i++) x[i] = 1ll * x[i] * ksm(org,MOD - 2) % MOD;
fill(c,c + len * 3,0);
getln(x,c,len);
for(int i = 0;i < len;i++) c[i] = 1ll * c[i] * pts % MOD;
getexp(c,y,len);
for(int i = 0;i < len;i++) y[i] = 1ll * y[i] * ksm(org,pts2) % MOD;
len += pos;
for(int i = len - 1;i >= pos * pts;i--) y[i] = y[i - pos * pts];
fill(y,y + pos * pts,0);
}
（没有包含特判部分）。
多项式开根
求\(G(x)\)，使\(G^2(x) \equiv F(x) \mod \ x^n\)。
倍增。
设\(H^2(x)\equiv F(x)\mod\ x^{\lceil \frac n2 \rceil}\)。
\[G^2(x) \equiv F(x)\mod\ x^{\lceil \frac n2 \rceil} \]\[G^2(x) - H^2(x) \equiv 0\mod\ x^{\lceil \frac n2 \rceil} \]\[G^4(x) - 2G^2(x)H^2(x) + H^4(x) \equiv 0\mod\ x^n \]\[G^4(x) + 2G^2(x)H^2(x) + H^4(x) \equiv 4G^2(x)H^2(x)\mod x^n \]\[G^2(x) + H^2(x) \equiv 2G(x)H(x) \mod\ x^n \]\[F(x) + H^2(x) \equiv 2G(x)H(x)\mod x^n \]\[G(x)\equiv \frac {F(x) + H^2(x)}{2H(x)} \]递归计算即可，边界值\(G(0)\)用二次剩余计算。
inline void sqrt(int *x,int *y,int len)
{
if(len == 1)
{
y[0] = getsq(x[0]);
return;
}
sqrt(x,y,(len + 1) >> 1);
fill(c + 0,c + len * 2,0);
getinv(y,c,len);
int iv2 = ksm(2,MOD - 2);
getrev(len << 1);
for(int i = 0;i < len;i++) help[i] = x[i];
fill(help + len,help + tt,0);
NTT(help,1);NTT(c,1);
for(int i = 0;i < tt;i++) help[i] = 1ll * help[i] * c[i] % MOD;
NTT(help,-1);
for(int i = 0;i < len;i++) y[i] = 1ll * (y[i] + help[i]) % MOD * iv2 % MOD;
fill(y + len,y + tt,0);
}
任意模数：三模数NTT
对于任意模数取模，如果\(n \leq 10^5,a_i \leq 10^9\)，我们发现最终系数小于\(10^{23}\)。所以我们取3个便于计算的模数，并且相乘大于\(10^{23}\)，最后CRT合并，就可以得到正确答案。
证明：咕。
#include<bits/stdc++.h>
using namespace std;
const int N = 4e5 + 5,MOD[3] = {998244353,1004535809,469762049},g = 3,gi[3] = {332748118,334845270,156587350};
inline int ksm(int base,int pts,int i)
{
int ret = 1;
for(;pts > 0;pts >>= 1,base = 1ll * base * base % MOD[i])
if(pts & 1)
ret = 1ll * ret * base % MOD[i];
return ret;
}
inline void exgcd(__int128 a,__int128 b,__int128 &x,__int128 &y)
{
if(b == 0)
{
x = 1,y = 0;
return;
}
exgcd(b,a % b,x,y);
__int128 t = y;
y = x - (a / b) * y;
x = t;
}
int rev[N],tt,tw,F[N],G[N],R[N],a[N],b[N],ans[4][N],n,m,p;
inline void getrev(int len)
{
tt = 1,tw = 0;
while(tt < len) tt <<= 1,tw++;
for(int i = 1;i < tt;i++) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (tw - 1));
}
inline void NTT(int *x,int type,int mod)
{
for(int i = 0;i < tt;i++) if(i > rev[i]) swap(x[i],x[rev[i]]);
for(int mid = 1;mid < tt;mid <<= 1)
{
int omega = ksm((type == 1) ? g : gi[mod],(MOD[mod] - 1) / (mid << 1),mod);
for(int i = 0,R = mid << 1;i < tt;i += R)
{
int w = 1;
for(int j = i;j < i + mid;j++,w = 1ll * w * omega % MOD[mod])
{
int X = x[j],Y = 1ll * x[j + mid] * w % MOD[mod];
x[j] = (X + Y) % MOD[mod];
x[j + mid] = (X - Y + MOD[mod]) % MOD[mod];
}
}
}
if(type == -1)
{
int iv = ksm(tt,MOD[mod] - 2,mod);
for(int i = 0;i < tt;i++) x[i] = 1ll * x[i] * iv % MOD[mod];
}
}
inline int crt(int x,int y,int z)
{
__int128 M = MOD[0] * 1ll * MOD[1];
M *= MOD[2];
__int128 iv10,iv20,iv01,iv21,iv02,iv12;
iv10 = ksm(MOD[1],MOD[0] - 2,0);
iv20 = ksm(MOD[2],MOD[0] - 2,0);
iv01 = ksm(MOD[0],MOD[1] - 2,1);
iv21 = ksm(MOD[2],MOD[1] - 2,1);
iv02 = ksm(MOD[0],MOD[2] - 2,2);
iv12 = ksm(MOD[1],MOD[2] - 2,2);
__int128 ret = 0;
ret = (1ll * x * MOD[1] % M * MOD[2] % M * iv10 % M * iv20 % M + 1ll * y * MOD[0] % M * MOD[2] % M * iv01 % M * iv21 % M) % M;
ret = (ret + 1ll * z * MOD[0] % M * MOD[1] % M * iv02 % M * iv12 % M) % M;
return ret % p;
}
int main()
{
cin>>n>>m>>p;++n;++m;
for(int i = 0;i < n;i++) cin>>F[i];
for(int i = 0;i < m;i++) cin>>G[i];
getrev(m + n);
for(int i = 0;i < 3;i++)
{
memset(a,0,sizeof(a));memset(b,0,sizeof(b));memset(ans[i],0,sizeof(ans[i]));
for(int j = 0;j < n;j++) a[j] = F[j];
for(int j = 0;j < m;j++) b[j] = G[j];
NTT(a,1,i);NTT(b,1,i);
for(int j = 0;j < tt;j++) ans[i][j] = 1ll * a[j] * b[j] % MOD[i];
NTT(ans[i],-1,i);
}
for(int i = 0;i < n + m;i++) ans[3][i] = crt(ans[0][i],ans[1][i],ans[2][i]);
for(int i = 0;i < n + m - 1;i++) cout<<ans[3][i]<<" ";
return 0;
}
分治NTT
给定\(n\)次多项式\(G(x)\)，\(G(0) = 0\)，求\(F(x)\)。
其中\(f_i = \sum_{j = 1}^if_{i - j}g_j\)。边界为\(f_0 = 1\)。
这个\(F\)要依赖前项进行运算，不能直接卷积，我们考虑类cdq分治的方法。
假设我们已经计算了\(f_{l,mid}\)，考虑\(f_{l,mid}\)对\(f_{mid + 1,r}\)的贡献。考虑取\(g\)的多少项，\(f_l + g_{r - l} \to f_r\)。所以我们复制\(g_0\)到\(g_{r - l}\)，将\(f_{l,mid}\)和\(g_{0,r - l}\)做加法卷积贡献即可。
inline void cdqNTT(int *x,int *y,int l,int r)//[l,r)
{
if(l == r - 1)
{
if(!l) x[l] = 1;
return;
}
int mid = (l + r) >> 1;
cdqNTT(x,y,l,mid);
getrev(r - l);
for(int i = 0;i < mid - l;i++) a[i] = x[i + l];
fill(a + mid - l,a + tt,0);
for(int i = 0;i < r - l;i++) b[i] = y[i];
fill(b + r - l,b + tt,0);
NTT(a,1);NTT(b,1);
for(int i = 0;i < tt;i++) a[i] = 1ll * a[i] * b[i] % MOD;
NTT(a,-1);
for(int i = mid;i < r;i++) x[i] = (x[i] + a[i - l]) % MOD;
cdqNTT(x,y,mid,r);
}
完结撒花