【BZOJ5020】【THUWC2017】在美妙的数学王国中畅游 LCT 泰勒展开

题目大意

　　给你一棵树，每个点有一个函数$f(x)$

• 正弦函数 $\sin(ax+b) (a\in[0,1],b\in[0,\pi],a+b\in[0,\pi])$
• 指数函数 $e^{ax+b} (a\in[−1,1],b\in[−2,0],a+b\in[−2,0])$

• 一次函数 $ax+b (a\in[−1,1],b\in[0,1],a+b\in[0,1])$

  　　还有一些操作：

• 操作1：连接两个点（保证连接完后还是森林）

• 操作2：断开两个点之间的边
• 操作3：修改某一个点的函数

• 操作4：询问两个点路径上的所有函数的$f(v)$的和。

  　　$n\leq 100000,m\leq 200000,0\leq v\leq 1,0\leq f(v)\leq 1$

题目描述

　　前面两个操作就是link和cut。

　　我们需要在每个点上维护这个函数，但这些都不是多项式函数。我们可以发现$0\leq v\leq 1$，所以可以暴力维护这些函数在$x=0$处的泰勒展开式（就是生成函数）。我维护了前面$15$项。

$$\begin{align} e^x&=\sum_{i=0}^\infty \frac{x^i}{i!}\\ \sin(x)&=\sum_{i=0}^\infty{(-1)}^i\frac{x^{2i+1}}{(2i+1)!} \end{align}$$
　　每个函数的 $x=av+b$，直接用二项式定理暴力展开就可以了。

　　时间复杂度：$O(len\times m\log n)$

代码

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cstdlib>
#include<ctime>
#include<utility>
#include<cmath>
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
typedef pair<int,int> pii;
typedef pair<ll,ll> pll;
void sort(int &a,int &b)
{
    if(a>b)
        swap(a,b);
}
void open(const char *s)
{
#ifdef DEBUG
    char str[100];
    sprintf(str,"%s.in",s);
    freopen(str,"r",stdin);
    sprintf(str,"%s.out",s);
    freopen(str,"w",stdout);
#endif
}
const int len=15;
namespace lct
{
    double v[100010][16];
    double s[100010][16];
    int a[100010][2];
    int f[100010];
    int r[100010];
    void cp(int x)
    {
        memcpy(s[x],v[x],sizeof(double)*(len+1));
    }
    int root(int x)
    {
        return !f[x]||(a[f[x]][0]!=x&&a[f[x]][1]!=x);
    }
    void reverse(int x)
    {
        swap(a[x][0],a[x][1]);
        r[x]^=1;
    }
    void push(int x)
    {
        if(r[x])
        {
            if(a[x][0])
                reverse(a[x][0]);
            if(a[x][1])
                reverse(a[x][1]);
            r[x]=0;
        }
    }
    void mt(int x)
    {
        int i;
        int ls=a[x][0];
        int rs=a[x][1];
        for(i=0;i<=len;i++)
            s[x][i]=v[x][i]+s[ls][i]+s[rs][i];
    }
    void rotate(int x)
    {
        int p=f[x];
        int q=f[p];
        int ps=(x==a[p][1]);
        int qs=(p==a[q][1]);
        int ch=a[x][ps^1];
        if(!root(p))
            a[q][qs]=x;
        a[x][ps^1]=p;
        a[p][ps]=ch;
        if(ch)
            f[ch]=p;
        f[p]=x;
        f[x]=q;
        mt(p);
        mt(x);
    }
    void pushdown(int x)
    {
        if(!root(x))
            pushdown(f[x]);
        push(x);
    }
    void splay(int x)
    {
        pushdown(x);
        while(!root(x))
        {
            int p=f[x];
            if(!root(p))
            {
                int q=f[p];
                if((x==a[p][1])==(p==a[q][1]))
                    rotate(p);
                else
                    rotate(x);
            }
            rotate(x);
        }
    }
    void access(int x)
    {
        int y=x,t=0;
        while(x)
        {
            splay(x);
            a[x][1]=t;
            mt(x);
            t=x;
            x=f[x];
        }
        splay(y);
    }
    void change(int x)
    {
        access(x);
        reverse(x);
    }
    void link(int x,int y)
    {
        change(x);
        f[x]=y;
        splay(x);
    }
    void cut(int x,int y)
    {
        change(x);
        access(y);
        f[a[y][0]]=0;
        a[y][0]=0;
        mt(y);
    }
    int find(int x)
    {
        access(x);
        while(a[x][0])
            x=a[x][0];
        splay(x);
        return x;
    }
    double query(int x,int y,double v)
    {
        double res=0;
        change(x);
        access(y);
        int i;
        for(i=len;i>=0;i--)
        {
            res*=v;
            res+=s[y][i];
        }
        return res;
    }
}
double c[16][16];
double fac[16];
void init()
{
    int i,j;
    fac[0]=1;
    for(i=1;i<=len;i++)
        fac[i]=fac[i-1]*i;
    for(i=0;i<=len;i++)
    {
        c[i][0]=1;
        for(j=1;j<=i;j++)
            c[i][j]=c[i-1][j-1]+c[i-1][j];
    }
}
inline double get(int t,int x)
{
    if(t==1)
    {
        if(x%4==1)
            return 1/fac[x];
        else if(x%4==3)
            return -1/fac[x];
        return 0;
    }
    else if(t==2)
    {
        return 1/fac[x];
    }
    else
    {
        if(x==1)
            return 1;
        return 0;
    }
}
void change(double *s,int t,double a,double b)
{
    static double pa[16],pb[16];
    int i,j;
    for(i=0;i<=len;i++)
        s[i]=0;
    pa[0]=pb[0]=1;
    for(i=1;i<=len;i++)
    {
        pa[i]=pa[i-1]*a;
        pb[i]=pb[i-1]*b;
    }
    double v;
    if(t==1)
    {
        for(i=1;i<=len;i+=2)
        {
            v=get(t,i);
            for(j=0;j<=i;j++)
                //(ax+b)^t=C(t,i)a^ib^(t-i)
                s[j]+=v*c[i][j]*pa[j]*pb[i-j];
        }
    }
    else if(t==2)
    {
        for(i=0;i<=len;i++)
        {
            v=get(t,i);
            for(j=0;j<=i;j++)
                //(ax+b)^t=C(t,i)a^ib^(t-i)
                s[j]+=v*c[i][j]*pa[j]*pb[i-j];
        }
    }
    else
    {
        s[0]=b;
        s[1]=a;
    }
}
int main()
{
    open("bzoj5020");
    init();
    int n,m;
    char type[100];
    scanf("%d%d%s",&n,&m,type);
    double a,b,v;
    int x,y,t;
    int i;
    for(i=1;i<=n;i++)
    {
        scanf("%d%lf%lf",&t,&a,&b);
        change(lct::v[i],t,a,b);
        lct::cp(i);
    }
    for(i=1;i<=m;i++)
    {
        scanf("%s",type);
        if(type[0]=='a')
        {
            scanf("%d%d",&x,&y);
            x++;
            y++;
            lct::link(x,y);
        }
        else if(type[0]=='d')
        {
            scanf("%d%d",&x,&y);
            x++;
            y++;
            lct::cut(x,y);
        }
        else if(type[0]=='m')
        {
            scanf("%d%d%lf%lf",&x,&t,&a,&b);
            x++;
            lct::access(x);
            change(lct::v[x],t,a,b);
            lct::cp(x);
        }
        else
        {
            scanf("%d%d%lf",&x,&y,&v);
            x++;
            y++;
            if(lct::find(x)!=lct::find(y))
            {
                printf("unreachable\n");
                continue;
            }
            double ans=lct::query(x,y,v);
            printf("%.8le\n",ans);
        }
    }
    return 0;
}