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

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

前置芝士

快速数论变换NTT

约定NTT前对于一定长度的范围处理和rev数组初始化函数为 $g e t r e v ()$ 。

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 (m o d x^{n})$ 。系数对 $998244353$ 取模。

考虑倍增法，假设已经有一个函数 $H (x), F (x) * H (x) \equiv 1 (m o d x^{⌈ \frac{n}{2} ⌉})$ 。我们又知道， $F (x) * G (x) \equiv 1 (m o d x^{⌈ \frac{n}{2} ⌉})$ 。

作差得：

F (x) * (G (x) - H (x)) \equiv 0 (m o d x^{⌈ \frac{n}{2} ⌉})

假定F不为0，则有：

G (x) - H (x) \equiv 0 (m o d x^{⌈ \frac{n}{2} ⌉})

平方得：

G^{2} (x) - 2 G (x) H (x) + H^{2} (x) \equiv 0 (m o d x^{n})

同时乘上F(x)：

G (x) - 2 H (x) + H^{2} (x) F (x) \equiv 0 (m o d x^{n})

G (x) \equiv 2 H (x) + F (x) H^{2} (x) (m o d x^{n})

递归进行计算即可，边界值为 $G (0) = F (0)^{M O D - 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{1}{x}) x^{n}

用这个转化式子：

F (x) = Q (x) * G (x) + R (x)

F (\frac{1}{x}) = Q (\frac{1}{x}) * G (\frac{1}{x}) + R (\frac{1}{x})

x^{n} F (\frac{1}{x}) = x^{n - m} Q (\frac{1}{x}) * x^{m} G (\frac{1}{x}) + x^{n - m + 1} * x^{m - 1} R (\frac{1}{x})

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)} (m o d 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

事实上多项式求 $l n$ 值基本是不会用到的（哪个闲人会让你求多项式的ln值），但是在例如快速幂这些算法中会用到，所以有一定价值。

给定 $F (x)$ ，求 $l n F (x) m o d x^{n}$ ，保证 $f_{0} = 1$ 。

有定理：多项式 $F (x)$ 有对数多项式当且仅当 $f_{0} = 1$ （证明请读者自行了解）。

设 $T (x) = l n F (x)$ ，我们发现多项式很好求导，而 $(l n x)^{'} = \frac{1}{x}$ ，很好计算，考虑求导：

T^{'} (x) = F^{'} (x) * \frac{1}{F (x)}

（复合函数求导）

发现这个只需要求逆就好了，容易计算，再将函数做不定积分就可以求出答案。

求导方法：

设一个多项式为 $\sum_{i = 0}^{n - 1} f_{i} x^{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 m o d x^{n}$ ，则 $F (H (x)) \equiv 0 m o d x^{⌈ \frac{n}{2} ⌉}$ 。

考虑将这个代入exp，我们发现：

G (x) \equiv e^{F (x)} m o d x^{n}

则：

l n G (x) - F (x) \equiv 0 m o d x^{n}

假设 $T (G (x)) = l n G (x) - F (x)$ ，则 $T (G (x)) \equiv 0 m o d x^{n}$ 。

对这个函数求导：

T^{'} (G (x)) = \frac{1}{G (x)}

因为：

(f - g)^{'} = f^{'} - g^{'}

(l n x)^{'} = \frac{1}{x}

c^{'} (常 数) = 0

在这里，由于已知，我们将 $F (x)$ 看作一个“常数”

代入：

G (x) \equiv H (x) - \frac{l n H (x) - F (x)}{\frac{1}{H (x)}} m o d x^{n}

G (x) \equiv H (x) (1 - l n H (x) - F (x)) m o d 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 l n F (x)}$

套 $l n$ 和 $e x p$ 即可。

有结论：次数 $k$ 可以对模数取模而没有影响。

如果 $a_{0} \neq 1$ ，就没有办法求 $l n$ 了。这时我们直接将整个多项式乘上 $a_{0}$ 的逆元，快速幂后再乘上原来 $a_{0}$ 值的 $k m o d ϕ (M O D)$ 倍即可。

如果 $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^{⌈ \frac{n}{2} ⌉}$ 。

G^{2} (x) \equiv F (x) \mod x^{⌈ \frac{n}{2} ⌉}

G^{2} (x) - H^{2} (x) \equiv 0 \mod x^{⌈ \frac{n}{2} ⌉}

G^{4} (x) - 2 G^{2} (x) H^{2} (x) + H^{4} (x) \equiv 0 \mod x^{n}

G^{4} (x) + 2 G^{2} (x) H^{2} (x) + H^{4} (x) \equiv 4 G^{2} (x) H^{2} (x) \mod x^{n}

G^{2} (x) + H^{2} (x) \equiv 2 G (x) H (x) \mod x^{n}

F (x) + H^{2} (x) \equiv 2 G (x) H (x) \mod x^{n}

G (x) \equiv \frac{F (x) + H^{2} (x)}{2 H (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}^{i} f_{i - j} g_{j}$ 。边界为 $f_{0} = 1$ 。

这个 $F$ 要依赖前项进行运算，不能直接卷积，我们考虑类cdq分治的方法。

假设我们已经计算了 $f_{l, m i d}$ ，考虑 $f_{l, m i d}$ 对 $f_{m i d + 1, r}$ 的贡献。考虑取 $g$ 的多少项， $f_{l} + g_{r - l} \to f_{r}$ 。所以我们复制 $g_{0}$ 到 $g_{r - l}$ ，将 $f_{l, m i d}$ 和 $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 阅读(39) 评论(0) 编辑收藏举报

刷新页面返回顶部

登录后才能查看或发表评论，立即登录或者逛逛博客园首页

相关博文：

· 2023.2.9【数学】快速傅里叶变换

· 【数学】拉格朗日插值

· 多项式基本操作

· 多项式笔记

· 【瞎口胡】多项式操作

阅读排行：
· 地球OL攻略 —— 某应届生求职总结
· 周边上新：园子的第一款马克杯温暖上架
· Open-Sora 2.0 重磅开源！
· 提示词工程——AI应用必不可少的技术
· .NET周刊【3月第1期 2025-03-02】

公告

友链

jeefy
hfjh
cyantist
ricky_lin
神犇Smallbasic

昵称： The_Last_Candy
园龄： 2年2个月
粉丝： 6
关注： 6

+加关注

2025年3月

日

一

二

三

四

五

六

The_Last_Candy

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

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

前置芝士

多项式求逆

多项式除法/多项式取模

多项式ln

多项式exp

多项式快速幂

多项式开根

任意模数：三模数NTT

分治NTT

公告

友链

搜索

常用链接

随笔分类

随笔档案

相册

阅读排行榜

评论排行榜

推荐排行榜

最新评论