LOJ #150. 挑战多项式/多项式全家桶

Infinite_Jerry

于 2020-06-18 09:37:22 发布

阅读量465

点赞数

分类专栏： # 多项式

本文为博主原创文章，未经博主允许不得转载。一经发现，必将追究法律责任。

本文链接：https://blog.csdn.net/qq_42886072/article/details/106823704

版权

多项式专栏收录该内容

5 篇文章 0 订阅

订阅专栏

题目

这个题堪称简易模板全家桶。

关于各个函数的实现，大多数都是利用牛顿迭代公式( $x=x_0-\dfrac {f(x)}{f'(x)}$ )+倍增。

以下为多项式全家桶证明:

下面的 $f, b, a$ 分别为已知多项式，当前所求多项式和上一个状态(规模小一半)的多项式.(简记 $y = a (x)$ )

$\text{sqrt}$

$b(x)^2=f(x) ~~->b(x)^2-f(x)=0$
$设\varphi(b(x))=b(x)^2-f(x)$
$则\varphi'(b(x))=2*b(x)$
$此时, 我们求函数零点$
$b(x)=a(x)-\dfrac{\varphi(a(x))}{\varphi'(a(x))}=a(x)-\dfrac{a^2(x)-f(x)}{2a(x)}=\dfrac{a^2(x)+f(x)}{2a(x)}$

$\dfrac {y^2+f}{2y}$

$\text{inv}$

$\varphi(b(x))=b(x)f(x)-1$

$\varphi'(b(x))=f(x)$

$b(x)=a(x)-\dfrac{a(x)f(x)-1}{f(x)}$

$b (x) = a (x) - b (x) (a (x) f (x) - 1)$

由于 $a(x)f(x)-1\equiv 0(\mod x^{n/2})$ .

所以 $b (x)$ 的高 $n / 2$ 位乘上 $a (x) f (x) - 1$ 在 $mod x^n$ 意义下为0.

所以 $b (x)$ 与 $a (x)$ 在当前等价.

即 $b (x) = a (x) - a (x) (a (x) f (x) - 1) = a (x) (2 - a (x) f (x))$

$y (2 - y f)$

$\ln$

$\ln 'x=\dfrac 1 x$

$b(x)=\ln f(x)$

对两边求导: $b^{'} (x) = i n v (f (x)) f^{'} (x)$ .(链式反应)

积分 $b^{'} (x)$ 即可得到 $b (x)$ .

$\int \dfrac {f'}f$

$\exp$

因为 $\exp,\ln$ 为逆运算,所以可得:

$b(x)=\exp f(x)\rightarrow \ln b(x)=f(x)$

设 $\varphi(b(x))=\ln b(x)-f(x)$

$\varphi'(b(x))=\dfrac 1{b(x)}$

则可得到: $b(x)=a(x)(1-\ln a(x)+f(x))$ .

$y(1-\ln y+f)$

$\text{pow}$

$f(x)^k=e^{ln(f(x))*k}$

$\text{f(x)/g(x),f(x)mod g(x)}$

多项式带余除法.

$f(x)=q(x)g(x)+r(x)(\mod x^n)$ .

已知 $f, g$ 求 $q, r$ . $f$ 为 $n$ 次多项式, $g$ 为 $m$ 次多项式, $r$ 的次数 $< m$ .

可以发现 $x^n f(x^{-1})$ 的系数为 $f$ 系数的翻转,我们简记这个多项式为 $F (x)$ .

则有: $F(x)=x^n q(x^{-1})g(x^{-1})+x^nr(x^{-1})=x^{n-m} q(x^{-1}) x^m g(x^{-1})+x^{n-m+1}x^{m-1}r(x^{-1})= Q(x)G(x)+x^{n-m+1}R(x)(\mod x^n)$

我们可以发现 $F(x)=Q(x)G(x)(\mod x^{n-m+1})$ .

这样求出 $Q$ 后即可得到 $q, r$ .

代码:

upd on 2021/2/3,大改码风和修改一些实现可能的问题.

#include <bits/stdc++.h>
#define gc getchar()
#define TP template<class o>
#define TPP template<typename t1, typename t2>
#define rep(i, a, b) for (int i = a; i <= (int)(b); i++)
#define REP(i, a, b) for (int i = b; i >= (int)(a); i--)
#define FOR(i, n) rep(i, 1, n) 
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
const int N = 1 << 19 | 10, mod = 998244353;

TP void qr(o &x) {
	char c = gc; int f = 1;
	while(!isdigit(c)) { if(c == '-') f = -1; c = gc;}
	while( isdigit(c)) x = x * 10 + c - '0', c = gc;
	x *= f;
}
template<class o, class ... O> void qr(o& x, O& ... y) {qr(x); qr(y...);}
TP void qw(o x) {
	if(x < 0) putchar('-'), x = -x;
	if(x / 10) qw(x / 10);
	putchar(x % 10 + '0');
}
TP void pr1(o x) {qw(x); putchar(' ');}

TPP void ad(t1 &x, t2 y) { x += y; if(x >= mod) x -= mod;}
TPP void dl(t1 &x, t2 y) { x -= y; if(x < 0) x += mod;}

ll power(ll a, ll b = mod - 2, ll p = mod) {
	ll c = 1; while(b) {
		if(b & 1) 	c = c * a % p;
		b /= 2; 	a = a * a % p;
	}
	return c;
}

namespace Cipolla {
	ll p = mod, w;
	struct CP {
		ll x, y;
		CP(ll a = 1, ll b = 0) : x(a), y(b){};
		CP operator *(CP b) const { return CP((x * b.x + y * b.y % p * w) % p, (x * b.y + y * b.x) % p); }
	};
	CP power(CP a, ll b = (p + 1) / 2) {
		CP c; while(b) {
			if(b & 1) c = c * a;
			b /= 2; a = a * a;
		}
		return c;
	}
	ll solve(ll n) {
		if(n <= 1) return n;
		ll a; do {
			a = ((ll)rand() << 30 | rand()) % p;
			w = (a * a - n % p + p) % p;
		} while(::power(w, (p - 1) / 2, p) != p - 1);
		ll x = power(CP(a, 1)).x; return min(x, p - x);
	}
}

namespace poly {
	const int g = 3; const ll inv2 = (mod + 1) / 2;
	int R[N], w[N], Inv[N];
	int calc(int x) {
		if((x & -x) == x) return x;
		int n = 1; while(n < x) n *= 2;
		return n;
	}
	int pre(int m) {
		int n = calc(m);
		FOR(i, n - 1) R[i] = (R[i >> 1] >> 1) | (i & 1 ? n >> 1 : 0);
		return n;
	}
	void init(int m) {
		int n = calc(m) * 2;
		Inv[0] = Inv[1] = 1; rep(i, 2, n) Inv[i] = (ll)Inv[mod % i] * (mod - mod / i) % mod;
		for(int i = 1; i < n; i *= 2) {
			ll s = 1, d = power(g, (mod - 1) / (2 * i));
			rep(j, 0, i - 1) w[i + j] = s, s = s * d % mod;
		}
	}
	void DFT(int *f, int n) {
		static ull p[N]; rep(i, 0, n - 1) p[R[i]] = f[i];
		for(int i = 1, t; i < n; i *= 2) for(int j = 0; j < n; j += 2 * i)
			rep(k, 0, i - 1) t = p[j + k + i] * w[i + k] % mod, p[j + k + i] = p[j + k] + mod - t, p[j + k] += t;
		rep(i, 0, n - 1) f[i] = p[i] % mod;
	}
	void IDFT(int *f, int n) {
		reverse(f + 1, f + n); DFT(f, n); ll inv = Inv[n];
		rep(i, 0, n - 1) f[i] = f[i] * inv % mod;
	}
	const int sz = sizeof(int);
	void copy(int *a, int *b, int n) {memcpy(a, b, sz * n);}
	void clear(int *a, int n) { memset(a + n, 0, sz * n);}
	void clear(int *a, int x, int y) {if(x < y) memset(a + x, 0, sz * (y - x)); }
	
	void dao(int *a, int *b, int n) {
		FOR(i, n - 1) a[i - 1] = (ll)b[i] * i % mod;
		a[n - 1] = 0;
	}
	void ji(int *a, int *b, int n) {
		REP(i, 1, n - 1) a[i] = (ll)b[i - 1] * Inv[i] % mod;
		a[0] = 0;
	}
	void mult(int *a, int *b, int n, int m) {
		static int c[N]; int len = pre(n + m);
		copy(c, b, m); clear(c, m, len); DFT(a, len); DFT(c, len);
		rep(i, 0, len - 1) c[i] = (ll)a[i] * c[i] % mod;
		IDFT(c, len); copy(a, c, n + m);
	}
	void getinv(int *a, int *b, int n) {
		static int c[N]; clear(a, 0, 2 * n); clear(b, n); clear(c, 0, 2 * n); a[0] = power(b[0]);
		for(int p = 2; p <= n; p *= 2) {
			int len = pre(p * 2); copy(c, b, p); clear(c, p); DFT(a, len); DFT(c, len);
			rep(i, 0, len - 1) a[i] = (ll)a[i] * (mod + 2 - (ll)a[i] * c[i] % mod) % mod;
			IDFT(a, len); clear(a, p);
		}
	}
	void getsqrt(int *a, int *b, int n) {
		static int c[N];
		clear(a, 0, 2 * n); clear(b, n);
		clear(c, 0, 2 * n);
		a[0] = Cipolla::solve(b[0]);
		for(int p = 2; p <= n; p *= 2) {
			getinv(c, a, p); mult(c, b, p, p);
			rep(i, 0, p - 1) a[i] = (ll)(a[i] + c[i]) * inv2 % mod;
		}
	}
	void getln(int *a, int *b, int n) {
		static int c[N]; dao(a, b, n); getinv(c, b, n);
		mult(c, a, n, n); ji(a, c, n);
	}
	void getexp(int *a, int *b, int n) {
		static int c[N]; clear(c, 0, 2 * n); a[0] = 1;
		for(int p = 2; p <= n; p *= 2) {
			getln(c, a, p); dl(c[0], 1);
			rep(i, 0, p - 1)  dl(c[i] = -c[i] + b[i], 0);
			mult(c, a, p, p); copy(a, c, p);
		}
	}
}

int n, m, t, f[N], g[N], h[N];

int main() {
	qr(n); qr(m); t = poly::calc(++n); poly::init(t);
	rep(i, 0, n - 1) qr(f[i]);
	poly::getsqrt(g, f, t);
	poly::getinv (h, g, t);
	poly::ji(g, h, t);
	poly::getexp(h, g, t);
	FOR(i, t - 1) dl(f[i], h[i]);
	f[0] = 1; poly::getln(g, f, t); g[0]++;
	poly::getln(f, g, t); rep(i, 0, t - 1) f[i] = (ll)f[i] * m % mod;
	poly::getexp(g, f, t); poly::dao(f, g, t); n--;
	rep(i, 0, n - 1) pr1(f[i]);
	return 0;
}