多项式基础操作学习笔记

多项式简单操作学习笔记

前置知识与约定

• 多项式乘法,我们可以用FFT/NTT解决.
• 简单求导/微积分(简单的意思是如果这都不会就不要说自己学过了)
• ([P])表示艾弗森括号,当(P)为真时表达式的值为(1)，(P)为假时表达式的值为(0)
• (f^{(n)})表示函数(f(x))的(n)阶导函数

多项式乘法

我们可以用FFT/NTT在(mathcal O(n log n))的时间内解决, 但两者优劣势不尽相同.

FFT的缺点在于使用复数,精度差,使用单位根,常数大等.

而NTT的缺点在于必须在模数的原根较特殊的情况下才能使用(如(998244353)等).

Code

inline void NTT(int *f, int type) {
	for(int i = 0; i < n; ++i) if(i < rev[i]) std::swap(f[i], f[rev[i]]);
	for(int i = 1; i < n; i <<= 1) {
		int tk = quick_power(3, (mod - 1) / (i << 1));
		for(int j = 0; j < n; j += (i << 1)) {
			int t = 1, x, y;
			for(int k = 0; k < i; ++k, t = 1ll * t * tk % mod) {
				x = f[j + k], y = 1ll * t * f[j + k + i] % mod;
				f[j + k] = (x + y) % mod; f[j + k + i] = (x - y + mod) % mod;
			}
		}
	}
	if(type == 1) return ;
	int inv = quick_power(n, mod - 2);
	std::reverse(f + 1, f + n);
	for(int i = 0; i < n; ++i) f[i] = 1ll * f[i] * inv % mod;
}

多项式泰勒展开

如果函数(f(x))在(x_0)处存在(n)阶导数,那么定义(f(x))在(x = x_0)处的( m Taylor)展开式为:

[f(x) = sum_{i=0}^n frac {f^{(n)}(x_0)}{Gamma(i)} (x - x_0)^i + R_n(x) ]

其中(R_n(x))表示拉格朗日余项,即((x - x_0)^n)的高阶无穷小.

一个例子:(e^x = 1 + frac {x}{Gamma(1)} + frac {x^2}{Gamma(2)} + cdots)

多项式牛顿迭代

牛顿迭代可以用来求函数的零点,而多项式牛顿迭代可以求多项式函数的零点.

即对于一个多项式(F(x)), 求一个多项式(G_k(x))满足(F(G_k(x)) equiv 0 pmod{2^k})

假设我们已经求出了(G_{k-1}(x)),考虑如何推得(G_k(x))

我们把(F(G_k(x)))在(x = G_{k-1}(x))处进行( m Taylor)展开,可以得到:

[F(G_k(x)) = F(G_{k-1}(x)) + frac {F'(G_{k-1}(x))}{Gamma(1)} (G_k(x) - G_{k-1}(x)) ]

化简,移项得:

[G_k(x) = G_{k-1}(x) - frac {F(G_{k-1}(x))}{F'(G_{k-1}(x))} ]

多项式求逆

已知(F(x)),求(G(x))满足(F(x)G(x) equiv 1 pmod{x^n})

和推牛顿迭代的式子一样,我们假设已经求出了(G_0(x))满足:

[F(x)G_0(x) equiv 1 pmod{x^{lfloor frac {2}{n} floor}} ]

由题目条件可得:

[F(x)G(x) equiv 1 pmod{x^{lfloor frac {2}{n} floor}} ]

上下两式相减可得:

[F(x)G_0(x) - F(x)G(x) equiv 0 pmod{x^{lfloor frac {2}{n} floor}} ]

同余式两边同时除以(F(x))可得:

[G_0(x) - G(x) equiv 0 pmod{x^{lfloor frac {2}{n} floor}} ]

同余式两遍平方可得:

[G_0(x)^2 - 2G(x)G_0(x) + G(x)^2 equiv 0 pmod{x^n} ]

由题目条件,同余时两边同时乘(F(x))得:

[F(x)G_0(x)^2 - 2G_0(x) + G(x) equiv 0 pmod{x^n} ]

把式子移项,可得:

[G(x) = 2G_0(x) - F(x)G_0(x)^2 equiv 0 pmod{x^n} ]

递归求解即可.

Code

inline void PolyInv(int *f, int *g, int n) {
    static int c[MAXN];
    if(n == 1) {g[0] = quick_power(f[0], mod - 2); return ;}
    PolyInv(f, g, (n + 1) >> 1);
    int len = 1ll, lgm = 0ll; while(len < (n << 1)) {len <<= 1; ++lgm;}
    for(register int i = 1; i < len; ++i) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (lgm - 1));
    for(register int i = 0; i < n; ++i) c[i] = f[i];
    for(register int i = n; i < len; ++i) c[i] = 0; NTT(c, len, 1); NTT(g, len, 1);
    for(register int i = 0; i < len; ++i) g[i] = 1ll * (2 - c[i] * g[i] % mod + mod) % mod * g[i] % mod;
    NTT(g, len, -1); for(register int i = n; i < len; ++i) g[i] = 0;
}

多项式求导

没啥好说的吧, 给出公式:

[x^{a^{prime}} = ax^{a-1} ]

Code

inline void PolyDer(int *f, int *g, int len) {
	for(int i = 1; i < len; ++i) g[i - 1] = 1ll * i * f[i] % mod;
	g[len - 1] = 0;
}

多项式积分

根据牛顿-莱布尼兹公式:

[int x^a { m d}x = frac {1}{a+1} x^{a+1} ]

Code

inline void PolyInt(int *f, int *g, int len) {
	for(int i = 1; i < len; ++i) g[i] = 1ll * f[i - 1] * quick_power(i, mod - 2) % mod;
	g[0] = 0;
}

多项式对数函数

即多项式求(ln)

已知(F(x)),求(G(x))满足(G(x) equiv ln F(x) pmod{x^n})

设函数(H(x) = ln x),那么同余式右边即为(H(F(x))).

考虑将同余式两边求导,可得:

[G'(x) equiv H'(x) pmod{x^n} ]

即:

[G'(x) equiv frac{F'(x)}{F(x)} pmod{x^n} ]

那么我们把(F(x))求个导,再求个逆,然后把得到的两个多项式求一下卷积.

这样我们得到的是(G'(x)),然后再把他积回去即可.

Code

inline void PolyLn(int *f, int *g, int n) {
    static int a[MAXN], b[MAXN];
    PolyDer(f, a, n); PolyInv(f, b, n);
    int len = 1ll, lgm = 0ll; while(len < (n << 1)) {len <<= 1; ++lgm;}
    for(register int i = 1; i < len; ++i) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (lgm - 1));
    NTT(a, len, 1); NTT(b, len, 1);
    for(register int i = 0; i < len; ++i) a[i] = a[i] * b[i] % mod;
    NTT(a, len, -1); PolyInt(a, g, n);
    for(register int i = 0; i < len; ++i) a[i] = b[i] = 0;
}

多项式指数函数

即多项式求( m exp)

已知(F(x)),求(G(x))满足(G(x) equiv exp (F(x)) pmod{x^n}).

考虑同余式两边同时取自然对数:

[ln G(x) equiv F(x) pmod{x^n} ]

移项,可得:

[ln G(x) - F(x) equiv 0 pmod{x^n} ]

将这个式子套进 多项式牛顿迭代 可得:

[G(x) equiv G_0(x) - frac {ln G_0(x) - F(x)}{frac {1}{G_0(x)}} ]

化简可得:

[G(x) equiv G_0(x) - (ln G_0(x) - F(x)) imes G_0(x) pmod {x^n} ]

[G(x) equiv G_0(x) imes (1 - ln G_0(x) + F(x)) pmod {x^n} ]

递归+多项式Ln计算即可.

Code

inline void PolyExp(int *f, int *g, int n) {
    static int a[MAXN], b[MAXN];
    if(n == 1) {g[0] = 1ll; return ;} PolyExp(f, g, (n + 1) >> 1);
    int len = 1ll, lgm = 0ll; while(len < (n << 1)) {len <<= 1; ++lgm;}
    for(register int i = 1; i < len; ++i) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (lgm - 1));
    for(register int i = 0; i < (n << 1); ++i) a[i] = b[i] = 0;
    PolyLn(g, a, n); for(register int i = 0; i < n; ++i) b[i] = f[i];
    NTT(g, len, 1); NTT(a, len, 1); NTT(b, len, 1);
    for(register int i = 0; i < len; ++i) g[i] = (1ll - a[i] + b[i] + mod) % mod * g[i] % mod;
    NTT(g, len, -1); for(register int i = n; i < len; ++i) g[i] = 0;
}

多项式快速幂

已知(F(x)),求(G(x))满足(G(x) equiv F(x)^k pmod{x^n})

首先我们在初中阶段就应该学过这样一个式子:

[log_a b^k = k log_a b ]

这个式子对多项式也成立,

那么我们只需要把(F(x))求一个(ln),然后把它乘上(k),最后再(exp)回去即可.

Code

inline void PolyPow(int *f, int *g, int k, int len){
        static a[MAXN];
	PolyLn(f, a, len);
	for(int i = 0; i < len; ++i) a[i] = 1ll * a[i] * k % mod;
	PolyExp(a, g, len);
	for(int i = 0; i < len; ++i) a[i] = 0;
}

多项式开根

已知(F(x)),求(G(x))满足(G(x)^2 equiv F(x) pmod{x^n})

和多项式求逆类似,利用那个结论,得到:

[2 ln G(x) equiv ln F(x) pmod{x^n} ]

那么我们先把(F(x))求个(ln),然后把它的系数乘以(2)的逆元,最后再(exp)回去即可.

Code

inline void PolySqr(int *f, int *g, int len) {
        static a[MAXN];
	PolyLn(f, a, len);
	int invt = quick_power(2, (mod - 2));
	for(int i = 0; i < len; ++i) a[i] = 1ll * a[i] * invt % mod;
	PolyExp(a, g, len);
	for(int i = 0; i < len; ++i) a[i] = 0;
}

多项式除法 & 多项式取模

已知多项式(F(x),G(x)), 求(Q(x),R(x)), 满足(F(x) equiv Q(x)G(x) + R(x) pmod{x^n})

又是一种神仙思路。 考虑设(F_R(x))表示把(F(x))的系数 reverse 一遍的多项式。

显然对于(n)次多项式(F(x)), 我们有(F_R(x) = x^nF(frac 1x))。

然后我们来推式子：

[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}(frac 1x) cdot x^{m}G(frac 1x) + x^{n-m+1}x^{m-1}R(x) ]

[F_R(x) = Q_R(x)G_R(x) + x^{n-m+1}R_R(x) ]

然后把他换成(pmod {x^{n-m+1}})意义下的式子

[Q_R(x) equiv F_R(x)G_R^{-1}(x) pmod{x^{n-m+1}} ]

通过这个式子求出(Q_R(x))，然后代回上面的结果算(R(x))即可。

Code

inline void PolyDiv(int *f, int *g, int *q, int *r, int n, int m) {
    static a[MAXN], b[MAXN], c[MAXN];
    for(int i = 0; i < n; ++i) a[i] = f[n - i - 1];
    for(int i = 0; i < m; ++i) b[i] = g[m - i - 1];
    PolyInv(a, c, n - m + 2);
    int len = 1ll, lgm = 0ll; while(len < (n << 1)) {len <<= 1; ++lgm;}
    for(int i = 1; i < len; ++i) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (lgm - 1));
    NTT(a, len, 1); NTT(c, len, 1);
    for(int i = 0; i < len; ++i) q[i] = a[i] * c[i] % mod;
    NTT(q, len, -1); std::reverse(q, q + n - m + 1);
    for(int i = 0; i < len; ++i) a[i] = 0;
    for(int i = n - m + 1; i < len; ++i) q[i] = 0;
    for(int i = m; i < len; ++i) g[i] = 0;
    NTT(q, len, 1); NTT(g, len, 1);
    for(int i = 0; i < len; ++i) a[i] = q[i] * g[i] % mod;
    NTT(q, len, -1); NTT(g, len, -1); NTT(a, len, -1);
    for(int i = 0; i < m - 1; ++i) r[i] = (f[i] - a[i] + mod) % mod;
    for(int i = 0; i < len; ++i) a[i] = b[i] = c[i] = 0;
}

更新日志&引用

• 2019/9/27, 更新至多项式求逆
• 2019/9/28, 更新至多项式开根
• 2020/7/25, 更新至多项式除法

作者是完完全全跟着Venus和M-sea大爷学的多项式,所以本篇文章有很多引用来自她们的博客.

再次感谢这两篇博客给予本菜鸡在多项式学习方面的帮助!

摘要来自NaCly_Fish的Luogu个人空间.

To be continue