数学:多项式

参考文章:

博主的理解并不算深刻,有不恰当之处欢迎指出。


拉格朗日插值

通常用于解决如下问题:

  • 已知一个 n1 次函数 F(x) 图像上的 n 个点 (xi,yi),试求出这个函数 F(x)x=k 处的取值。

一般拉格朗日插值法

现在让我们来用一些比较 MO 的方式解决这个问题。考虑能不能用一些奇技淫巧构造这么一个函数。

考虑一个跟之前 CRT 比较类似的思路。对每个点 xi 构造一个函数,使得在 n 个点中,这个函数在 xi 上恰好取到 yi,而其它位置,即那些 jixj 上都恰好取到 0。将 n 个函数加起来即是我们要的答案多项式。

直接给出这个答案多项式的形式 F(k)=i=1nyijikxjxixj。和号去掉就是每个点对应的多项式。可以发现它确实在 k=xi 的时候取到 yik=xj 的时候分子乘积一定会有一项是 0,没了。

直接实现是 O(n2) 的,代码实现是简单的。

upd:拉插可以 O(n4) 反向根据点值求系数。感兴趣可以了解一下 CF1874E 这题。

xi 取值为连续正整数时

钦定 xi=i

那么原式变为 F(k)=i=1nyijikjij。其中 jikjij=ji(kj)ji(ij)=ji(kj)(i1)!(ni)!(i1)ni

预处理 j=1n(kj) 的前缀积、后缀积,阶乘逆元,用的时候直接前缀后缀乘起来中间空一个即可。O(n)

xi=i+c 时分母不变,分子预处理几乎没有影响。

重心拉格朗日插值法

这是什么,没用到过。

应用

i=1nikCF622F The Sum of the k-th Powers

(i=1nik)mod(109+7)

1n109,0k106

根据一些懂的都懂不懂的不懂的数学知识,这个玩意是一个关于 nk+1 次多项式。

代入 1,2,,k+2 把这个多项式在 n 处的取值求出来即可。因为自变量是连续正整数,可以使用上面的优化做到 O(n)

ll n,k,xx[N],yy[N];
ll fac[N],fnv[N],f[N],b[N];
ll qp(ll x,ll y) { return y?(y&1?x:1)*qp(x*x%P,y>>1)%P:1; }
ll lag(ll n,ll *x,ll *y,ll k) {
	fac[0]=1; for (ll i=1;i<=n;i++) fac[i]=fac[i-1]*i%P;
	fnv[n]=qp(fac[n],P-2); for (ll i=n-1;~i;i--) fnv[i]=fnv[i+1]*(i+1)%P;
	ll ret=0;
	f[0]=b[n+2]=1;
	for (ll i=1;i<=n+1;i++) f[i]=f[i-1]*(k-x[i])%P;
	for (ll i=n+1;i;i--) b[i]=b[i+1]*(k-x[i])%P;
	for (ll i=1;i<=n+1;i++) (ret+=y[i]*(f[i-1]*b[i+1]%P)%P*(fnv[i-1]*fnv[n+1-i]*(n+1-i&1?-1:1)%P))%=P;
	return (ret+P)%P;
}
void mian() {
	scanf("%lld%lld",&n,&k);
	for (ll i=1;i<=k+2;i++) xx[i]=i,yy[i]=(yy[i-1]+qp(i,k))%P;
	cout<<lag(k+1,xx,yy,n);
}

求只有一个自变量的次数已知的多项式的取值(P3270 [JLOI2016] 成绩比较

1rin1001ui109

跳过一堆步骤,总之就是中间有一步要求出每个 ti=l=0uilnri(nl)ri1,但是 ui 巨大,怎么办呢。

你惊喜地发现全程只有 l 在变,而 l 的最高次数恒为 n1,所以整个和式应该是一个关于 uin 次多项式,直接代入 n+1 个值把它在 ui 的取值求出来即可。

快速傅里叶变换(FFT)

已经不知道被这玩意劝退了多少次了。但理解之后其实不算很阴间的东西?

前置知识

多项式

我们平常使用的什么 f(x)=a0+a1x1+a2x2+ 什么的都属于系数表示法。用 a 这个系数序列就能表示整个多项式。

还有一种叫点值表示法的东西。因为一个 n 次多项式可以由 n+1 处已知的取值确定,所以用多项式上的 n+1 处取值来表示这个多项式。比如 (0,1),(1,2),(2,1) 表示 f(x)=x2+2x+1

点值表示法看上去非常反人类,但是你用系数表示法的话就反 FFT 了。

向量、三角函数、复数

每一个第一次学的人:?多项式是怎么跟复数扯上关系的?

简单来说就是实数域没有满足 FFT 需要的性质的数,所以就跑到复数域来乱搞了,因为多项式 f(x) 里的 x 也可以是复数。不要去尝试想象函数图像……

假设你有高中数学中向量、三角函数以及此方面的相关知识。大概就是向量的定义、相加;三角函数的定义;复数可表示为 z=a+bi,而以 a 为横坐标 b 为纵坐标定义一个叫复平面的东西的话,每一个复数都可以表示为原点到 (a,b) 的一条向量。

单位根

我们以复平面的原点画一个半径为 1 的圆,再将这个圆 n 等分,原点向每一个等分点连一条向量。那么幅角最小且不为 0 的那条向量代表的复数就是 n 次单位根(ωn)。由于复数的乘法是模长相乘,幅角相加,那么 ωn0,ωn1,,ωnn1 即为上面的 n 条向量。

把三角函数和向量相加的技术放到复平面上,有 wnk=cos2πkn+isin2πkn。又由几何意义易知 wnk=w2n2kwnk=wnk+n2。这两条在 FFT 里有大用。

快速傅里叶变换(FFT)

首先我们把听了一万年的多项式乘法之类的东西暂时忘掉,不要管它们。现在我们研究的是傅里叶变换,不是多项式乘法。在 OI 里你可以把它简单地理解为就是把一个多项式的系数表示法转为点值表示法。FFT 就是快速地实现这个过程。好像在别的领域的应用没有这么简单,但我们暂时用不到。

假设现在我们要求出一个 n1 次(n 项)多项式 F(x)=a0+a1x1+a2x2++an1xn1n 个点值,以进行点值表示。为方便处理,钦定 n2 的幂次方,不是就补几项。

考虑将多项式的奇偶项分开。

F(x)=(a0+a2x2++an2xn2)+(a1x1+a3x3++an1xn1)=(a0+a2x2++an2xn2)+x(a1+a3x2++an1xn2)

记两个新的多项式 F0(x)=a0+a2x++an2xn21F1(x)=a1+a3x++an1xn21

那么有 F(x)=F0(x2)+xF1(x2)

如果我们要求出一个点值,我们需要带入 x 为某值,求一遍 F(x),然后又要求一遍 F0F1。看起来还是没有前途。

但是如果我们把眼光放到复数域,事情就有点意思了。

对于每一个 0k<n2,我们代入单位根 ωnk。(前面提到的单位根的两条性质:wnk=w2n2kwnk=wnk+n2

F(ωnk)=F0(ωn2k)+ωnkF1(ωn2k)=F0(ωn2k)+ωnkF1(ωn2k)

再代入 ωnk+n2

F(ωnk+n2)=F0(ωn2(k+n2))+ωnk+n2F1(ωn2(k+n2))=F0(ωn2k+n)ωnkF1(ωn2k+n)=F0(ωn2k)ωnkF1(ωn2k)=F0(ωn2k)ωnkF1(ωn2k)

你发现,这时如果你算出了 F0(ωn2k)F1(ωn2k),你改一下正负性就可以一次性求出 F(x) 的两个点值 F(ωnk)F(ωnk+n2)

那我们考虑,令我们要求的 n 个点值就为 F(ωn0),F(ωn1),,F(ωnn1)

在求的过程中,我们枚举 0k<n2,将 F(ωnk),F(ωnk+n2) 分到一组,通过求出 F0(ωn2k)F1(ωn2k) 的值,将它们两个同时算出来。

那么整体来看,求出 F(ωn0),F(ωn1),,F(ωnn1),相当于需要求 F0(ωn20),F0(ωn21),,F0(ωn2n21)F1(ωn20),F1(ωn21),,F1(ωn2n21) 的值。把这两部分东西求出来后,按上面的方法组合一下就可以得到所有点值了。好像这个组合的方式叫蝴蝶操作。

你发现其实这两部分跟现在求解的问题没有任何区别。对于前一部分,只要令 n 变为 n2F 变为 F0,就是跟现在一模一样的问题。后一部分,让 F 变为 F1。每次,我们可以将现在的问题分解为两个和现在本质相同,但是规模减半的子问题。跟线段树建树有几分神韵相似。

一直这样对半劈下去,没过多久就会进入 n=1 的子问题,这时只有一个常数项,直接拿出来就是这个子问题的答案。然后把这个子问题的答案一层层带回去,一路把更大的子问题的答案也求出来即可。于是对一个 n 次多项式进行 FFT 的时间复杂度就是 O(nlogn) 的。

总之就是利用单位根的性质实现分治求解。非常绝妙!

快速傅里叶逆变换(IFFT)

因为点值表示法实在太过反人类,我们需要一个把多项式的点值表示法转成系数表示法的东西,即 IFFT。因为实现方式实在太过相似,通常不跟 FFT 作区分。

假设我们已知多项式 F(x)=a0+a1x1++an1xn1 在 FFT 之后得到的 n 个点值 b0=F(ωn0),b1=F(ωn1),bn1=F(ωnn1),那么有结论:多项式 G(x)=b0+b1x1++bn1xn1a0=G(ωn0)n,a1=G(ωn1)n,,an1=G(ωn(n1))n。即 G(x)n 个点值 G(ωn0),G(ωn1),,G(ωn(n1)) 在除以 n 之后就为 a0,a1,,an1

也就是我们把点值作为多项式 G(x) 的系数,把单位根 ωni 换为其倒数(将虚部取负)跑一遍 FFT,最后把得到的点值除以 n 就还原出系数了。

证明有点离谱,晚点补(

代码实现

因为涉及到复数运算,建议写一个复数类,或者直接用 STL 的 complex。

这就不可避免地给 FFT 带来了一定的常数和精度问题。

递归实现

即直接模拟上面那个东西。FFT 的过程本质上就是一个递归分治,所以比较容易按照这个思路实现。

但是递归实现的常数放到 FFT 上简直是雪上加霜,所以我们考虑把递归去掉。

迭代实现

也就是说我们需要一种方法直接知道 logn 层之后到 n=1 的状态时,每一个系数都在哪个地方。然后从最底层倒回去,一层一层地处理每个子问题的变化。

先跑一遍?可以,但不太优雅,也有点慢。

盗你谷第一篇题解的图。

你发现把原序列的二进制位翻转之后就得到了新序列的位置。其实可以理解,因为我们在第 i 层决策每个位置的数应该被放到哪一边的子问题的时候,看的就是这个原位置的从低到高数的第 i 个二进制位是不是 1,而这样带来的后果就是把它在新序列中的位置的从高到低数的第 i 个二进制位变成了这个值。这样有没有办法快速算出来呢?

直接给出递推式:设 r 是新序列每一项在原序列中的位置,那么 r[i]=(r[i>>1]>>1)|((i&1)<<(m-1))m 是层数。理解很简单,就是 r[i]=r[i>>1]*2+(i%2),用位运算技巧变成 r[i]=(r[i>>1]<<1)|(i&1),然后把二进制位反转了一下。

底层常数项搞定了。而对于从最底层往最上层处理的过程,回忆一下之前的式子 F(ωnk)=F0(ωn2k)+ωnkF1(ωn2k)F(ωnk+n2)=F0(ωn2k)ωnkF1(ωn2k),直接把两半子问题的元素依次拿出来乘系数加减组合一下放回原处即可。

注意如果是 IFFT 的话最后最好四舍五入一下,FFT 的精度问题还是不可忽视的。

下面代码中,com 是手写的复数类,ini 是求出第一个大于 a 的初始项数 _n2 的幂次方作为 FFT 点值数,并求出位逆序,mid 是此时子问题的规模的一半,w 是单位根,x 是单位根的幂,c 用于区分 FFT 和 IFFT,1 为 FFT,1 为 IFFT。

struct com {
	double a,b;
	com(double _a=0,double _b=0) { a=_a,b=_b; }
	friend com operator + (com x,com y) { return {x.a+y.a,x.b+y.b}; }
	friend com operator - (com x,com y) { return {x.a-y.a,x.b-y.b}; }
	friend com operator * (com x,com y) { return {x.a*y.a-x.b*y.b,x.a*y.b+x.b*y.a}; }
} a[N],b[N];
ll n=1,m=0,r[N];
void ini(ll _n) {
	while (n<=_n) n<<=1,m++;
	for (ll i=0;i<n;i++) r[i]=r[i>>1]>>1|(i&1)<<m-1;
}
void fft(com *a,ll c) {
	for (ll i=0;i<n;i++) if (i<r[i]) swap(a[i],a[r[i]]);
	for (ll mid=1;mid<n;mid<<=1) {
		com w(cos(pi/mid),c*sin(pi/mid));
		for (int i=0;i<n;i+=mid<<1) {
			com x(1,0);
			for (int j=0;j<mid;j++,x=x*w) {
				com y=a[i+j],z=x*a[i+mid+j];
				a[i+j]=y+z,a[i+mid+j]=y-z;
			}
		}
	}
	if (c==-1) for (ll i=0;i<n;i++) a[i].x=(ll)(a[i].x/n+0.5);
}

快速数论变换(NTT/FNTT)

明显 FFT 常数和精度问题都很大。

众所周知精度的最高境界是取模,所以 NTT 来了。即在模意义下快速实现系数表示法和点值表示法之间的转换。

其实准确来说这玩意应该叫 FNTT,NTT 是数论变换,但是没人在意。

前置知识:原根

若一个数 g 在模 P 意义下的阶(即 gimodP 的不同取值个数)为 ϕ(P),那么称 gP 的原根。

P 为质数时,可以理解为 g 在模 P 时的一个周期覆盖了 [1,P1] 里的所有元素。

快速数论变换(NTT)

原根这玩意逆天的地方在于,在模意义下,FFT 时用到的单位根的性质(ωnk=ω2n2kωnk=ωnk+n2 等)它都满足。证明待补。

于是,若 gP 的原根,ωngP1n(modP)

也就是说,NTT 就是把 FFT 的单位根换成上面的原根幂。从此我们就非常轻松地学会了在取模的条件下飞快地进行系数和点值表示之间的转换!

……吗?注意这个指数。因为 n 在模 φ(P) 意义下不一定有逆元,所以我们很可能找不到一个同余的 P1n,只能求出这个玩意的具体值了,那要是具体值不是整数又要怎么办呢……

那就不办了,考虑 P1n 什么时候才是整数。因为计算时 n 在我们的调整下变成了一个 2 的幂次方,所以如果 P=p2k+1 而且 n2k 的话,P1n 就是整数。这时直接除一下算出 P1n 的值然后幂一下求出单位根即可。

那就没有多少个质数模数满足条件了。一个著名模数是 998244353=119×223+1,它的原根是 3

当然人类的智慧是无穷的,确实有一个叫作任意模数 NTT 的东西,在下面。

实现

跟 FFT 差不多。

ll n=1,m=0,a[N],r[N];
void ini(ll _n) {
	while (n<=_n) n<<=1,m++;
	for (ll i=0;i<n;i++) r[i]=r[i>>1]>>1|(i&1)<<m-1;
}
void ntt(ll c) {
	for (ll i=0;i<n;i++) if (i<r[i]) swap(a[i],a[r[i]]);
	for (ll o=1;o<n;o<<=1) {
		ll w=qp(c==1?3:qp(3),(P-1)/(o<<1));
		for (ll j=0,x,y;j<n;j+=o<<1) for (ll k=j,v=1;k<j+o;k++,v=v*w%P)
			x=a[k],y=v*a[k|o]%P,a[k]=(x+y)%P,a[k|o]=(x+P-y)%P;
	}
	if (c==-1) for (ll i=0,iv=qp(n);i<n;i++) a[i]=a[i]*iv%P;
}

任意模数 NTT

见下方应用。

应用

多项式相关操作

乘法:P3803 【模板】多项式乘法(FFT)

给定一个 n 次多项式 F(x),和一个 m 次多项式 G(x)。请求出 F(x)G(x) 的卷积。

1n,m106,输入系数 [0,9],2s,500MB。

点值表示法简直是浑然天成地适合用来做卷积,不难知道,拿足够的点数,然后把两个多项式在这些点上的点值乘起来即是卷积后的对应点值。

所以我们先令 FFT 出的点值个数为大于 n+m1 的最小的 2 的幂次方,求出 F(x),G(x) 的这么多个点值,然后把点值相乘,再 IFFT 回去即是乘积结果了。钦定 n,m 同阶,复杂度 O(nlogn)

因为系数实在太小,这题把 998244353 当模数跑三遍 NTT 也是可以的。好像有个什么三次变两次优化,但博主还不会。

代码实现在后面。

任意模数乘法:P4245 任意模数多项式乘法

给定两个多项式 F(x),G(x),求 F(x)G(x),系数对 P 取模。不保证 P 是 NTT 模数。

1n,m1052P109+9,2s,500MB。

注意到卷积后得到的系数在不取模的结果下最大值也就 nP21023。那么就有一个很暴力的思路:拿出三个 NTT 模数分别跑一遍,然后用 CRT(中国剩余定理)合出最终输入的模数模意义下的结果。

这首先要求我们记住三个 NTT 模数。一般使用 469762049=7×226+1,998244353=119×223+1,1004535809=479×221+1,这三位原根都是 3。它们的 lcm1026 左右,可以轻松覆盖所有可能的系数。

最后就是 CRT 的技术了。这是简单的。同样是单次 O(nlogn),但是三个模数不可避免地带来了超级巨大的常数。

求逆:P4238 【模板】多项式乘法逆

给定一个 n 次多项式 F(x),求一个多项式 G(x) 使得 F(x)G(x)1(modxn),系数对 998244353 取模。

1n105,1s,125MB。

也是一个不知道前人是怎么想出来的东西。同样令 n2 的幂次方。

GF 在模 xn2 意义下的逆元,假设我们已经求出来了这个东西。

那么有

FG1(modxn2)FG1(modxn)FG1(modxn2)F(GG)0(modxn2)GG0(modxn2)

然后骚操作来了。两边平方得

(GG)0(modxn)G22GG+G20(modxn)

两边同乘 F

G2G+FG20(modxn)G2GFG2(modxn)

那么从 n=1 的一个常数项的逆元开始一层层迭代回去即可。NTT 优化乘法下可以做到单次求逆 O(nlogn)

注意不要在算对 xi 取模的结果的时候把 Fi 位之后的值也拿进去算,会假成双 log

lnP4725 【模板】多项式对数函数(多项式 ln)

给定一个 n 次多项式 F(x),求一个多项式 G(x) 使得 lnF(x)G(x)(modxn),系数对 998244353 取模。

1n105,2s,125MB。

expP4726 【模板】多项式指数函数(多项式 exp)

给定一个 n 次多项式 F(x),求一个多项式 G(x) 使得 eF(x)G(x)(modxn),系数对 998244353 取模。

1n105,2s,125MB。

快速幂:P5245 【模板】多项式快速幂

给定一个 n 次多项式 F(x)k,求一个多项式 G(x) 使得 (F(x))kG(x)(modxn),系数对 998244353 取模。

1n1050<k10105,1.5s,250MB。

多项式封装

一种来自 czk 的比较符合直觉而且兼容性还可以的封装写法,重载了括号,跟用普通数组没什么区别,常数还行。

目前只封了 NTT 和加减乘。

namespace poly {
	vector<ll> r;
	struct ply {
		ll n; vector<ll> a;
		void sz(ll x=0) {
			if (x) n=x;
			else while (n>1&&!a[n-1]) n--;
			a.resize(n);
		}
		ll& operator [] (ll x) {
			if (x>=n) sz(x+1);
			return a[x];
		}
		ply(ll x=0,ll y=0) { sz(x+1),a[x]=y; }
		void ntt(ll c) {
			for (ll i=0;i<n;i++) if (i<r[i]) swap(a[i],a[r[i]]);
			for (ll o=1;o<n;o<<=1) {
				ll w=qp(c==1?3:qp(3),(P-1)/(o<<1));
				for (ll j=0,x,y;j<n;j+=o<<1) for (ll k=j,v=1;k<j+o;k++,v=v*w%P)
					x=a[k],y=v*a[k|o]%P,a[k]=(x+y)%P,a[k|o]=(x+P-y)%P;
			}
			if (c==-1) for (ll i=0,iv=qp(n);i<n;i++) a[i]=a[i]*iv%P;
		}
		friend ply operator * (ply a,ply b) {
			ply ret; ll t=1,c=0;
			while (t<a.n+b.n-1) t<<=1,c++;
			a.sz(t),b.sz(t),ret.sz(t),r.resize(t);
			for (ll i=0;i<t;i++) r[i]=r[i>>1]>>1|(i&1)<<c-1;
			a.ntt(1),b.ntt(1);
			for (ll i=0;i<t;i++) ret[i]=a[i]*b[i]%P;
			ret.ntt(-1);
			return ret;
		}
		friend ply operator + (ply a,ply b) {
			ply ret(max(a.n-1,b.n-1));
			for (ll i=0;i<ret.n;i++) ret[i]=(a[i]+b[i])%P;
			return ret;
		}
		ply operator - () { ply ret; for (ll i=0;i<n;i++) ret[i]=P-a[i]; return ret; }
		friend ply operator - (ply a,ply b) { return a+(-b); }
	} a,b,c;
} using namespace poly;
ll n,m;
void mian() {
	scanf("%lld%lld",&n,&m);
	for (ll i=0;i<=n;i++) scanf("%lld",&a[i]);
	for (ll i=0;i<=m;i++) scanf("%lld",&b[i]);
	c=a*b;
	for (ll i=0;i<=n+m;i++) cout<<c[i]<<" ";
}
posted @   CarroT1212  阅读(34)  评论(4编辑  收藏  举报
相关博文:
阅读排行:
· 物流快递公司核心技术能力-地址解析分单基础技术分享
· 单线程的Redis速度为什么快?
· 展开说说关于C#中ORM框架的用法!
· Pantheons:用 TypeScript 打造主流大模型对话的一站式集成库
· SQL Server 2025 AI相关能力初探
点击右上角即可分享
微信分享提示