zoukankan      html  css  js  c++  java
  • THUWC2017 在美妙的数学王国中畅游

    题目链接


    这个提示貌似没有什么用.
    学过泰勒展开就会了.没学过可能不太行.
    第一第二个操作就是(LCT)板子.
    第三个就在(LCT)上改改就好了.
    第四个好像很麻烦的样子.
    首先,我们要知道它要求的是什么.
    这个(e^x),(sin(x))是什么玩意儿啊
    好像不资瓷合并啊
    但是我们发现,它不要求很高的精度.
    而且下面给了一个提示.
    因此学过泰勒展开的同学就知道该怎么做了.
    具体而言,我们将其转化为一个多项式.
    由于后面的项中,(x)的次数很高,系数很小,因此对答案没有什么影响.
    如果我们要求(10^{-7})的精度的话,大约只要记录(13)~(14)项.


    如何转换成多项式呢?
    先泰勒展开,然后再暴力二项式定理就好了.
    其实这个公式我们学校在上(FFT)的时候讲过

    [e^x=sum_{i=0}^{infty}frac{x^i}{i!}\ e^{ax+b}=sum_{i=0}^{infty}frac{(ax+b)^i}{i!}\ =sum_{i=0}^{infty}frac{sum_{j=0}^iC_i^ja^jb^{i-j}x^j}{i!}\ ]

    这两个公式应该还挺有名的.

    [sin(x)=sum_{i=0}^{infty}(-1)^ifrac{x^{2i+1}}{(2i+1)!}\ sin(ax+b)=sum_{i=0}^{infty}(-1)^ifrac{(ax+b)^{2i+1}}{(2i+1)!}\ =sum_{i=0}^{infty}(-1)^ifrac{(ax+b)^{2i+1}}{(2i+1)!}\ =sum_{i=0}^{infty}(-1)^ifrac{sum_{j=0}^{2i+1}C_{2i+1}^ja^jb^{2i+1-j}x^j}{(2i+1)!} ]

    这个和上面应该差不多.
    不过求(sin)的时候注意要多算一点,避免炸精度.


    代码如下
    但是貌似不够清真
    而且跑得超级慢

    #include<iostream>
    #include<cstdio>
    #include<algorithm>
    #include<cstring>
    #include<cmath>
    #include<vector>
    #include<queue>
    #include<map>
    #include<bitset>
    #include<set>
    #define N (200010)
    #define P ()
    #define M ()
    #define inf (0x7f7f7f7f)
    #define rg register int
    #define Label puts("NAIVE")
    #define spa print(' ')
    #define ent print('
    ')
    #define rand() (((rand())<<(13))^(rand()))
    typedef long double ld;
    typedef long long LL;
    typedef unsigned long long ull;
    using namespace std;
    namespace fastIO2{
    	template<class T>
    	inline void read(T &x){
    		static bool iosig;
    		static char c;
    		for(iosig=false,c=getchar();!isdigit(c);c=getchar()){
    			if(c=='-')iosig=true;
    			if(c==-1)return;
    		}
    		for(x=0;isdigit(c);c=getchar())x=((x+(x<<2))<<1)+(c^'0');
    		if(iosig)x=-x;
    	}
    }
    using namespace fastIO2;
    ld jc[31],C[31][31];
    struct poly{
    	ld a[14];
    	int n=13,op;
    	void change(int tp,ld x,ld y){
    		op=tp;
    		ld mia[28],mib[28];
    		if(tp==3){
    			for(int i=2;i<=n;i++)
    			a[i]=0; a[1]=x,a[0]=y;
    		}
    		else if(tp==2){
    			for(int i=0;i<=n;i++)a[i]=0;
    			mia[0]=1,mib[0]=1;
    			for(int i=1;i<=n;i++)
    			mia[i]=mia[i-1]*x,mib[i]=mib[i-1]*y;
    			for(int i=0;i<=n;i++)
    			for(int j=0;j<=i;j++)
    			a[j]=a[j]+mib[i-j]*mia[j]/jc[i]*C[i][j];
    		}
    		else if(tp==1){
    			for(int i=0;i<=n;i++)a[i]=0;
    			mia[0]=1,mib[0]=1;
    			for(int i=1;i<=n*2+1;i++)
    			mia[i]=mia[i-1]*x,mib[i]=mib[i-1]*y;
    			for(int i=0;i<=13;i++)
    			for(int j=0;j<=2*i+1;j++){
    				if(j>13)break;
    				ld k=mia[j]*mib[2*i+1-j]/jc[2*i+1]*C[2*i+1][j];
    				if(i&1)k=-k; a[j]+=k;
    			}
    		}
    	}
    	ld query(ld val){
    		ld t=1,ans=0;
    		for(int i=0;i<=n;i++)
    		ans+=t*a[i],t*=val;
    		return ans;
    	}
    	friend poly operator +(poly a,poly b){
    		poly res;
    		for(int i=0;i<=13;i++)
    		res.a[i]=a.a[i]+b.a[i];
    		return res;
    	}
    };
    struct node{
    	bool rev,rt;
    	int s[2],fa;
    	poly sum,val;
    }T[N];
    bool son(int fa,int x){return T[fa].s[1]==x;}
    void rev(int x){if(x)swap(T[x].s[0],T[x].s[1]),T[x].rev^=1;}
    void pushdown(int x){if(T[x].rev)rev(T[x].s[0]),rev(T[x].s[1]),T[x].rev=0;}
    void pushup(int x){T[x].sum=T[T[x].s[0]].sum+T[T[x].s[1]].sum+T[x].val;}
    void rotate(int x){
        if(T[x].rt)return;
        int a=T[x].fa,b=T[a].fa,p=son(a,x),q=p^1;
        pushdown(a);pushdown(x),T[a].s[p]=T[x].s[q];
        if(T[x].s[q])T[T[x].s[q]].fa=a;T[x].s[q]=a,T[a].fa=x,T[x].fa=b;
        if(!T[a].rt)T[b].s[son(b,a)]=x;else T[x].rt=1,T[a].rt=0;
        pushup(a),pushup(x);
    }
    void push(int x){if(!T[x].rt)push(T[x].fa);pushdown(x);}
    void Splay(int x){
        push(x);
        for(int fa;!T[x].rt;rotate(x))
        if(!T[fa=T[x].fa].rt)rotate((son(T[x].fa,x)==son(T[fa].fa,fa))?fa:x);
        pushup(x);
    }
    void access(int x){
        int y=0;
        while(x){
            Splay(x),T[T[x].s[1]].rt=1;
            T[T[x].s[1]=y].rt=false,pushup(x);
            y=x,x=T[x].fa;
        }
    }
    void mroot(int x){access(x),Splay(x),rev(x);}
    void Link(int x,int y){mroot(x),T[x].fa=y,Splay(x);}
    int find(int x){access(x),Splay(x);for(;T[x].s[0];x=T[x].s[0]);return x;}
    void Cut(int x,int y){
        mroot(x),access(y),Splay(y),pushdown(y);
        T[x].fa=T[y].s[0]=0,T[x].rt=1,pushup(x),pushup(y);
    }
    void modify(int x,int tp,ld a,ld b){
        access(x),Splay(x);
        T[x].val.change(tp,a,b),pushup(x);
    }
    bool check(int x,int y){
        if(find(x)==find(y))return 1;
        else return puts("unreachable"),0;
    }
    void query(int x,int y,ld v){
        if(!check(x,y))return;
        mroot(x),access(y),Splay(y);
        printf("%.10Lf
    ",T[y].sum.query(v));
    }
    int n,zz,m;
    int main(){
    	read(n),read(m),read(zz),C[0][0]=1;
    	jc[0]=1; for(int i=1;i<=27;i++)jc[i]=jc[i-1]*i;
    	for(int i=1;i<=27;i++){
    		C[i][0]=1;
    		for(int j=1;j<=i;j++)
    		C[i][j]=C[i-1][j]+C[i-1][j-1];
    	}
    	for(int i=1;i<=n;i++)T[i].rt=1;
    	for(int i=1;i<=n;i++){
    		ld a,b; int tp;
    		read(tp),scanf("%Lf%Lf",&a,&b);
    		modify(i,tp,a,b);
    	}
    	while(m--){
    		char ch=getchar();
    		while(ch<'a'||ch>'z')ch=getchar();
    		if(ch=='a'){
    			int x,y; read(x),read(y),Link(++x,++y);
    		}
    		else if(ch=='d'){
    			int x,y; read(x),read(y),Cut(++x,++y);
    		}
    		else if(ch=='m'){
    			int tp,x; ld a,b;
    			read(x),read(tp),scanf("%Lf%Lf",&a,&b);
    			modify(++x,tp,a,b);
    		}
    		else if(ch=='t'){
    			int x,y; ld v;
    			read(x),read(y),scanf("%Lf",&v);
    			query(++x,++y,v);
    		}
    	}
    }
    
  • 相关阅读:
    Web前端开发的一点记录
    关于表情符号与UTF-8的探讨
    微信公众号H5游戏本地搭建环境和调试(花生壳+wampserver+微信公众号测试号)
    ES6下的数组去重
    关于JS数组循环删除元素出现下标不对的问题
    好资源收藏
    ftp 操作,支持断点续传或者继续下载。
    推荐 greenrobot eventbus,简化安卓开发,提高安卓维护性,优化安卓性能
    T4模板TextTemplatingFileGenerator和TextTemplatingFilePreprocessor区别
    .net使用FluentValidation进行服务端验证。
  • 原文地址:https://www.cnblogs.com/Romeolong/p/10278437.html
Copyright © 2011-2022 走看看