【数学】简单的多项式技巧汇总

【数学】简单的多项式技巧汇总

下面对一些多项式常见操作进行总结

前置芝士

快速数论变换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) - 2H(x) + H^2(x)F(x) \equiv 0\ (mod\ x^n) \]

\[G(x) \equiv 2H(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) * G(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

首先要从一点说起——牛顿迭代法。

牛顿迭代常用于求一个函数的零点，假定我们有一个曲线：

image

我们任取一个点，做函数的切线，取其与\(x\)轴的交点，然后再以这个点的横坐标为\(x\)，继续以上过程，就可以快速接近一个零点。

image

假设切点为\(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);
}

完结撒花

posted @ 2023-07-30 21:02  The_Last_Candy  阅读(27)  评论(0)  编辑  收藏  举报