BZOJ5020: [THUWC 2017]在美妙的数学王国中畅游
其实题面好像有点不全，建议去洛谷：
P4546 [THUWC2017]在美妙的数学王国中畅游
这里还是$BZOJ$的题面。
Description
Input
Output
Sample Input
1 1 0
3 0.5 0.5
3 -0.5 0.7
appear 0 1
travel 0 1 0.3
appear 0 2
travel 1 2 0.5
disappear 0 1
appear 1 2
travel 1 2 0.5
Sample Output
1.67942554e+000
1.20000000e+000
题解Here!
数字和数学规律主宰着这个世界。。。
学好数理化，走遍天下都不怕。。。
其实原题面有一个提示，我们称它为——泰勒展开。
如下：（选自洛谷）
若函数$f(x)$的$n$阶导数在$[a,b]$区间内连续，则对$f(x)$在$x_0(x_0in[a,b])$处使用$n$次拉格朗日中值定理可以得到带拉格朗日余项的泰勒展开式：
$$f(x)=f(x_0)+frac{f'(x_0)(x-x_0)}{1!}+frac{f''(x_0)(x-x_0)^2}{2!}+ cdots +frac{f^{(n-1)}(x_0)(x-x_0)^{n-1}}{(n-1)!}+frac{f^{(n)}(xi)(x-x_0)^n}{n!},xin[a,b]$$
其中，当$x>x_0$时，$xiin[x_0,x]$；当$x<x_0$时，$xiin[x,x_0]$。
$f^{(n)}(x)$表示函数$f(x)$的$n$阶导数。
其实别看这么多乱七八糟的名词，关键就在于那个式子。
我们注意到，等号的左边是$f(x)$！
也就是说，我们可以不用搞其它的东西，只要把$f(x)$的$n$阶导数搞出来就好。
什么？你不知道导数？
回去学《数学 选修2-2》去吧。。。
这里给出此题要用到的导数计算：
基本公式：
$$f(x)=xquadRightarrowquad f'(x)=1$$
$$f(x)=sin(x)quadRightarrowquad f'(x)=cos(x)$$
$$f(x)=cos(x)quadRightarrowquad f'(x)=-sin(x)$$
$$f(x)=e^xquadRightarrowquad f'(x)=e^x$$
计算公式：
$$[f(x)+g(x)]'=f'(x)+g'(x)$$
$$[f(x)g(x)]'=f'(x)g(x)+f(x)g'(x)$$
$$ ext{链式法则：}[f(g(x))]'=f'(g(x)) imes g'(x)$$
于是我们可以愉快地推式子了！
以$f(x)=sin(x)$为例：
泰勒展开式中选取了一个$x_0$，我们为了方便，直接$x_0=0$。
设$f^n(x)$为$f(x)$的$n$阶导数。
那么：$$f(x)=f(0)+frac{f^1(0)(x-0)}{1!}+frac{f^2(0)(x-0)^2}{2!}+ cdots +frac{f^{n-1}(0)(x-0)^{n-1}}{(n-1)!}+frac{f^{n}(xi)(x-0)^n}{n!}$$
化简一下：$$f(x)=f(0)+frac{f^1(0)}{1!}x+frac{f^2(0)}{2!}x^2+ cdots +frac{f^{n-1}(0)}{(n-1)!}x^{n-1}+frac{f^{n}(xi)}{n!}x^n$$
但是这个$n$到底是多大呢？
其实，我们发现$n!$在$n$比较小时就是暴增函数。
也就是说，越往后的项对答案贡献越小，而精度也有$SPJ$，所以我们并不需要末尾那么多项。
这里我选取了$n=16$这个值，这是可以过的。
于是我们只要对每个点维护$16$项多项式就好。
但是又有人问了：题目中不是求$f(x)=sin(ax+b)$吗？
其实根据上面那个链式法则我们可以得出：$$f(x)=sin(ax+b)quadRightarrowquad f'(x)=acos(ax+b)$$
同理对于另外两个式子有：
$$f(x)=ax+bquadRightarrowquad f'(x)=a$$
$$f(x)=e^{ax+b}quadRightarrowquad f'(x)=ae^{ax+b}$$
但是这只是一阶导数啊，高阶导数怎么办？
对于$f(x)=ax+b$，只有一阶导数$f'(x)=a$，高阶导数全部为$0$。
对于$f(x)=e^{ax+b}$，一阶导数$f'(x)=ae^{ax+b}$，然后每高一阶多乘一个$a$，如二阶导数为$f''(x)=a^2e^{ax+b}$。
对于$f(x)=sin(ax+b)$，一阶导数$f'(x)=acos(ax+b)$，然后每高一阶多乘一个$a$，$ ext{阶数}mod 2==0$的导数为$sin$，否则为$cos$，并且$ ext{阶数}mod 4leq1$的导数为正，否则为负。
由于精度会出锅，再加上每次计算$sin(x),cos(x),e^x$会比较慢，于是每次计算之前预处理一下即可。
至于阶乘。。。我一开始写了个逆元，发现模数不知道。。。
我这个大制杖好久才反应过来$16!=2.0922789888 imes 10^{13}<10^{18}$。。。
所以直接开$double$预处理阶乘就好。
$LCT$应该不用多说，维护多项式的每一项单独的和。
答案就是多项式每一项之和。
至于那个什么科学计数法。。。大多数人直接$\%.8lf$。。。
其实我们可以直接用$C++$的自带版本——$\%.8e$！
自动将小数转成科学计数法！
还有啥细节那就看代码吧。。。
附代码：
#include<iostream>
#include<algorithm>
#include<cstdio>
#include<cmath>
#include<cstring>
#define MAXN 100010
#define MAXM 20
using namespace std;
int n,m;
double fact[MAXM],val[MAXN][MAXM];
int top,stack[MAXN];
struct Link_Cut_Tree{
int f,flag,son[2];
double v[MAXM];
}a[MAXN];
inline int read(){
int date=0,w=1;char c=0;
while(c<'0'||c>'9'){if(c=='-')w=-1;c=getchar();}
while(c>='0'&&c<='9'){date=date*10+c-'0';c=getchar();}
return date*w;
}
inline bool isroot(int rt){
return a[a[rt].f].son[0]!=rt&&a[a[rt].f].son[1]!=rt;
}
inline void pushup(int rt){
if(!rt)return;
for(int i=0;i<=15;i++)a[rt].v[i]=a[a[rt].son[0]].v[i]+a[a[rt].son[1]].v[i]+val[rt][i];
}
inline void pushdown(int rt){
if(!rt||!a[rt].flag)return;
a[a[rt].son[0]].flag^=1;a[a[rt].son[1]].flag^=1;a[rt].flag^=1;
swap(a[rt].son[0],a[rt].son[1]);
}
inline void turn(int rt){
int x=a[rt].f,y=a[x].f,k=a[x].son[0]==rt?1:0;
if(!isroot(x)){
if(a[y].son[0]==x)a[y].son[0]=rt;
else a[y].son[1]=rt;
}
a[rt].f=y;a[x].f=rt;a[a[rt].son[k]].f=x;
a[x].son[k^1]=a[rt].son[k];a[rt].son[k]=x;
pushup(x);pushup(rt);
}
void splay(int rt){
top=0;
stack[++top]=rt;
for(int i=rt;!isroot(i);i=a[i].f)stack[++top]=a[i].f;
while(top)pushdown(stack[top--]);
while(!isroot(rt)){
int x=a[rt].f,y=a[x].f;
if(!isroot(x)){
if((a[y].son[0]==x)^(a[x].son[0]==rt))turn(rt);
else turn(x);
}
turn(rt);
}
}
void access(int rt){
for(int i=0;rt;i=rt,rt=a[rt].f){
splay(rt);
a[rt].son[1]=i;
pushup(rt);
}
}
inline void makeroot(int rt){access(rt);splay(rt);a[rt].flag^=1;}
int findroot(int rt){
access(rt);splay(rt);
while(a[rt].son[0])rt=a[rt].son[0];
return rt;
}
inline void split(int x,int y){makeroot(x);access(y);splay(y);}
inline void link(int x,int y){makeroot(x);a[x].f=y;}
inline void cut(int x,int y){
split(x,y);
if(a[y].son[0]==x&&a[x].f==y&&!a[x].son[1])a[x].f=a[y].son[0]=0;
}
void change(int x,int f,double a,double b){
memset(val[x],0,sizeof(val[x]));
if(f==1){
double t=1,Sin=sin(b),Cos=cos(b);
for(int i=0;i<=15;i++){
val[x][i]=((i%2)?Cos:Sin)*((i%4<=1)?1:-1)*t/fact[i];
t*=a;
}
}
else if(f==2){
double t=1,Exp=exp(b);
for(int i=0;i<=15;i++){
val[x][i]=Exp*t/fact[i];
t*=a;
}
}
else{
val[x][0]=b;
val[x][1]=a;
}
}
double query(int x,int y,double k){
split(x,y);
double s=0,t=1;
for(int i=0;i<=15;i++){
s+=a[y].v[i]*t;
t*=k;
}
return s;
}
void work(){
int u,v;
double x,y,k;
char ch[20];
while(m--){
scanf("%s",ch);u=read();v=read();
switch(ch[0]){
case 'a':{
u++;v++;
link(u,v);
break;
}
case 'd':{
u++;v++;
cut(u,v);
break;
}
case 'm':{
u++;
access(u);splay(u);
scanf("%lf%lf",&x,&y);
change(u,v,x,y);
pushup(u);
break;
}
case 't':{
u++;v++;scanf("%lf",&k);
if(findroot(u)!=findroot(v))printf("unreachable
");
else printf("%.8e
",query(u,v,k));
break;
}
}
}
}
void init(){
int f;
double x,y;
char ch[5];
n=read();m=read();scanf("%s",ch);
fact[0]=1;
for(int i=1;i<=16;i++)fact[i]=fact[i-1]*i;
for(int i=1;i<=n;i++){
f=read();
scanf("%lf%lf",&x,&y);
change(i,f,x,y);
}
}
int main(){
init();
work();
return 0;
}