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; }