[矩阵 点分治] BZOJ 4623 Styx

首先我们可以发现$g=(x*\phi)*1=(\phi*1)*x=x*x$ 所以 $g(n)=n*d_0(n)$ 其中 $d_0(n)$ 表示 n <script type="math/tex" id="MathJax-Element-267">n</script>的约数个数
然后就是树上的问题了
我们知道叉乘不满足结合律 打完之后才知道 汗 但是满足反交换律
然后我们就可以把叉乘表示成矩阵的形式 这是有结合律的
然后就是卡常历程 最后几乎就是照着Claris打了

#include<cstdio>
#include<cstdlib>
#include<algorithm>
#include<vector>
#include<cstring>
#define cl(x) memset(x,0,sizeof(x))
using namespace std;
typedef long long ll;

inline char nc(){
  static char buf[100000],*p1=buf,*p2=buf;
  return p1==p2&&(p2=(p1=buf)+fread(buf,1,100000,stdin),p1==p2)?EOF:*p1++;
}
inline void read(int &x){
  char c=nc(),b=1;
  for (;!(c>='0' && c<='9');c=nc()) if (c=='-') b=-1;
  for (x=0;c>='0' && c<='9';x=x*10+c-'0',c=nc()); x*=b;
}
inline void write(int x) {
  if (!x) return (void) putchar('0');
  if (x < 0) putchar('-'), x = -x;
  static short s[12], t;
  while (x) s[++t] = x % 10, x /= 10;
  while (t) putchar('0' + s[t--]);
}

const int P=1000000007;
inline ll Pow(ll a,int b){
  ll ret=1; for (;b;b>>=1,a=a*a%P) if (b&1) ret=ret*a%P; return ret;
}
inline ll Inv(ll a){
  return Pow(a,P-2);
}
namespace FG{
  const int maxn=1000001;
  int p[maxn+5],vst[maxn+5],num;
  int back[maxn+5];
  inline void Pre(){
    for (int i=2;i<=maxn;i++){
      if (!vst[i]) p[++num]=i,back[i]=num;
      for (int j=1;j<=num && (ll)i*p[j]<=maxn;j++){
    vst[i*p[j]]=1;
    if (i%p[j]==0) break;
      }
    }
  }
  int q[maxn+5]; ll sum=1;
  inline void add(int i,int j){
    sum=sum*Inv(q[i]+1)%P;
    q[i]+=j;
    sum=sum*(q[i]+1)%P;
  }
  inline void Add(int n,int r){
    int y=n;
    for (int i=1;i<=num && (ll)p[i]*p[i]<=n;i++)
      if (y%p[i]==0){
    int t=0; while (y%p[i]==0) t++,y/=p[i];
    add(i,t*r);
      }
    if (y!=1) add(back[y],r);
  }
}

const int N=100005;

struct edge{
  int u,v,next;
}G[N<<1];
int head[N],inum;
inline void add(int u,int v,int p){
  G[p].u=u; G[p].v=v; G[p].next=head[u]; head[u]=p;
}
#define V G[p].v

struct Mat{
  int a[3][3];
  int *operator[](int x){ return a[x]; }
  friend Mat operator * (Mat A,Mat B){
    Mat C; cl(C.a);
    for (int i=0;i<3;i++)
    for (int j=0;j<3;j++)
    for (int k=0;k<3;k++)
    (C[i][j]+=(ll)A[i][k]*B[k][j]%P)%=P;
    return C;
    }
  void read(int *x){
    a[0][0]=0; a[0][1]=P-x[2]; a[0][2]=x[1];
    a[1][0]=x[2]; a[1][1]=0; a[1][2]=P-x[0];
    a[2][0]=P-x[1];a[2][1]=x[0]; a[2][2]=0;
  }
};

inline void mul(Mat &a,Mat &b,Mat &c){
  c[0][0]=((ll)a[0][0]*b[0][0]+(ll)a[0][1]*b[1][0]+(ll)a[0][2]*b[2][0])%P;
  c[0][1]=((ll)a[0][0]*b[0][1]+(ll)a[0][1]*b[1][1]+(ll)a[0][2]*b[2][1])%P;
  c[0][2]=((ll)a[0][0]*b[0][2]+(ll)a[0][1]*b[1][2]+(ll)a[0][2]*b[2][2])%P;
  c[1][0]=((ll)a[1][0]*b[0][0]+(ll)a[1][1]*b[1][0]+(ll)a[1][2]*b[2][0])%P;
  c[1][1]=((ll)a[1][0]*b[0][1]+(ll)a[1][1]*b[1][1]+(ll)a[1][2]*b[2][1])%P;
  c[1][2]=((ll)a[1][0]*b[0][2]+(ll)a[1][1]*b[1][2]+(ll)a[1][2]*b[2][2])%P;
  c[2][0]=((ll)a[2][0]*b[0][0]+(ll)a[2][1]*b[1][0]+(ll)a[2][2]*b[2][0])%P;
  c[2][1]=((ll)a[2][0]*b[0][1]+(ll)a[2][1]*b[1][1]+(ll)a[2][2]*b[2][1])%P;
  c[2][2]=((ll)a[2][0]*b[0][2]+(ll)a[2][1]*b[1][2]+(ll)a[2][2]*b[2][2])%P;
}

int n,m;

int val[N],_val[N];
int size[N];

int a[N][3]; Mat Val[N];
inline void pdfs(int u,int fa){
  using namespace FG;
  size[u]=1; Add(val[u],1);
  _val[u]=(ll)val[u]*_val[fa]%P;
  a[u][0]=sum*_val[u]%P; a[u][1]=4*u;
  for (int p=head[u];p;p=G[p].next)
    if (V!=fa)
      pdfs(V,u),size[u]+=size[V];
  a[u][2]=size[u]; Add(val[u],-1);
  Val[u].read(a[u]);
}

Mat A[N],B[N],C;
int ans[N][3];

int del[N];
int sum,minv,rt;

inline void Root(int u,int fa){
  size[u]=1; int maxv=0;
  for (int p=head[u];p;p=G[p].next)
    if (V!=fa && !del[V])
      Root(V,u),size[u]+=size[V],maxv=max(maxv,size[V]);
  maxv=max(maxv,sum-size[u]);
  if (maxv<minv) minv=maxv,rt=u;
}

int depth[N],fat[N],col[N];
//int ncnt=0,_cnt[N];
int cnt;
void dfs(int u,int fa,int z){
  //++ncnt; _cnt[u]++;
  ++cnt;
  col[u]=z; fat[u]=fa; depth[u]=depth[fa]^1;
  mul(Val[u],A[fa],A[u]); mul(B[fa],Val[u],B[u]);
  //A[u]=Val[u]*A[fa]; B[u]=B[fa]*Val[u];
  for (int p=head[u];p;p=G[p].next)
    if (!del[V] && V!=fa)
      dfs(V,u,z);
}

vector<int> Vec[N];

int _u[N],_v[N];
int _rt[N];
inline void Divide(int u){
  if (Vec[u].size()==0) return;
  del[u]=1; col[u]=u; depth[u]=0; fat[u]=0;
  cl(A[u].a); cl(B[u].a);
  A[u][0][0]=A[u][1][1]=A[u][2][2]=1;
  for (int i=0;i<3;i++) for (int j=0;j<3;j++) B[u][i][j]=Val[u][i][j];
  for (int p=head[u];p;p=G[p].next)
    if (!del[V])
      cnt=0,dfs(V,u,V),size[V]=cnt;
  for (int p=head[u];p;p=G[p].next)
    if (!del[V]){
      sum=size[V],minv=1<<30;
      Root(V,u); _rt[V]=rt;
    }
  for (int i=0;i<(int)Vec[u].size();i++){
    int z=Vec[u][i],x=_u[z],y=_v[z];
    if (col[x]!=col[y]){
      //C=A[y]*B[fat[x]];
      mul(A[y],B[fat[x]],C);
      ans[z][0]=((ll)C[0][0]*a[x][0]+(ll)C[0][1]*a[x][1]+(ll)C[0][2]*a[x][2])%P;
      ans[z][1]=((ll)C[1][0]*a[x][0]+(ll)C[1][1]*a[x][1]+(ll)C[1][2]*a[x][2])%P;
      ans[z][2]=((ll)C[2][0]*a[x][0]+(ll)C[2][1]*a[x][1]+(ll)C[2][2]*a[x][2])%P;
      if (depth[x]^depth[y])
    ans[z][0]=(P-ans[z][0])%P,ans[z][1]=(P-ans[z][1])%P,ans[z][2]=(P-ans[z][2])%P;
    }else{
      Vec[_rt[col[x]]].push_back(Vec[u][i]);
    }      
  }
  for (int p=head[u];p;p=G[p].next)
    if (!del[V])
      Divide(_rt[V]);
}

int main(){
  int x,y;
  freopen("styx.in","r",stdin);
  freopen("styx.out","w",stdout);
  read(n); read(m); FG::Pre();
  cl(B[0].a); B[0][0][0]=B[0][1][1]=B[0][2][2]=1;
  for (int i=1;i<=n;i++) read(val[i]);
  for (int i=1;i<n;i++) read(x),read(y),add(x,y,++inum),add(y,x,++inum);
  _val[0]=1; pdfs(1,0);
  sum=n,minv=1<<30,Root(1,0);
  for (int i=1;i<=m;i++){
    read(x); read(y); _u[i]=x,_v[i]=y;
    if (x==y)
      for (int j=0;j<3;j++)
    ans[i][j]=a[x][j];
    else
      Vec[rt].push_back(i);
  }
  Divide(rt);
//  for (int i=1;i<=m;i++) printf("%d %d %d\n",ans[i][0],ans[i][1],ans[i][2]);
  for (int i=1;i<=m;i++) write(ans[i][0]),putchar(' '),write(ans[i][1]),putchar(' '),write(ans[i][2]),putchar('\n');
  //printf("%d\n",ncnt);
  //for (int i=1;i<=n;i++)
//  printf("%d\n",_cnt[i]);
  return 0;
}