[BZOJ5020][THUWC 2017]在美妙的数学王国中畅游（LCT + 一点数学知识）

Address

• 洛谷 P4546
• BZOJ 5020
• LOJ #2289

Solution

• 如果只有一次函数 $ax+b$ ，那么这是非常裸的 LCT ，维护 $a$ 之和与 $b$ 之和即可
• 然后你会发现
• $$e^{a_{1}x+b_1}\times e^{a_{2}x+b_2}=e^{(a_{1}x+b_1)(a_{2}x+b_2)}=e^{a_1{a}_2{x}^2+(a_1{b}_2+a_2{b}_1)x+b_1{b}_2}$$
• 幂中有 $x^2$ 的项，照这样的话如果一条路径很长，那么使用到的 $x$ 的最高次幂与路径上的点数相同，无法简单维护
• $\sin$ 甚至没有可加性
• 回到 $ax+b$ ，我们可以很容易地推广到，如果不是一次函数而是 $m$ 次多项式（ $m$ 较小），那么我们同样可以通过维护多项式每一次项系数的和，做到 $O(nm\log n)$ 的优秀复杂度
• 这让我们思考能不能把 $\sin$ 和 $\exp$ 转化成多项式
• 先科普下泰勒展开
• $$f(x)=\sum _{i=0}^{\infty }\frac{f^{(i)}(x_0)(x-x_0{)}^i}{i!}$$
• 其中 $f^{(i)}(x)$ 表示 $f(x)$ 的 $i$ 阶导数
• 一般地，为了让函数转成多项式的形式，我们通常取 $x_0=0$
• 而 $\sin (ax+b)$ 和 $e^{ax+b}$ 是复合函数
• 可以按照 $h(x)=f(g(x)),h^{\prime }(x)=f^{\prime }(g(x))g^{\prime }(x)$ 进行求导
• 于是我们得到
• $$\sin (ax+b)=\sum _{i=0}^{\infty }\frac{a^i\text{orz}(i,b)}{i!}{x}^i$$
• 其中
• $$\text{orz}(i,x)=\{\begin{array}{ll}\sin x& i\mathrm{m}\mathrm{o}\mathrm{d}4=0\\ \cos x& i\mathrm{m}\mathrm{o}\mathrm{d}4=1\\ -\sin x& i\mathrm{m}\mathrm{o}\mathrm{d}4=0\\ -\cos x& i\mathrm{m}\mathrm{o}\mathrm{d}4=3\end{array}$$
• $$e^{ax+b}=\sum _{i=0}^{\infty }\frac{a^i{e}^b}{i!}{x}^i$$
• 注意到当 $i>16$ 时，上面两个多项式的 $i$ 次项系数已经被 $\frac{1}{i!}$ 压得非常小，其误差可以忽略
• 所以我们只需要维护出泰勒展开后的多项式的次数 $\le 16$ 的项的系数和即可
• 如果维护多项式的前 $m$ 项，则复杂度为 $O(nm\log n)$
• 将 $m$ 减小能优化复杂度，但误差会变大
• $m$ 取 $16$ 左右较合适

Code

// luogu-judger-enable-o2
#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
#define For(i, a, b) for (i = a; i <= b; i++)
#define Rof(i, a, b) for (i = a; i >= b; i--)
#define LCT(y, x) for (y = 0; x; y = x, x = fa[x])

template <class T>
inline void Swap(T &a, T &b) {T t = a; a = b; b = t;}

const int N = 1e5 + 5, M = 18;

int n, m, lc[N], rc[N], fa[N], rev[N], len, que[N];
double val[N][M], sum[N][M];
char op[M];

void getsin(double a, double b, double *val)
{
	int i; double tmp = 1, _sin = sin(b), _cos = cos(b);
	For (i, 0, 16)
	{
		if (i) tmp /= i;
		double st = i & 1 ? _cos : _sin;
		if ((i >> 1) & 1) val[i] = -st * tmp;
		else val[i] = st * tmp;
		tmp *= a;
	}
}

void getexp(double a, double b, double *val)
{
	int i; double tmp = 1, _exp = exp(b);
	For (i, 0, 16)
	{
		if (i) tmp /= i;
		val[i] = _exp * tmp;
		tmp *= a;
	}
}

void getlin(double a, double b, double *val)
{
	int i;
	val[0] = b; val[1] = a;
	For (i, 2, 16) val[i] = 0;
}

int which(int x) {return rc[fa[x]] == x;}

bool isroot(int x)
{
	return !fa[x] || (lc[fa[x]] != x && rc[fa[x]] != x);
}

void down(int x)
{
	if (rev[x])
	{
		rev[x] = 0;
		Swap(lc[x], rc[x]);
		if (lc[x]) rev[lc[x]] ^= 1;
		if (rc[x]) rev[rc[x]] ^= 1;
	}
}

void upt(int x)
{
	int i;
	For (i, 0, 16)
	{
		sum[x][i] = val[x][i];
		if (lc[x]) sum[x][i] += sum[lc[x]][i];
		if (rc[x]) sum[x][i] += sum[rc[x]][i];
	}
}

void rotate(int x)
{
	int y = fa[x], z = fa[y], b = lc[y] == x ? rc[x] : lc[x];
	if (!isroot(y)) (lc[z] == y ? lc[z] : rc[z]) = x;
	fa[x] = z; fa[y] = x; if (b) fa[b] = y;
	if (lc[y] == x) rc[x] = y, lc[y] = b;
	else lc[x] = y, rc[y] = b; upt(y); upt(x);
}

void splay(int x)
{
	int i, y = x;
	que[len = 1] = x;
	while (!isroot(y)) que[++len] = fa[y], y = fa[y];
	Rof (i, len, 1) down(que[i]);
	while (!isroot(x))
	{
		if (!isroot(fa[x]))
		{
			if (which(x) == which(fa[x])) rotate(fa[x]);
			else rotate(x);
		}
		rotate(x);
	}
}

void access(int x)
{
	int y;
	LCT(y, x)
	{
		splay(x); rc[x] = y;
		if (y) fa[y] = x; upt(x);
	}
}

int findroot(int x)
{
	access(x); splay(x);
	while (down(x), lc[x]) x = lc[x];
	return splay(x), x;
}

void makeroot(int x)
{
	access(x); splay(x);
	rev[x] ^= 1;
}

void link(int x, int y)
{
	makeroot(x); fa[x] = y;
}

void cut(int x, int y)
{
	makeroot(x); access(y); splay(y);
	lc[y] = fa[x] = 0; upt(y);
}

double query(int x, int y, double val)
{
	int i; double res = 0, tmp = 1;
	makeroot(x); access(y); splay(y);
	For (i, 0, 16) res += tmp * sum[y][i], tmp *= val;
	return res;
}

int main()
{
	int i, u, v, x; double s, t;
	scanf("%d%d%*s", &n, &m);
	For (i, 1, n)
	{
		scanf("%d%lf%lf", &x, &s, &t);
		if (x == 1) getsin(s, t, val[i]);
		else if (x == 2) getexp(s, t, val[i]);
		else getlin(s, t, val[i]);
		upt(i);
	}
	while (m--)
	{
		scanf("%s%d%d", op + 1, &u, &v);
		if (op[1] == 'a') link(u + 1, v + 1);
		else if (op[1] == 'd') cut(u + 1, v + 1);
		else if (op[1] == 'm')
		{
			scanf("%lf%lf", &s, &t);
			u++; splay(u);
			if (v == 1) getsin(s, t, val[u]);
			else if (v == 2) getexp(s, t, val[u]);
			else getlin(s, t, val[u]);
			upt(u);
		}
		else
		{
			scanf("%lf", &s); u++; v++;
			if (findroot(u) != findroot(v))
				puts("unreachable");
			else printf("%.8lf\n", query(u, v, s));
		}
	}
	return 0;
}