「学习笔记」多项式相关

序

学多项式也有好久了，可是我自己还没怎么认认真真推过柿子，导致啥都不会，然后被吊打。

看来再不回顾一下就不行了啊。

多项式乘法

写了一个好看一点的 NTT 板子，仅供参考。

inline int Add(int x, int y) { return (x + y) % Mod; }
inline int Sub(int x, int y) { return (x - y + Mod) % Mod; }
inline int Mul(int x, int y) { return 1ll * x * y % Mod; }
int fastpow(int x, int y)
{
	int ans = 1;
	for (; y; y >>= 1, x = 1ll * x * x % Mod)
		if (y & 1) ans = 1ll * ans * x % Mod;
	return ans;
}

int r[maxn], w[maxn];
void FFT(int *p, int N)
{
	for (int i = 0; i < N; i++) if (i < r[i]) std::swap(p[i], p[r[i]]);
	for (int i = 1, s = 2, t = N >> 1; i < N; i <<= 1, s <<= 1, t >>= 1)
		for (int j = 0; j < N; j += s) for (int k = 0, o = 0; k < i; ++k, o += t)
		{
			int x = p[j + k], y = 1ll * w[o] * p[i + j + k] % Mod;
			p[j + k] = (x + y) % Mod, p[i + j + k] = (x - y + Mod) % Mod;
		}
}

template<typename C>
void PolyMul(int *a, int *b, int N, int P, C cal)
{
	w[0] = 1, w[1] = fastpow(3, (Mod - 1) / N);
	for (int i = 0; i < N; i++) r[i] = (r[i >> 1] >> 1) | ((i & 1) << P);
	for (int i = 2; i < N; i++) w[i] = 1ll * w[i - 1] * w[1] % Mod;
	FFT(a, N), FFT(b, N); for (int i = 0; i < N; i++) b[i] = cal(a[i], b[i]);
	FFT(b, N); std::reverse(b + 1, b + N); int invn = fastpow(N, Mod - 2);
	for (int i = 0; i < N; i++) b[i] = 1ll * b[i] * invn % Mod;
}

泰勒展开

如果(f(x))在(x_0)处存在(n)阶导，那么有：

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

其中(xi)是余项，当(n)趋近于无穷大时，(xi)趋近于高阶无穷小。

比如说(e ^ x = 1 + frac x{1!} + frac {x^2}{2!} + cdots)

牛顿迭代

首先可以知道多项式的任何一个运算都可以表示成对于一个多项式(B(x))以及一个给定的函数(F(x))，求(F(B(x)) equiv 0 pmod {x ^ n})

设(B_n(x))表示当模数是(x ^ n)的合法解。那么当(n = 1)是我们很容易可以得到结果，考虑如何用(B_n(x))推到(B_{2n}(x))。

对(F(B_{2n}(x)))在(B_n(x))处泰勒展开，我们得到(F(B_{2n}(x)) = F(B_n(x)) + F'(B_n(x))(B_{2n}(x) - B_n(x)))。

那么我们化简一下就是：

[B_{2n}(x) = B_n(x) - frac {F(B_n(x))} {F'(B_n(x))} ]

这样我们就可以倍增求解。

多项式运算

接下来均假设我们要做运算的多项式是(A(x))。

多项式求逆

令(F(B_n(x)) = A(x) * B_n(x) - 1 equiv 0)

于是：

[egin{aligned} B_{2n} &= B_n(x) - frac {A(x) * B_n(x) - 1} {A(x)} \ &= B_n(x) - B_n(x)(A(x) * B_n(x) - 1) \ &= 2B_n(x) - A(x) * B_n^2(x) end{aligned} ]

注意第一步推到第二步是因为(B_n(x))是(A(x))的逆。

void Inv(int *a, int *b, int N)
{
	static int c[maxn]; if (N == 1) return (void) (*b = fastpow(*a, Mod - 2));
	Inv(a, b, (N + 1) >> 1); int L = 1, P = -1; while (L < (N << 1)) L <<= 1, ++P;
	std::copy(a, a + N, c), std::fill(c + N, c + L, 0);
	PolyMul(c, b, L, P, [] (int a, int b) { return Mul(Sub(2, Mul(a, b)), b); });
	std::fill(b + N, b + L, 0);
}

多项式开根

ln一下再exp一下

令(F(B_n(x)) = B_n^2(x) - A(x) equiv 0)

于是有：

[B_{2n} = B_n(x) - frac {B_n^2(x) - A(x)} {2B_n(x)} = frac 12left(B_n(x) + frac {A(x)} {B_n(x)} ight) ]

可以看出多项式开根中需要套用多项式求逆。

void Sqrt(int *a, int *b, int N)
{
	static int c[maxn], d[maxn]; if (N == 1) return (void) (*b = 1);
	Sqrt(a, b, (N + 1) >> 1); int L = 1, P = -1; while (L < (N << 1)) L <<= 1, ++P;
	std::copy(a, a + N, c), std::fill(c + N, c + L, 0);
	std::fill(d, d + L, 0), Inv(b, d, N), PolyMul(c, d, L, P, Mul);
	for (int i = 0; i < N; i++) b[i] = Mul(Add(b[i], d[i]), 499122177);
}

多项式(ln)

不需要牛顿迭代。

[egin{aligned} ln(A(x)) &= B(x) \ Rightarrowfrac {A'(x)}{A(x)} &= B'(x) end{aligned} ]

然后(A(x))就求导求逆，乘起来再积分一下就可以了。

void Ln(int *f, int *g, int N)
{
	static int A[maxn], B[maxn]; Inv(f, B, N), A[N - 1] = 0;
	for (int i = 1; i < N; i++) A[i - 1] = Mul(f[i], i);
	int L = 1, P = -1; while (L < (N << 1)) L <<= 1, ++P;
	PolyMul(A, B, L, P, Mul), g[0] = 0;
	for (int i = 1; i < N; i++) g[i] = Mul(B[i - 1], fastpow(i, Mod - 2));
	std::fill(A, A + L, 0), std::fill(B, B + L, 0);
}

多项式(exp)

令(F(B_n(x)) = ln(B_n(x)) - A(x) equiv 0)

推下式子可得：

[egin{aligned} B_{2n}(x) &= B_n(x) - frac {ln(B_n(x)) - A(x)} {frac 1 {B_n(x)}} \ &= B_n(x)(1 - ln(B_n(x)) + A(x)) end{aligned} ]

void Exp(int *a, int *b, int N)
{
	static int c[maxn]; if (N == 1) return (void) (*b = 1);
	Exp(a, b, (N + 1) >> 1), Ln(b, c, N);
	int L = 1, P = -1; while (L < (N << 1)) L <<= 1, ++P;
	for (int i = 0; i < N; i++) c[i] = Sub(a[i], c[i]); c[0] = Add(c[0], 1);
	PolyMul(c, b, L, P, Mul), std::fill(b + N, b + L, 0);
}

多项式除法

给定一个(n)次多项式(A(x))和一个(m)次多项式(B(x))，求解一个(n - m)次的多项式(Q(x))以及一个小于(n - m)次的多项式(R(x))，使得(A(x) = Q(x)B(x) + R(x))

定义运算(R)使得(A^R(x) = x ^ nA(frac 1x))，也就是说将(A(x))的系数翻转。

那么可以得到：

[egin{aligned} A(x) &= Q(x)B(x) + R(x) \ Aleft(frac 1x ight) &= Qleft(frac 1x ight)Bleft(frac 1x ight) + Rleft(frac 1x ight) \ x ^ nAleft(frac 1x ight) &= left(x ^ m Bleft(frac 1x ight) ight) * left( x ^ {n - m} Qleft(frac 1x ight) ight) + x ^ n Rleft(frac 1x ight) \ A^R(x) &= Q^R(x)B^R(x) + x^{n-m+1}R^R(x) \ A^R(x) &equiv Q^R(x)B^R(x) pmod{x ^ {n - m + 1}} \ Q^R(x) &equiv frac {A^R(x)} {B^R(x)} pmod{x ^ {n - m + 1}} end{aligned} ]

那么这样就可以用多项式求逆求出(Q)，再用(R(x) = A(x) - Q(x)B(x))即可求出(R)。