BZOJ4623 : Styx
g是积性函数,可以通过分解质因数在O(nlognloglogn)的时间内求出。
对于((A×B)×C)×D,可以转化为D×(C×(B×A)),并视向量个数的奇偶性取反答案。
对于D×(C×(B×A)),可以将D×,C×,B×用3个3×3的矩阵表示,然后对树进行点分治即可。
时间复杂度O(nlognloglogn)。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 | #include<cstdio> const int N=100010,E=N*20,M=1000010,P=1000000007; int n,m,i,j,x,y,a[N],mv; int g[N],v[N<<1],nxt[N<<1],ok[N<<1],ed; int all,f[N],son[N],now,pos[N]; int G[N],V[E],NXT[E],ED,p[N],d[N],cnt,ans[N][3]; struct Q{ int x,y;}q[N]; inline void read( int &a){ char c; while (!(((c= getchar ())>= '0' )&&(c<= '9' )));a=c- '0' ; while (((c= getchar ())>= '0' )&&(c<= '9' ))(a*=10)+=c- '0' ;} inline void add( int x, int y){v[++ed]=y;nxt[ed]=g[x];ok[ed]=1;g[x]=ed;} inline void ADD( int x, int y){V[++ED]=y;NXT[ED]=G[x];G[x]=ED;} namespace Num{ int inv[M],p[N],tot,mul,m[M],pm[M],g[M]; bool v[M]; void init(){ for (i=2;i<=mv;i++){ if (!v[i])p[tot++]=i; for (j=0;j<tot;j++){ if (i*p[j]>mv) break ; v[i*p[j]]=1; if (i%p[j]==0) break ; } } for (inv[0]=inv[1]=1,i=2;i<=mv;i++)inv[i]=1LL*(P-inv[P%i])*(P/i)%P; for (mul=1,i=2;i<=mv;i++)pm[i]=g[i]=1; } int exgcd( int a, int b){ if (!b) return x=1,y=0,a; int d=exgcd(b,a%b),t=x; return x=y,y=t-a/b*y,d; } inline int rev( int a){exgcd(a,P); while (x<0)x+=P; return x%P;} inline void up( int &a, int b){a+=b; if (a>=P)a-=P;} inline void add( int x){ m[x]++; g[x]=(1LL*(x-1)*pm[x]%P*m[x]+g[x])%P; pm[x]=1LL*pm[x]*x%P; up(g[x],pm[x]); } inline void del( int x){ up(g[x],P-pm[x]); pm[x]=1LL*pm[x]*inv[x]%P; g[x]=(1LL*(P-x+1)*pm[x]%P*m[x]+g[x])%P; m[x]--; } inline void divadd( int n){ for ( int i=0;p[i]*p[i]<=n;i++) if (n%p[i]==0){ int j=p[i],old=rev(g[j]); while (n%j==0)n/=j,add(j); mul=1LL*mul*old%P*g[j]%P; } if (n==1) return ; int old=rev(g[n]); add(n); mul=1LL*mul*old%P*g[n]%P; } inline void divdel( int n){ for ( int i=0;p[i]*p[i]<=n;i++) if (n%p[i]==0){ int j=p[i],old=rev(g[j]); while (n%j==0)n/=j,del(j); mul=1LL*mul*old%P*g[j]%P; } if (n==1) return ; int old=rev(g[n]); del(n); mul=1LL*mul*old%P*g[n]%P; } } int mat[N][3][3],vec[N][3],A[N][3][3],B[N][3][3],C[3][3]; inline void mul( int a[][3], int b[][3], int c[][3]){ 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; } void dfs( int x, int y){ Num::divadd(a[x]); vec[x][0]=Num::mul; vec[x][1]=x<<2; vec[x][2]=1; for ( int i=g[x];i;i=nxt[i]) if (v[i]!=y)dfs(v[i],x),vec[x][2]+=vec[v[i]][2]; mat[x][0][1]=P-vec[x][2]; mat[x][0][2]=vec[x][1]; mat[x][1][0]=vec[x][2]; mat[x][1][2]=P-vec[x][0]; mat[x][2][0]=P-vec[x][1]; mat[x][2][1]=vec[x][0]; Num::divdel(a[x]); } void findroot( int x, int y){ son[x]=1;f[x]=0; for ( int i=g[x];i;i=nxt[i]) if (ok[i]&&v[i]!=y){ findroot(v[i],x); son[x]+=son[v[i]]; if (son[v[i]]>f[x])f[x]=son[v[i]]; } if (all-son[x]>f[x])f[x]=all-son[x]; if (f[x]<f[now])now=x; } void dfs1( int x, int y, int z){ pos[x]=z;f[x]=y;d[x]=d[y]^1; mul(mat[x],A[y],A[x]); for ( int i=g[x];i;i=nxt[i]) if (ok[i]&&v[i]!=y)dfs1(v[i],x,z); } void dfs2( int x, int y){ mul(B[y],mat[x],B[x]); for ( int i=g[x];i;i=nxt[i]) if (ok[i]&&v[i]!=y)dfs2(v[i],x); } void solve( int x){ if (!G[x]) return ; f[0]=all=son[x],findroot(x,now=0); int i; pos[now]=now; for ( int i=0;i<3;i++) for ( int j=0;j<3;j++)A[now][i][j]=0; A[now][0][0]=A[now][1][1]=A[now][2][2]=1; mul(A[now],mat[now],B[now]); f[now]=d[now]=0; for (i=g[now];i;i=nxt[i]) if (ok[i])dfs1(v[i],now,v[i]),dfs2(v[i],now); for (cnt=0,i=G[x];i;i=NXT[i])p[++cnt]=V[i];G[x]=0; for (i=1;i<=cnt;i++) if (pos[q[p[i]].x]==pos[q[p[i]].y])ADD(pos[q[p[i]].x],p[i]); else { int X=q[p[i]].x,Y=q[p[i]].y; mul(A[Y],B[f[X]],C); ans[p[i]][0]=(1LL*C[0][0]*vec[X][0]+1LL*C[0][1]*vec[X][1]+1LL*C[0][2]*vec[X][2])%P; ans[p[i]][1]=(1LL*C[1][0]*vec[X][0]+1LL*C[1][1]*vec[X][1]+1LL*C[1][2]*vec[X][2])%P; ans[p[i]][2]=(1LL*C[2][0]*vec[X][0]+1LL*C[2][1]*vec[X][1]+1LL*C[2][2]*vec[X][2])%P; if (d[X]^d[Y]){ ans[p[i]][0]=P-ans[p[i]][0]; ans[p[i]][1]=P-ans[p[i]][1]; ans[p[i]][2]=P-ans[p[i]][2]; } } for (i=g[now];i;i=nxt[i]) if (ok[i])ok[i^1]=0,solve(v[i]); } int main(){ B[0][0][0]=B[0][1][1]=B[0][2][2]=1; read(n),read(m); for (ed=i=1;i<=n;i++){ read(a[i]); if (a[i]>mv)mv=a[i]; } for (i=1;i<n;i++)read(x),read(y),add(x,y),add(y,x); Num::init(); dfs(1,0); for (i=1;i<=m;i++){ read(q[i].x),read(q[i].y); if (q[i].x==q[i].y){ ans[i][0]=vec[q[i].x][0]; ans[i][1]=vec[q[i].x][1]; ans[i][2]=vec[q[i].x][2]; } else ADD(1,i); } son[1]=n;solve(1); for (i=1;i<=m;i++) printf ( "%d %d %d\n" ,ans[i][0],ans[i][1],ans[i][2]); return 0; } |
【推荐】编程新体验,更懂你的AI,立即体验豆包MarsCode编程助手
【推荐】博客园携手 AI 驱动开发工具商 Chat2DB 推出联合终身会员
【推荐】抖音旗下AI助手豆包,你的智能百科全书,全免费不限次数
【推荐】轻量又高性能的 SSH 工具 IShell:AI 加持,快人一步
· 一个超经典 WinForm,WPF 卡死问题的终极反思
· ASP.NET Core - 日志记录系统(二)
· .NET 依赖注入中的 Captive Dependency
· .NET Core 对象分配(Alloc)底层原理浅谈
· 聊一聊 C#异步 任务延续的三种底层玩法
· 互联网不景气了那就玩玩嵌入式吧,用纯.NET开发并制作一个智能桌面机器人(一):从.NET IoT入
· .NET 开发的分流抢票软件,不做广告、不收集隐私
· 一个超经典 WinForm,WPF 卡死问题的终极反思
· 开箱你的 AI 语音女友「GitHub 热点速览」
· 前端实现 HTML 网页转 PDF 并导出