[BZOJ4623][树分治][矩阵][数学]Styx

%%%%tlzmybm

大概就是把向量积转化成行列式表示，然后就可以用反交换律计算两点间的叉积，然后树分治一下

#include <cstdio>
#include <iostream>
#include <cstring>
#include <string>
#include <vector>
#include <algorithm>
#define N 100010
#define P 1000000007

using namespace std;

typedef pair<int,int> paris;
typedef vector<int>::iterator itr;

int n,m,u,v,cnt,Root,isize,Max;
int a[N],G[N],vis[N],pos[N],dpt[N],fa[N];
struct edge{
    int t,nx;
}E[N<<1];
struct Vector{
    int a[3];
    Vector(int x=0,int y=0,int z=0){ a[0]=x; a[1]=y; a[2]=z; }
    int &operator [](int x){ return a[x]; }
    inline void reverse(){ a[0]=P-a[0]; a[1]=P-a[1]; a[2]=P-a[2]; }
}V[N],Ans[N];
struct Matrix{
    int a[3][3];
    int *operator [](int x){ return a[x]; }
    Matrix(){ memset(a,0,sizeof(a)); }
    Matrix(Vector A){
        a[0][0]=0; a[0][1]=P-A[2]; a[0][2]=A[1];
        a[1][0]=A[2]; a[1][1]=0; a[1][2]=P-A[0];
        a[2][0]=P-A[1]; a[2][1]=A[0]; a[2][2]=0;
    }
    friend Matrix operator *(Matrix A,Matrix B){
        Matrix C;
        C[0][0]=(1ll*A[0][0]*B[0][0]+1ll*A[0][1]*B[1][0]+1ll*A[0][2]*B[2][0])%P;
        C[0][1]=(1ll*A[0][0]*B[0][1]+1ll*A[0][1]*B[1][1]+1ll*A[0][2]*B[2][1])%P;
        C[0][2]=(1ll*A[0][0]*B[0][2]+1ll*A[0][1]*B[1][2]+1ll*A[0][2]*B[2][2])%P;
        C[1][0]=(1ll*A[1][0]*B[0][0]+1ll*A[1][1]*B[1][0]+1ll*A[1][2]*B[2][0])%P;
        C[1][1]=(1ll*A[1][0]*B[0][1]+1ll*A[1][1]*B[1][1]+1ll*A[1][2]*B[2][1])%P;
        C[1][2]=(1ll*A[1][0]*B[0][2]+1ll*A[1][1]*B[1][2]+1ll*A[1][2]*B[2][2])%P;
        C[2][0]=(1ll*A[2][0]*B[0][0]+1ll*A[2][1]*B[1][0]+1ll*A[2][2]*B[2][0])%P;
        C[2][1]=(1ll*A[2][0]*B[0][1]+1ll*A[2][1]*B[1][1]+1ll*A[2][2]*B[2][1])%P;
        C[2][2]=(1ll*A[2][0]*B[0][2]+1ll*A[2][1]*B[1][2]+1ll*A[2][2]*B[2][2])%P;
        return C;
    }
}A[N],B[N],C,Mat[N];
struct PrimeTable{
    #define Maxx 1000010
    int g;
    int p[Maxx+10],c[Maxx+10],Inv[Maxx+10];
    PrimeTable(){ g=1; }
    inline void Init(){
        for(int i=2;i<=Maxx;i++){
            if(!p[i]) p[++*p]=i;
            for(int j=1;j<=*p&&1ll*p[j]*i<=Maxx;j++)
                if(p[p[j]*i]=1,i%p[j]==0) break;
        }
        Inv[0]=Inv[1]=1;
        for(int i=2;i<=Maxx;i++) 
            Inv[i]=1ll*(P-P/i)*Inv[P%i]%P;
    }
    inline void Add(int x,int w){
        g=1ll*g*Inv[c[x]+1]%P;
        g=1ll*g*((c[x]+=w)+1)%P;
    }
    inline void add(int x,int w){
        if(w>0) g=1ll*g*x%P;
        else g=1ll*g*Inv[x]%P;
        for(int i=1;i<=*p&&1ll*p[i]*p[i]<=x;i++)
            while(x%p[i]==0)
                Add(p[i],w),x/=p[i];
        if(x>1) Add(x,w);
    }
    #undef Maxx
}PT;
vector<int> Q[N],work;
paris q[N];

inline char nc(){
    static char buf[100000],*p1=buf,*p2=buf;
    if(p1==p2){ p2=(p1=buf)+fread(buf,1,100000,stdin); if(p1==p2) return EOF; }
    return *p1++;
}

inline void reaD(int &x){
    char c=nc(); x=0;
    for(;c>57||c<48;c=nc());for(;c>=48&&c<=57;x=x*10+c-48,c=nc());
}

inline void Insert(int x,int y){
    E[++cnt].t=y; E[cnt].nx=G[x]; G[x]=cnt;
    E[++cnt].t=x; E[cnt].nx=G[y]; G[y]=cnt;
}

int dfs(int x,int f){
    PT.add(a[x],1);
    int Size=1;
    for(int i=G[x];i;i=E[i].nx)if(E[i].t!=f) Size+=dfs(E[i].t,x);
    V[x]=Vector(PT.g,4*x,Size);
    PT.add(a[x],-1);
    if(PT.g==0)
        int s=0;
    Mat[x]=Matrix(V[x]);
    return Size;
}

int checking(int x,int f){
    int Size=1;
    for(int i=G[x];i;i=E[i].nx)
        if(E[i].t!=f&&!vis[E[i].t]) Size+=checking(E[i].t,x);
    return Size;
}

int explore(int x,int f){
    int Size=1,iMax=0;
    for(int i=G[x];i;i=E[i].nx)
        if(E[i].t!=f&&!vis[E[i].t]){
            int k=explore(E[i].t,x);
            iMax=max(iMax,k);
            Size+=k;
        }
    iMax=max(iMax,isize-Size);
    if(iMax<Max) Max=iMax,Root=x;
    return Size;
}

void dfs1(int x,int f,int par){
    pos[x]=par; dpt[x]=dpt[f]^1;
    A[x]=Mat[x]*A[f]; fa[x]=f;
    for(int i=G[x];i;i=E[i].nx)
        if(E[i].t!=f&&!vis[E[i].t]) dfs1(E[i].t,x,par);
}

void dfs2(int x,int f){
    B[x]=B[f]*Mat[x];
    for(int i=G[x];i;i=E[i].nx)
        if(E[i].t!=f&&!vis[E[i].t]) dfs2(E[i].t,x);
}

void solve(int x){
    isize=checking(x,0); Max=1<<30;
    explore(x,0); 
    int iRoot=Root;
    vis[iRoot]=1; pos[iRoot]=iRoot;
    for(int i=0;i<3;i++)
        for(int j=0;j<3;j++) A[iRoot][i][j]=0;
    A[iRoot][0][0]=A[iRoot][1][1]=A[iRoot][2][2]=1;
    B[iRoot]=A[iRoot]*Mat[iRoot]; dpt[iRoot]=0; fa[iRoot]=0;
    for(int i=G[iRoot];i;i=E[i].nx)
        if(!vis[E[i].t]) dfs1(E[i].t,iRoot,E[i].t),dfs2(E[i].t,iRoot);
    work=Q[x]; Q[x].clear();
    for(int i=0;i<(int)work.size();i++)
        if(pos[q[work[i]].first]==pos[q[work[i]].second]) Q[pos[q[work[i]].first]].push_back(work[i]);
        else{
            u=q[work[i]].first; v=q[work[i]].second;
            C=A[v]*B[fa[u]]; int x=work[i];
            Ans[work[i]][0]=(1ll*C[0][0]*V[u][0]+1ll*C[0][1]*V[u][1]+1ll*C[0][2]*V[u][2])%P;
            Ans[work[i]][1]=(1ll*C[1][0]*V[u][0]+1ll*C[1][1]*V[u][1]+1ll*C[1][2]*V[u][2])%P;
            Ans[work[i]][2]=(1ll*C[2][0]*V[u][0]+1ll*C[2][1]*V[u][1]+1ll*C[2][2]*V[u][2])%P;
            if(dpt[u]^dpt[v]) Ans[work[i]].reverse();
        }
    for(int i=G[iRoot];i;i=E[i].nx)
        if(!vis[E[i].t]) solve(E[i].t);
}

int main(){
    B[0][0][0]=B[0][1][1]=B[0][2][2]=1;
    reaD(n); reaD(m); PT.Init();
    for(int i=1;i<=n;i++) reaD(a[i]);
    for(int i=1;i<n;i++)
        reaD(u),reaD(v),Insert(u,v);
    dfs(1,0);
    for(int i=1;i<=m;i++){
        reaD(q[i].first); reaD(q[i].second);
        if(q[i].first==q[i].second){ Ans[i]=V[q[i].first]; continue; }
        Q[1].push_back(i);
    }
    solve(1);
    for(int i=1;i<=m;i++)
        printf("%d %d %d\n",Ans[i][0],Ans[i][1],Ans[i][2]);
    return 0;
}