多项式(Poly)笔记

开头先扔板子:多项式板子们

定义

多项式(polynomial)是形如 P(x)=i=0naixi 的代数表达式。其中 x 是一个不定元。

(P(x)) 称为这个多项式的次数。

多项式的基本运算

  • 多项式的加减法

A(x)=i=0naixi,B(x)=i=0mbixi

A(x)±B(x)=i=0max(n,m)(ai±bi)xi

  • 多项式的乘法

A(x)×B(x)=k=0nmi+j=kaibjxk

  • 多项式除法

这里讨论多项式的带余除法。

可以证明,一定存在唯一的 q(x),r(x) 满足 A(x)=q(x)B(x)+r(x),且 (r(x))<(B(x))

q(x) 称为 B(x)A(x) 的商式,r(x) 称为 B(x)A(x) 的余式。记作:

A(x)r(x)(mod B(x))

特别的,若 r(x)=0,则称 B(x) 整除 A(x)A(x)B(x) 的一个倍式,B(x)A(x) 的一个因式。记作 B(x)|A(x)

有关多项式的引理

  • 对于 n+1 个点可以唯一确定一个 n 次多项式。

  • 如果 f(x),g(x) 都是不超过 n 次的多项式,且它对于 n+1 个不同的数 α1,α2αn 有相同的值,即 i[1,n+1],iZ,f(αi)=g(αi)。则 f(x)=g(x)

多项式的点值表示

如果选取 n+1 个不同的数 x0,x1,xn 对多项式进行求值,得到 A(x0),A(x1)A(xn),那么就称 {(xi,A(xi)) | 0in,iZ}A(x) 的点值表示。

快速傅里叶变换(FFT)

快速傅里叶变换是借助单位根进行求值和插值,从而快速进行多项式乘法的算法。

单位根

将复平面上的单位圆均分成 n 份,从 x 轴数,第 i 条线与单位圆的交点称为 i 次单位根,记作 ωni

根据定义,可以得到:ωn1=isinα+cosα,α=2πn

根据欧拉恒等式,可以得到:ωn1=e2πin

由此那么可以得到下面的性质:

  • ωni×ωnj=ωni+j

  • ωni+n2=ωni

  • ωni+n=ωni

离散傅里叶变换(DFT)

离散傅里叶变换,是将 ωnk,0k<n 代入 f(x)g(x) 中求值的过程。

对于朴素的方法,每次代入一个单位根,然后用 O(n) 的复杂度计算函数值。时间复杂度 O(n2)

离散傅里叶变换利用了单位根的性质巧妙优化了这个过程。离散傅里叶变换过程如下:

首先将 f(x) 根据次数的奇偶性拆成两部分,分别分为:

{e(x)=i=0n22a2ix2i=a0+a2x2+a4x4an2xn2o(x)=i=0n22a2i+1x2i+1=a1x+a3x3+a5x5an1xn1

{e(x)=i=0n22a2ixi=a0+a2x+a4x2an2xn22o(x)=i=0n22a2i+1xi=a1+a3x1+a5x2an1xn22

f(x)=e(x2)+xo(x2)

ωnk 代入得到:

{f(ωnk)=e((ωnk)2)+ωnko((ωnk)2)=e(ωn2k)+ωnko(ωn2k)f(ωnk+n2)=e((ωnk+n2)2)+ωnk+n2o((ωnk+n2)2)=e(ωn2k)ωnko(ωn2k)

然后你发现,f(ωnk)f(ωnk+n2) 仅仅差了一个符号!!!

所以只需要求出 e(x)o(x)ωnk0kn2)上的取值,就可以推出 f(x) 在两倍点数上的取值。

每次问题规模缩小一半,因此时间复杂度 O(nlogn)

离散傅里叶逆变换(IDFT)

假设对于两个多项式都得到了 t=n+m1 个点值,设为 {(xi,A(xi)) | 0i<t,iZ},{(xi,B(xi)) | 0i<t,iZ}

那么可以知道,多项式 C(x)=A(x)×B(x) 的点值表示一定为:

{(xi,A(xi)B(xi)) | 0i<t,iZ}

现在,只需要将这 t 个点插值回去,就可以得到 A(x)B(x) 了。

先设这 t 个点值分别是:{(xi,vi) | 0i<t,iZ},设最后的多项式为 C(x)=i=0n+m2cixi,这里直接给出结论:

ck=1ni=0t1viωtki

由此可见,IDFT 和 DFT 仅仅有一个负号的区别。只要将所有的单位根从 ωnk 变成 ωnk 即可。

void FFT(cp a[], int n, int op) {
	if (n == 1) return;
	cp a1[n], a2[n];
	rop(i, 0, n >> 1) a1[i] = a[i << 1], a2[i] = a[(i << 1) + 1];
	FFT(a1, n >> 1, op), FFT(a2, n >> 1, op);
	cp Wn = {cos(2 * pi / n), op * sin(2 * pi / n)};
	cp Wk = {1, 0};
	rop(i, 0, n >> 1) {
		a[i] = a1[i] + Wk * a2[i];
		a[i + (n >> 1)] = a1[i] - Wk * a2[i];
		Wk = Wk * Wn;
	}
}
void FFT(cp a[], cp b[], int n, int m) {
	m = n + m; n = 1;
	while (n <= m) n <<= 1;
	FFT(a, n, 1); FFT(b, n, 1);
	rop(i, 0, n) a[i] = a[i] * b[i];
	FFT(a, n, -1);
	rep(i, 0, m) a[i].x = a[i].x / n;
}

FFT 优化

  • 三次变两次优化

原本的朴素 FFT,将 {a},{b} 两个序列分别求值,乘起来再 IDFT 插值一下,一共跑了三次 FFT。这是不好的。

三次变两次优化是这样的:将原序列合并成一个复数:{ai+bi×i}。做一遍 DFT 把求出的每个函数值平方。因为 (a+bi)2=(a2b2)+(2abi)。因此把虚部取出来以后除以 2 就是答案。

void FFT(cp a[], int n, int op) {
	if (n == 1) return;
	cp a1[n], a2[n];
	rop(i, 0, n >> 1) a1[i] = a[i << 1], a2[i] = a[(i << 1) + 1];
	FFT(a1, n >> 1, op), FFT(a2, n >> 1, op);
	cp Wn = {cos(2 * pi / n), op * sin(2 * pi / n)};
	cp Wk = {1, 0};
	rop(i, 0, n >> 1) {
		a[i] = a1[i] + a2[i] * Wk;
		a[i + (n >> 1)] = a1[i] - a2[i] * Wk;
		Wk = Wk * Wn;
	}
}
int main() {
	read(n, m);
	rep(i, 0, n) scanf("%lf", &a[i].x);
	rep(i, 0, m) scanf("%lf", &a[i].y);
	m = n + m; n = 1;
	while (n <= m) n <<= 1;
	FFT(a, n, 1);
	rop(i, 0, n) a[i] = a[i] * a[i];
	FFT(a, n, -1);
	rep(i, 0, m) printf("%d ", (int)(a[i].y / 2 / n + 0.5));
}
  • 蝴蝶变换优化

后面再补吧。其实本人感觉这个优化不是那么必要,因为三次变两次实在太快了。

FFT 例题

  1. A * B problem(plus)

可以设 x=10,把 a 写成 A(x)=i=0naixi 的形式(n=log10a)。同理可以把 b 转化为多项式 B(x)

这样求两个数相乘就是求 A(x)×B(x) 啊。

所以直接 O(nlogn) 做完了。

  1. P3338 [ZJOI2014] 力

给出 n 个数 q1,q2,qn,定义

Fj = i=1j1qi×qj(ij)2  i=j+1nqi×qj(ij)2

Ei = Fiqi

1in,求 Ei 的值。

首先发现这个除以 qi 就是没用。所以可以化简成:

Ej=i=1j1qi(ij)2  i=j+1nqi(ij)2

先看前面这个式子。答案就是:

(j1)2q1+(j2)2q2+(j3)2q3

f(x)=qixi,g(x)=1i2xi。可以发现,Ej=(f×g)[j]

再看后面这一块的式子。我们把 f(x) 的系数翻转,变成 f(x)=qnj+1xj。那么可以发现 Ej=(f×g)[nj+1]

跑两次 FFT 就完事了。

  1. P3723 [AH2017/HNOI2017] 礼物

首先发现加减相对于两个手环是对称的。因此可以把对一个手环的减法转化成对另一个手环的加法。这样可以假设全是在第一个手环上执行的加减操作。

第一个手环执行了加 c 的操作,且旋转过之后的序列为 [x1,x2xn],第二个手环为 [y1,y2yn]。计算差异值并化简,可以得到差异值是:

x2+y2+c2n+2c(xy)2xy

可以发现,这个序列只有最后一项是不定的。

因此将 y 序列翻转后再复制一倍,与 x 卷积,答案就是卷积后序列的 n+12n 项系数的 max

直接暴力枚举 c,加上前面依托就行了。

AC code

快速数论变换(NTT)

快速数论变换就是基于原根的快速傅里叶变换。

首先考虑快速傅里叶变换用到了单位根的什么性质。

  1. ωnk,0k<n 互不相同。

  2. ωnk=ωnk+n2

  3. ωnk=ω2n2k

数论中,原根恰好满足这些性质。

对于一个素数的原根 g,设 gn=gp1n。那么:

gnn=gp11(mod p)

gnn2=gp121(mod p)

gα+gβ=gα+β

gαnαk=gnk

我们发现它满足 ωnk 的全部性质!

因此,只需要用 gnk=gp1nk 带替全部的 ωnk 就可以了。

tips:对于一个质数,只有当它可以表示成 p=2α+1,且需要满足多项式的项数 <α 时才能使用 NTT。p 后面有个加一,是因为 gn 指数的分子上出现了 1p1 需要时二的整数次幂,是因为每次都要除以 2

bonus:常用质数及原根

p=998244353=119×223+1,g=3

p=1004535809=479×221+1,g=3

void NTT(int *a, int n, int op) {
	if (n == 1) return;
	int a1[n], a2[n];
	rop(i, 0, n >> 1) a1[i] = a[i << 1], a2[i] = a[(i << 1) + 1];
	NTT(a1, n >> 1, op), NTT(a2, n >> 1, op);
	int gn = qpow((op == 1 ? g : invg), (mod - 1) / n), gk = 1;
	rop(i, 0, n >> 1) {
		a[i] = (a1[i] + 1ll * gk * a2[i]) % mod;
		a[i + (n >> 1)] = (a1[i] - 1ll * gk * a2[i] % mod + mod) % mod;
		gk = 1ll * gk * gn % mod;
	}
}

NTT 优化:蝴蝶变换

盗大佬一张图

这是进行 NTT 的过程中数组下标的变化。

这样看似乎毫无规律。但是把他们写成二进制:

变换前:000 001 010 011 100 101 110 111

变换后:000 100 010 110 001 101 011 111

可以发现,就是对每个下标进行了二进制翻转。

因此可以先预处理出每个下标在递归底层对应的新下标。然后从底层往顶层迭代合并。

预处理每个下标的二进制翻转:

rop(i, 0, n) rev[i] = rev[i >> 1] >> 1 | (i & 1) << bit;

优化后的 NTT:

void NTT(int *a, int n, int op) {
	rop(i, 0, n) if (i < rev[i]) swap(a[i], a[rev[i]]);
	for (int mid = 1; (mid << 1) <= n; mid <<= 1) {
		int gn = qpow((op == 1 ? g : invg), (mod - 1) / (mid << 1));
		for (int i = 0, gk = 1; i < n; i += (mid << 1), gk = 1)
			for (int j = 0; j < mid; j ++ , gk = 1ll * gk * gn % mod) {
				int x = a[i + j], y = 1ll * a[i + j + mid] * gk % mod;
				a[i + j] = Mod(x + y), a[i + j + mid] = Mod(x - y);
			}
	}
}

当然了,FFT 也可以用蝴蝶变换来优化。实践证明,速度变快了至少二分之一。

FFT 的迭代实现

void FFT(cp *a, int n, int op) {
	rop(i, 0, n) if (i < rev[i]) swap(a[i], a[rev[i]]);
	for (int mid = 1; (mid << 1) <= n; mid <<= 1) {
		cp Wn = {cos(pi / mid), op * sin(pi / mid)};
		for (int i = 0; i < n; i += (mid << 1)) {
			cp Wk = {1, 0};
			for (int j = 0; j < mid; j ++ , Wk = Wk * Wn) {
				cp x = a[i + j], y = Wk * a[i + j + mid];
				a[i + j] = x + y, a[i + j + mid] = x - y;
			}
		}
	}
}

任意模数多项式乘法(MTT)

计算 f×g(mod p) 的结果(p 不一定能够表示成 2α+1 的形式)。

这个东西有两种做法,但是我只学会了三模 NTT。

首先,卷积之后每个系数最多达到 max{V}2×n 的大小。拿模板题举例,这个东西就是 1023。因此只需要找三个模数 p1,p2,p3,满足 p1p2p3>max{V}2×n,然后用他们分别 NTT,最后再 CRT / exCRT 合并即可。

int CRT(int a, int b, int c, int p) {
	int k = 1ll * Mod(b - a, p2) * inv1 % p2;
	LL x = 1ll * k * p1 + a;
	k = 1ll * Mod((c - x) % p3, p3) * inv2 % p3;
	x = (x + 1ll * k * p1 % p * p2 % p) % p;
	return x;
}
void MTT(int *a, int n, int *b, int m, int p) {
	int B = ((n + m) << 2) + 1;
	rev = new int [B]();
	int *a1 = new int [B](), *b1 = new int [B]();
	int *a2 = new int [B](), *b2 = new int [B]();
	int *a3 = new int [B](), *b3 = new int [B]();
	rop(i, 0, B) 
		a1[i] = a2[i] = a3[i] = a[i], b1[i] = b2[i] = b3[i] = b[i];
	NTT(a1, n, b1, m, p1); NTT(a2, n, b2, m, p2); NTT(a3, n, b3, m, p3);
	inv1 = inv(p1, p2); inv2 = inv(1ll * p1 * p2 % p3, p3);
	rep(i, 0, n + m) a[i] = CRT(a1[i], a2[i], a3[i], p);
}

这里选择的模数为:p1=998244353,p2=469762049,p3=1004535809。他们的原根都为 g=3

多项式求逆

求多项式 f(x) 的逆元 f1(x)f(x) 的逆元定义为满足 f(x)g(x)=1(mod xn) 的多项式 g(x)

使用倍增法即可求出多项式的逆元。

首先假设已经求出了 f(x)g(x)1(mod xn2)。再假设 f(x)mod xn 意义下逆元为 g(x)。那么有:

g(x)g(x)(mod xn2)

(g(x)g(x))0(mod xn2)

(g(x)g(x))20(mod xn)

g2(x)2g(x)g(x)+g2(x)0(mod xn)

两边同时乘以 f(x),可以得到:

g(x)2g(x)+f(x)g2(x)0(mod xn)

g(x)2g(x)f(x)g2(x)(mod xn)

g(x)g(x)(2f(x)g(x))(mod xn)

倍增求即可。

void Inv(int *f, int *g, int n) {
	if (n == 1) {
		g[0] = inv(f[0]); return;
	} Inv(f, g, (n + 1) >> 1);
	int m = 1, bit = 0;
	while (m < (n << 1)) m <<= 1, bit ++ ; bit -- ;
	rop(i, 0, n) tmp[i] = f[i];
	rop(i, n, m) tmp[i] = 0;
	rev = new int [m + 5]();
	rop(i, 1, m) rev[i] = (rev[i >> 1] >> 1) | (i & 1) << bit;
	NTT(tmp, m, 1); NTT(g, m, 1);
	rop(i, 0, m) g[i] = (2ll - 1ll * g[i] * tmp[i] % p + p) % p * g[i] % p;
	NTT(g, m, -1); int invn = inv(m);
	rop(i, 0, m) g[i] = 1ll * g[i] * invn % p;
	rop(i, n, m) g[i] = 0; 
	free(rev); rev = NULL;
}

分治 FFT

给定序列 g1n1,求序列 f0n1

其中 fi=j=1ifijgj,边界为 f0=1

答案对 998244353 取模。

其实这是个多项式求逆板子

首先考虑生成函数 F(x)=i=0+fixi,G(x)=i=0+gixi。然后可以发现:

F(x)G(x)=k=0+inftyxki+j=kfigj=F(x)f0x0

因此 F(x)(G(x)1)=f0,也就是:

F(x)f01G(x)(mod xn)

所以直接设 g0=0,然后把 1g 求逆就行了。

read(n);
rop(i, 1, n) read(g[i]);
rop(i, 1, n) g[i] = Poly::p - g[i];
(g[0] += 1) %= Poly::p; Poly::Inv(g, n);
rop(i, 0, n) write(' ', g[i]); return 0;

多项式对数函数(Polyln)

给定 f(x),求多项式 g(x) 使得 g(x)ln(f(x))(mod xn)

前置知识:简单的求导和积分,以及基本的多项式模板。

首先设 h(x)=ln(x),那么

g(x)h(f(x))(mod xn)

然后对同余方程两边同时求导,得到

g(x)[h(f(x))](mod xn)

根据复合函数求导公式 f(g(x))=f(g(x))g(x) 可得,

g(x)h(f(x))f(x)(mod xn)

然后又有 ln(x)=1x,因此

g(x)f(x)f(x)

计算 f(x)f1(x) 作为 a,b。计算 a×b(mod 998244353) 得到 g(x) ,然后求 g(x) 的积分就好了。

积分公式:xadx=1a+1xa+1

void der(int *f, int n) { rep(i, 1, n) f[i - 1] = 1ll * i * f[i] % p; f[n] = 0; }
void Int(int *f, int n) {dep(i, n, 1) f[i] = 1ll * inv(i) * f[i - 1] % p; f[0] = 0; }
void Ln(int *f, int n) {
	int B = (n << 2) + 5; int *_f = new int[B]();
	rep(i, 0, n) _f[i] = f[i];
	der(_f, n), Inv(f, n);
	NTT(f, n, _f, n); Int(f, n);
	free(_f); _f = NULL;
}

多项式指数函数(PolyExp)

给定多项式 f(x),求 g(x) 满足 g(x)ef(x)(mod xn)

前置知识:牛顿迭代。

牛顿迭代是用来求函数零点的有力工具。

例如,下面这个图是使用牛顿迭代法求 f(x)=x4+3x2+7x+3 零点的过程。

首先,随便选择一个点 (x0,f(x0)),过这个点做 f(x) 的切线。切线方程是 y=f(x0)(xx0)+f(x0)。它与 x 轴交于一点 A(0.56243,0)

接下来,再令 x1=0.56243,过点 (x1,f(x1)) 再做 f(x) 的切线,与 x 轴交于点 B(0.60088,0)

再令 x2=0.60088,如此迭代下去。可以发现,xn 会逐渐逼近零点。

刚才说切线方程为 y=f(x0)(xx0)+f(x0)。如果令 y=0,得到的 x 便是切线方程与 x 轴的交点,为

x=x0f(x0)f(x0)

运用于多项式,假设要求使 f(g(x))0(mod xn) 的多项式 g(x),并且已经知道 f(g0(x))0(mod xn2)

那么可以说,

g(x)=g0(x)f(g0(x))f(g0(x))

接下来解决多项式 Exp。所求为 g(x) 使得 g(x)ef(x)(mod xn)。两边都取 ln 得到:

lng(x)f(x)(mod xn)

移项得:

lng(x)f(x)0(mod xn)

F(g(x))=lng(x)f(x),那么所求就是 F(x) 的零点。

假设已经有 g0(x) 使得 F(g0(x))0(mod xn2),根据上面的牛顿迭代,有:

g(x)=g0(x)F(g0(x))F(g0(x))=g0(x)lng0(x)f(x)1g0(x)=g0(x)(1lng0(x)+f(x))

根据这个柿子倍增求即可。每次需要计算一个 ln,一个乘法,时间 O(nlogn)

写的有点丑,超级大肠数。

void Exp(int *f, int *g, int n) {
	if (n == 1) return void(g[0] = 1);
	Exp(f, g, (n + 1) >> 1);
	int B = (n << 2) + 5; int *lnb = new int[B]();
	rop(i, 0, n) lnb[i] = g[i]; Ln(lnb, n);
	tmp = new int[B](); rop(i, 0, n) tmp[i] = f[i];
	rop(i, 0, n) tmp[i] = (1ll * tmp[i] - lnb[i] + p) % p;
	tmp[0] ++ ; NTT(g, n, tmp, n);
	free(tmp); tmp = NULL; free(lnb); lnb = NULL;
}

多项式快速幂(PolyPower)

在模 xn 意义下计算 fk(x)

有了前面的知识铺垫,这部分就非常的简单。

根据对数恒等式,有 f(x)=elnf(x)

因此 fk(x)=eklnf(x)

f(x) 乘以 k,求多项式 ln,然后再 exp 回来就行了。

void Power(int *f, int n, int k) {
	Ln(f, n); rop(i, 0, n) f[i] = 1ll * f[i] * k % p; Exp(f, n);
}

多项式开根

在模 xn 意义下计算 f1k(x)

这个最简单。直接把 1kmodp,然后多项式快速幂。

posted @   Link-Cut-Y  阅读(316)  评论(0编辑  收藏  举报
相关博文:
阅读排行:
· 【.NET】调用本地 Deepseek 模型
· CSnakes vs Python.NET:高效嵌入与灵活互通的跨语言方案对比
· DeepSeek “源神”启动!「GitHub 热点速览」
· 我与微信审核的“相爱相杀”看个人小程序副业
· Plotly.NET 一个为 .NET 打造的强大开源交互式图表库
点击右上角即可分享
微信分享提示