[矩阵 点分治] BZOJ 4623 Styx

首先我们可以发现g=(x∗ϕ)∗1=(ϕ∗1)∗x=x∗x 所以 g(n)=n∗d0(n) 其中 d0(n) 表示n的约数个数
然后就是树上的问题了
我们知道叉乘不满足结合律 打完之后才知道 汗 但是满足反交换律
然后我们就可以把叉乘表示成矩阵的形式 这是有结合律的
然后就是卡常历程 最后几乎就是照着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].vstruct 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;}