莫比乌斯反演-题目
理论部分的证明实在太多了,而且很多题的做法就是个结论,如果全都放到理论部分看完后再说题不就没什么意思了吗?所以又单开了一篇题目:
理论上来说可能应该从易到难,但是有的题确实没什么太大难度,就放在强化版后面作为双倍经验了。开始做题前在网上收集了不少题目,结果概括一句话题意后发现双倍经验题真不少。
那么先从一道简单的题目开始吧:
1.Problem b:https://www.lydsy.com/JudgeOnline/problem.php?id=2301
题意概述:求$sum_{i=a}^bsum_{j=c}^d [(i,j)=k]$.$n$组询问,$n,a,b,c,d,k<=50000$.
这个询问方式很显然可以转化为二维平面上满足条件的点的个数,所以只需要找到一种求
$sum_{i=1}^nsum_{j=1}^m [(i,j)=k]$
的做法,然后用二维前缀和容斥一下即可.
一般用到这种给定$gcd$求数对的题目,都可以考虑枚举$gcd$然后找两个互质的数分别乘进去,但是对于这道$n,m$不同的题目,无法使用欧拉函数,那么先化一化式子:
$sum_{i=1}^{frac{n}{k}}sum_{j=1}^{frac{m}{k}}[(i,j)=1]$
把下标换的好看一点:
$N=frac{n}{k} ,M=frac{m}{k} ,sum_{i=1}^{N}sum_{j=1}^{M}[(i,j)=1]$
这个东西显然是非常好求的...理论部分已经证明过了,这里直接拿来用。
$sum_{d=1}^{min{N,M}}mu(d)frac{N}{d}frac{M}{d}$
两个整除部分可以除法分块,如果对两个分式的除法分块有疑问可以参考“模积和”那道题。
1 # include <cstdio> 2 # include <iostream> 3 # define ll long long 4 # define R register int 5 6 using namespace std; 7 8 const int maxn=50000; 9 int k,n,a,b,c,d,mu[maxn+5],pri[maxn+5],vis[maxn+5],s[maxn+5],h; 10 ll ans; 11 12 ll ask (int n,int m) 13 { 14 ll g=0; 15 int i=1,l; 16 if(n==0||m==0) return 0; 17 while(i<=n) 18 { 19 if(i>m) break; 20 l=min(n/(n/i),min(n,m/(m/i))); 21 g=g+1LL*(s[l]-s[i-1])*(n/i)*(m/i); 22 i=l+1; 23 } 24 return g; 25 } 26 27 int main() 28 { 29 scanf("%d",&n); 30 mu[1]=1; 31 for (R i=2;i<=maxn;++i) 32 { 33 if(!vis[i]) 34 pri[++h]=i,mu[i]=-1; 35 for (R j=1;j<=h&&i*pri[j]<=maxn;++j) 36 { 37 vis[ i*pri[j] ]=1; 38 if(i%pri[j]==0) break; 39 mu[ i*pri[j] ]=-mu[i]; 40 } 41 } 42 for (R i=1;i<=maxn;++i) 43 s[i]=s[i-1]+mu[i]; 44 for (R i=1;i<=n;++i) 45 { 46 scanf("%d%d%d%d%d",&a,&b,&c,&d,&k); 47 if(a%k) a=a/k+1; 48 else a=a/k; 49 b=b/k; 50 if(c%k) c=c/k+1; 51 else c=c/k; 52 d=d/k; 53 ans=ask(b,d)-ask(b,c-1)-ask(a-1,d)+ask(a-1,c-1); 54 printf("%lld ",ans); 55 } 56 return 0; 57 }
2.Zap:https://www.lydsy.com/JudgeOnline/problem.php?id=1101
题意概述:上一道题去掉二维前缀和的弱化版.
1 # include <cstdio> 2 # include <iostream> 3 # define ll long long 4 # define R register int 5 6 using namespace std; 7 8 const int maxn=50000; 9 int k,n,a,b,c,d,mu[maxn+5],pri[maxn+5],vis[maxn+5],s[maxn+5],h; 10 ll ans; 11 12 ll ask (int n,int m) 13 { 14 ll g=0; 15 int i=1,l; 16 if(n==0||m==0) return 0; 17 while(i<=n) 18 { 19 if(i>m) break; 20 l=min(n/(n/i),min(n,m/(m/i))); 21 g=g+1LL*(s[l]-s[i-1])*(n/i)*(m/i); 22 i=l+1; 23 } 24 return g; 25 } 26 27 int main() 28 { 29 scanf("%d",&n); 30 mu[1]=1; 31 for (R i=2;i<=maxn;++i) 32 { 33 if(!vis[i]) 34 pri[++h]=i,mu[i]=-1; 35 for (R j=1;j<=h&&i*pri[j]<=maxn;++j) 36 { 37 vis[ i*pri[j] ]=1; 38 if(i%pri[j]==0) break; 39 mu[ i*pri[j] ]=-mu[i]; 40 } 41 } 42 for (R i=1;i<=maxn;++i) 43 s[i]=s[i-1]+mu[i]; 44 for (R i=1;i<=n;++i) 45 { 46 scanf("%d%d%d",&b,&d,&k); 47 b=b/k; 48 d=d/k; 49 ans=ask(b,d); 50 printf("%lld ",ans); 51 } 52 return 0; 53 }
3.YY的GCD:https://lydsy.com/JudgeOnline/problem.php?id=2820
4.Gcd:https://www.lydsy.com/JudgeOnline/problem.php?id=2818
题意概述:上一道题的单组询问+(n=m)版.
这道题有两种做法!
(1).直接将上一道题改成单组;
(2).将第一道题搬过来,枚举质数做多次;
是不是有种被欺骗的感觉...
1 # include <cstdio> 2 # include <iostream> 3 # define ll long long 4 # define R register int 5 6 using namespace std; 7 8 const int maxn=10000007; 9 int n,mu[maxn],pri[maxn],vis[maxn],h,s[maxn]; 10 ll ans=0; 11 12 ll ask (int n) 13 { 14 ll ans=0; 15 int i=1,l; 16 while(i<=n) 17 { 18 l=min(n/(n/i),n); 19 ans=ans+1LL*(s[l]-s[i-1])*(n/i)*(n/i); 20 i=l+1; 21 } 22 return ans; 23 } 24 25 int main() 26 { 27 scanf("%d",&n); 28 mu[1]=1; 29 for (R i=2;i<=n;++i) 30 { 31 if(!vis[i]) pri[++h]=i,mu[i]=-1; 32 for (R j=1;j<=h&&i*pri[j]<=n;++j) 33 { 34 vis[ i*pri[j] ]=1; 35 if(i%pri[j]==0) break; 36 mu[ i*pri[j] ]=-mu[i]; 37 } 38 } 39 for (R i=1;i<=n;++i) s[i]=s[i-1]+mu[i]; 40 for (R i=1;i<=h;++i) 41 ans=ans+ask(n/pri[i]); 42 printf("%lld",ans); 43 return 0; 44 }
5.于神之怒加强版:https://www.lydsy.com/JudgeOnline/problem.php?id=4407
题意概述:
$$sum_{i=1}^nsum_{j=1}^m(i,j)^k$$
$T<=2000,1<=N,M,K<=5000000$
做这道题的时候感觉好难啊,各种奇怪的化式子.做了几道题之后发现...全是套路.
$sum_{i=1}^nsum_{j=1}^m(i,j)^k$
枚举约数
$sum_{d=1}^{n}d^k sum_{i=1}^{N}sum_{j=1}^{M}[(i,j)=1]$
改成莫比乌斯函数形式
$sum_{d=1}^{n}d^k sum_{i=1}^{N}sum_{j=1}^{M}sum_{k|i,k|j}mu(k)$
交换合适顺序
$sum_{d=1}^n d^k sum_{k=1}^{n}mu(k) frac{n}{dk}frac{m}{dk}$
交换合式顺序+约数变倍数
$sum_{t=1}^{n} frac{n}{t}frac{m}{t}sum_{d|t} d^kmu(frac{t}{d})$
这时候可以将后面的部分单独拿出来作为一个新的函数:
$h(n)=sum_{d|n}d^kmu(frac{n}{d})$
这个函数看起来不像一个很好求的函数.
但是这里我们可以想到一种从第一篇介绍了之后就从未通入实际应用的科技:狄利克雷卷积.刚刚那个函数可以表示为下列两个函数的卷积.
$f(n)=n^k,g(n)=mu(frac{n}{d})$
显然这两个函数都是单调的.都是积性的.
那么两个积性函数的卷积也是积性的,而且积性函数都可以线性筛,这让我们非常愉悦.其实并不愉悦,因为线筛的方法看起来不大好找.
这里假设已经掌握了筛积性函数的成套方法,显然筛完求个前缀和,分一下块就做好了.筛函数的方法(这道题作为例题)
1 # include <cstdio> 2 # include <iostream> 3 # define ll long long 4 # define mod 1000000007 5 # define R register int 6 using namespace std; 7 8 const int maxn=5000000; 9 int T,n,m,k,pri[maxn+6],h,mu[maxn+6],vis[maxn+6],ls[maxn+6]; 10 ll f[maxn+6],g[maxn+6],s[maxn+6]; 11 12 ll qui (ll a,ll b) 13 { 14 ll s=1; 15 while(b) 16 { 17 if(b&1LL) s=s*a%mod; 18 a=a*a%mod; 19 b=b>>1LL; 20 } 21 return s%mod; 22 } 23 24 ll ask (int n,int m) 25 { 26 ll ans=0; 27 int i=1,l; 28 while(i<=n) 29 { 30 l=min(min(n,n/(n/i)),m/(m/i)); 31 ans=(ans+1LL*(g[l]-g[i-1]+mod)*(n/i)%mod*(m/i)%mod)%mod; 32 i=l+1; 33 } 34 return ans; 35 } 36 37 int main() 38 { 39 scanf("%d%d",&T,&k); 40 for (R i=1;i<=maxn;++i) 41 f[i]=qui(i,k); 42 mu[1]=1; 43 g[1]=1; 44 for (R i=2;i<=maxn;++i) 45 { 46 if(!vis[i]) mu[i]=-1,pri[++h]=i,g[i]=(f[i]-1+mod)%mod,ls[i]=i; 47 for (R j=1;j<=h&&i*pri[j]<=maxn;++j) 48 { 49 vis[ i*pri[j] ]=1; 50 ls[ i*pri[j] ]=pri[j]; 51 if(i%pri[j]==0) 52 { 53 g[ i*pri[j] ]=g[i]*f[ pri[j] ]%mod; 54 break; 55 } 56 mu[ i*pri[j] ]=-mu[i]; 57 g[ i*pri[j] ]=(g[i]*g[ pri[j] ])%mod; 58 } 59 } 60 for (R i=1;i<=maxn;++i) 61 g[i]=(g[i]+g[i-1])%mod; 62 while(T--) 63 { 64 scanf("%d%d",&n,&m); 65 if(n>m) swap(n,m); 66 printf("%lld ",ask(n,m)); 67 } 68 return 0; 69 }
6.Crash的数字表格:https://www.lydsy.com/JudgeOnline/problem.php?id=2154
题意概述:求$$sum_{i=1}^nsum_{j=1}^m [n,m],n,m<=10^7$$.
其实还是比较套路的...
根据最小公倍数的性质:
$sum_{i=1}^n sum_{j=1}^mfrac{ i imes j}{(i,j)}$
套路1:枚举$gcd$.
$sum_{d=1}^n sum_{i=1}^n sum_{j=1}^m[(i,j)==d] frac{ i imes j}{d}$
套路2:变成$gcd=1$的形式.
$sum_{d=1}^n sum_{i=1}^N sum_{j=1}^M[(i,j)==1] { i imes j imes d}$
套路3:变成$sum_{d|i}mu(d)$的形式.
$sum_{d=1}^n d sum_{i=1}^N sum_{j=1}^M{ i imes j}sum_{k|i,k|j} mu(k) $
套路4:枚举约数改为枚举倍数,交换求和次序.
令$i=xk,j=yk$
$sum_{d=1}^n d sum_{k=1}^n mu(k) sum_{x=1}^{N'} sum_{y=1}^{M'}{ x imes y imes k^2} $
$sum_{d=1}^n d sum_{k=1}^n mu(k) k^2 sum_{x=1}^{N'} sum_{y=1}^{M'}{ x imes y } $
设$s(i)=sum_{i=1}^n i$
$sum_{d=1}^n d sum_{k=1}^n mu(k) k^2 s(frac{m}{dk}) imes s(frac{n}{dk})$
对于$d$首先进行数论分块,内部再对$k$分块,就可以通过本题.
1 // luogu-judger-enable-o2 2 # include <cstdio> 3 # include <iostream> 4 # define R register int 5 # define ll long long 6 # define mod 20101009 7 8 using namespace std; 9 10 const int maxn=10000000; 11 int T,n,m,pri[maxn+7],mu[maxn+7],h,vis[maxn+7],s[maxn+7]; 12 13 void init() 14 { 15 mu[1]=1; 16 s[1]=1; 17 for (R i=2;i<=m;++i) 18 { 19 s[i]=(s[i-1]+i)%mod; 20 if(!vis[i]) mu[i]=-1,pri[++h]=i; 21 for (R j=1;j<=h&&i*pri[j]<=m;++j) 22 { 23 vis[ i*pri[j] ]=true; 24 if(i%pri[j]==0) break; 25 mu[ i*pri[j] ]=-mu[i]; 26 } 27 } 28 for (R i=1;i<=m;++i) 29 { 30 mu[i]=(1LL*mu[i]*i*i%mod+mod)%mod; 31 mu[i]=(mu[i]+mu[i-1]+mod)%mod; 32 } 33 } 34 35 ll mul (int n,int m) 36 { 37 ll ans=0; 38 int i=1,l; 39 if(n==0||m==0) return 0; 40 while (i<=n) 41 { 42 l=min(n,min(n/(n/i),m/(m/i))); 43 ans=(ans+1LL*((mu[l]-mu[i-1])%mod+mod)%mod*s[n/i]%mod*s[m/i]%mod)%mod; 44 i=l+1; 45 } 46 return ans; 47 } 48 49 ll ask (int n,int m) 50 { 51 ll ans=0; 52 int i=1,l; 53 if(n==0||m==0) return 0; 54 while(i<=n) 55 { 56 l=min(n,min(n/(n/i),m/(m/i))); 57 ans=(ans+1LL*(s[l]-s[i-1]+mod)*mul(n/i,m/i))%mod; 58 i=l+1; 59 } 60 return ans; 61 } 62 63 int main() 64 { 65 scanf("%d%d",&n,&m); 66 if(n>m) swap(n,m); 67 init(); 68 printf("%lld",ask(n,m)); 69 return 0; 70 }
7.jzptab:https://www.lydsy.com/JudgeOnline/problem.php?id=2693
题意概述:同上,多组询问.
上边那个看起来肯定是要再化一化了,其实还是套路啊,这种看到两个数相乘的一般就是要改成一个数,然后某一个变量改为枚举约数,另一个直接算,然后线性筛.前半部分除法分块,后边用前缀和优化.这道题就非常的套路.
$sum_{d=1}^n d sum_{k=1}^n mu(k) k^2 s(frac{m}{dk}) imes s(frac{n}{dk})$
$t=dk$
$sum_{t=1}^nts(frac{m}{t})s(frac{n}{t})sum_{k|t}k mu(k)$
$f(n)=nsum_{d|n}dmu(d)$
现在的问题是$f$怎么筛,这个函数显然可以按照线性筛的套路来做,每个质因子出现多次时没有任何意义,所以有了奇妙的代码:
1 if(i%pri[j]==0) 2 f[ i*pri[j] ]=f[i];
1 # include <cstdio> 2 # include <iostream> 3 # define R register int 4 # define ll long long 5 # define mod 100000009 6 7 using namespace std; 8 9 const int maxn=10000000; 10 int T,n,m,pri[maxn+7],mu[maxn+7],h,vis[maxn+7],s[maxn+7],f[maxn+7]; 11 12 void init (int n) 13 { 14 mu[1]=s[1]=f[1]=1; 15 for (R i=2;i<=n;++i) 16 { 17 s[i]=(s[i-1]+i)%mod; 18 if(!vis[i]) mu[i]=-1,pri[++h]=i,f[i]=(1-i+mod)%mod; 19 for (R j=1;j<=h&&i*pri[j]<=n;++j) 20 { 21 vis[ i*pri[j] ]=true; 22 if(i%pri[j]==0) 23 { 24 f[ i*pri[j] ]=f[i]; 25 break; 26 } 27 mu[ i*pri[j] ]=-mu[i]; 28 f[ i*pri[j] ]=1LL*f[i]*f[ pri[j] ]%mod; 29 } 30 } 31 for (R i=1;i<=n;++i) 32 f[i]=(1LL*f[i]*i+f[i-1])%mod; 33 } 34 35 ll ask (int n,int m) 36 { 37 ll ans=0; 38 int l,i=1; 39 while(i<=n) 40 { 41 l=min(n,min(n/(n/i),m/(m/i))); 42 ans=(ans+1LL*(f[l]-f[i-1]+mod)%mod*s[n/i]%mod*s[m/i]%mod)%mod; 43 i=l+1; 44 } 45 return ans; 46 } 47 48 int main() 49 { 50 scanf("%d",&T); 51 init(10000000); 52 while(T--) 53 { 54 scanf("%d%d",&n,&m); 55 if(n>m) swap(n,m); 56 printf("%lld ",ask(n,m)); 57 } 58 return 0; 59 }
8.数字表格:https://lydsy.com/JudgeOnline/problem.php?id=4816
题意概述:多组询问,$T<=1000,1<=n,m<=10^6$
$f(n)=left{egin{matrix}0,n=0 \ 1,n=1 \f(n-1)+f(n-2),others end{matrix} ight. $
$$prod_{i=1}^nprod_{j=1}^mf[(i,j)]$$
以前做过的反演都是加法类的,这道题很新奇.
其实乘法甚至更简单,因为可以什么都不考虑地随意交换乘法顺序,便于化式子.
枚举gcd:$prod_{d=1}^nf(d)^{g(d)}$
$g(d)=sum_{i=1}^nsum_{j=1}^m[(i,j)==d]$
$g(d)=sum_{i=1}^nmu(i)frac{n}{id}frac{m}{id}$
对于连续的一段$n,m,g$的值是相同的,利用这个性质进行分块,求$g$的时候内部再分块,单组询问可以做到$O(N)$
$prod_{d=1}^nf(d)^{sum_{i=1}^nmu(i)frac{n}{id}frac{m}{id}}$
枚举id:$prod_{T=1}^nprod_{d|T}f(d)^{mu(frac{T}{d})frac{n}{t}frac{m}{t}}$
$prod_{T=1}^n(prod_{d|T}f(d)^{mu(frac{T}{d})})^{frac{n}{t}frac{m}{t}}$
里面那个东西线性筛好像有点困难,但是可以枚举$d$来计算.
运用调和级数计算一下会发现预处理的复杂度不大.
1 # include <cstdio> 2 # include <iostream> 3 # define R register int 4 # define mod 1000000007 5 # define ll long long 6 #define getchar() (S==T&&(T=(S=BB)+fread(BB,1,1<<15,stdin),S==T)?EOF:*S++) 7 8 using namespace std; 9 10 const int maxn=1000006; 11 char BB[1 << 18], *S = BB, *T = BB; 12 int mp[maxn],num,n,m,vis[maxn],mu[maxn],f[maxn],pri[maxn],h,F[maxn],inv[maxn]; 13 int N[maxn],M[maxn],maxx; 14 ll ans=1,t=1; 15 16 inline int read() 17 { 18 R x=0; 19 register char c=getchar(); 20 while(!isdigit(c)) c=getchar(); 21 while(isdigit(c)) x=(x<<3)+(x<<1)+(c^48),c=getchar(); 22 return x; 23 } 24 25 inline void init (int n) 26 { 27 mu[1]=mp[1]=F[1]=1; 28 for (R i=2;i<=n;++i) 29 { 30 if(!vis[i]) pri[++h]=i,mu[i]=-1,mp[i]=i; 31 for (R j=1;j<=h&&i*pri[j]<=n;++j) 32 { 33 vis[ i*pri[j] ]=1; 34 if(i%pri[j]==0) break; 35 mu[ i*pri[j] ]=-mu[i]; 36 } 37 } 38 for (R i=1;i<=n;++i) 39 for (R j=1;i*j<=n;++j) 40 { 41 if(mu[j]==1) F[ i*j ]=1LL*F[ i*j ]*f[i]%mod; 42 else if(mu[j]==-1) F[ i*j ]=1LL*F[ i*j ]*inv[i]%mod; 43 } 44 } 45 46 inline int qui (int a,int b) 47 { 48 if(!a) return 1; 49 int s=1; 50 while(b) 51 { 52 if(b&1) s=1LL*s*a%mod; 53 a=1LL*a*a%mod; 54 b>>=1; 55 } 56 return s; 57 } 58 59 inline int ask (int x) 60 { 61 int ans=1; 62 for (R i=1;i<=x;++i) 63 { 64 if(x%i) continue; 65 if(mu[x/i]==1) ans=1LL*ans*f[i]%mod; 66 else if(mu[x/i]==-1) ans=1LL*ans*inv[i]%mod; 67 } 68 return ans; 69 } 70 71 inline int ad (int a,int b) { a+=b; if(a>=mod) a-=mod; return a; } 72 73 int main() 74 { 75 num=read(); 76 for (R i=1;i<=num;++i) 77 { 78 N[i]=read(),M[i]=read(); 79 if(M[i]<N[i]) swap(N[i],M[i]); 80 maxx=max(maxx,N[i]); 81 } 82 f[1]=1; 83 inv[1]=1; 84 F[1]=1; 85 for (R i=2;i<=maxx;++i) 86 f[i]=ad(f[i-1],f[i-2]),inv[i]=qui(f[i],mod-2),F[i]=1; 87 init(maxx); 88 F[0]=1; 89 for (R i=1;i<=maxx;++i) F[i]=1LL*F[i]*F[i-1]%mod; 90 for (R i=1;i<=num;++i) 91 { 92 n=N[i],m=M[i]; 93 ans=1; 94 int l=1,r=1; 95 while(l<=n) 96 { 97 r=min(n,min(n/(n/l),m/(m/l))); 98 t=1LL*F[r]*qui(F[l-1],mod-2)%mod; 99 ans=ans*qui(t,1LL*(n/l)*(m/l)%(mod-1))%mod; 100 l=r+1; 101 } 102 printf("%lld ",ans); 103 } 104 return 0; 105 }
9.数表:https://lydsy.com/JudgeOnline/problem.php?id=3529
题意概述:求以下式子的值,多组询问.$T<=2 imes 10^4$,$n,m<=10^5$ $$sum_{i=1}^nsum_{j=1}^msigma( quad (i,j) quad)[ quadsigma(i,j)<=a quad]$$
首先考虑一下怎么求约数和函数,一种做法是枚举倍数算贡献,非常快,但是...如果想再优化一下也是可以的,见线性筛一文.
10.完全平方数:https://www.lydsy.com/JudgeOnline/problem.php?id=2440
11.约数个数和:https://www.lydsy.com/JudgeOnline/problem.php?id=3994
12.Surprise Me:http://codeforces.com/problemset/problem/809/E
题意概述:给定一棵 $n$ 个点的树,每个点有一个权值,保证其构成一个 $1-n$ 的排列,求$frac{1}{n(n-1)}sumlimits_{i=1}^nsumlimits_{j=1}^nvarphi(a_i imes a_j)dist(i,j)$的值 $\%1e9+7$.$n<=10^5$
发现两个毫无关系的数乘在一起再取phi是非常困难的,所以定义$a_{p_i}=i$,将式子变为:
$frac{1}{n(n-1)}sumlimits_{i=1}^nsumlimits_{j=1}^nvarphi(i imes j)dist(p_i,p_j)$
这样就可以愉快的化式子了。
$frac{1}{n(n-1)}sumlimits_{i=1}^nsumlimits_{j=1}^nfrac{varphi(i)varphi(j)d}{varphi(d)}dist(p_i,p_j)[d=(i,j)]$
$frac{1}{n(n-1)}sumlimits_{d=1}^n frac{d}{varphi(d)}sumlimits_{i=1}^nsumlimits_{j=1}^nvarphi(i)varphi(j)dist(p_i,p_j)[(i,j)=d]$
$frac{1}{n(n-1)}sumlimits_{d=1}^n frac{d}{varphi(d)}sumlimits_{i=1}^Nsumlimits_{j=1}^Nvarphi(id)varphi(jd)dist(p_{id},p_{jd})[(i,j)=1]$
$frac{1}{n(n-1)}sumlimits_{d=1}^n frac{d}{varphi(d)}sumlimits_{i=1}^Nsumlimits_{j=1}^Nvarphi(id)varphi(jd)dist(p_{id},p_{jd})sumlimits_{k|i,k|j}mu(k)$
$frac{1}{n(n-1)}sumlimits_{d=1}^n frac{d}{varphi(d)}sumlimits_{k=1}^Nmu(k)sumlimits_{i=1}^Nsumlimits_{j=1}^Nvarphi(ikd)varphi(jkd)dist(p_{ikd},p_{jkd})$
$frac{1}{n(n-1)}sumlimits_{T=1}^nsumlimits_{i=1}^nsumlimits_{j=1}^nvarphi(iT)varphi(jT)dist(p_{iT},p_{jT})sumlimits_{k|T}mu(frac{T}{k})frac{k}{varphi(k)}$
$frac{1}{n(n-1)}sumlimits_{T=1}^nf(T)sumlimits_{k|T}mu(frac{T}{k})frac{k}{varphi(k)}$
$f(x)=sumlimits_{i=1}^{n/T}sumlimits_{j=1}^{n/T}varphi(iT)varphi(jT)dist(p_{iT},p_{jT})$
那么f怎么求呢?首先将所有满足要求的点拿出来,要求的就是两两点对间的一些信息,可以建出虚树来做。对于任意两个点,如果有祖孙关系,可以在dfs时很方便的求出来;如果没有,考虑在lca处合并信息时求出。如果用 $d(x)$ 表示从 $x$ 到当前祖先的距离,那么合并后的信息应为 $varphi(x)varphi(y)(d(x)+d(y))$,只要在dfs过程中维护 $sum varphi(x)$ 以及 $sum varphi(x)d(x)$ 即可。复杂度的瓶颈其实在求LCA,如果使用树剖,可以得到一个小常数的两个log做法,如果使用更快的方法,复杂度则为 $nlogn$
1 # include <cstdio> 2 # include <iostream> 3 # include <cstring> 4 # include <algorithm> 5 # define ll long long 6 7 using namespace std; 8 9 const int maxn=200005; 10 const int mod=1000000007; 11 int n,v[maxn],h,firs[maxn],x,y,p[maxn],mu[maxn],phi[maxn],vis[maxn],p_cnt,pri[maxn],F[maxn]; 12 int siz[maxn],son[maxn],Top[maxn],f[maxn],ans,dep[maxn],inv_phi[maxn],a[maxn],m,id[maxn],id_cnt; 13 ll dp[maxn],sta[maxn],vv[maxn],va[maxn]; 14 struct edge { int too,nex,c; }; 15 struct gragh 16 { 17 edge g[maxn<<1]; 18 int h,firs[maxn]; 19 void add (int x,int y,int c) 20 { 21 g[++h].nex=firs[x]; 22 firs[x]=h; 23 g[h].too=y; 24 g[h].c=c; 25 g[++h].nex=firs[y]; 26 firs[y]=h; 27 g[h].too=x; 28 g[h].c=c; 29 } 30 }g1,g; 31 32 int qui (int a,int b) 33 { 34 int s=1; 35 while(b) 36 { 37 if(b&1) s=1LL*s*a%mod; 38 a=1LL*a*a%mod; 39 b>>=1; 40 } 41 return s; 42 } 43 44 bool cmp (int a,int b) { return id[a]<id[b]; } 45 46 void init() 47 { 48 mu[1]=1; phi[1]=1; 49 for (int i=2;i<=n;++i) 50 { 51 if(!vis[i]) mu[i]=-1,phi[i]=i-1,pri[++p_cnt]=i; 52 for (int j=1;j<=p_cnt&&i*pri[j]<=n;++j) 53 { 54 vis[ i*pri[j] ]=1; 55 if(i%pri[j]==0) { phi[ i*pri[j] ]=phi[i]*pri[j]; break; } 56 phi[ i*pri[j] ]=phi[i]*phi[ pri[j] ]; 57 mu[ i*pri[j] ]=-mu[i]; 58 } 59 } 60 for (int i=1;i<=n;++i) inv_phi[i]=qui(phi[i],mod-2); 61 for (int i=1;i<=n;++i) 62 for (int j=i;j<=n;j+=i) 63 F[j]=(F[j]+1LL*i*inv_phi[i]*mu[j/i]%mod+mod)%mod; 64 } 65 66 void dfs1 (int x) 67 { 68 int j,maxs=-1; 69 siz[x]=1; 70 for (int i=g1.firs[x];i;i=g1.g[i].nex) 71 { 72 j=g1.g[i].too; 73 if(dep[j]) continue; 74 dep[j]=dep[x]+1; f[j]=x; 75 dfs1(j); 76 siz[x]+=siz[j]; 77 if(siz[j]>=maxs) maxs=siz[j],son[x]=j; 78 } 79 } 80 81 void dfs2 (int x,int Tp) 82 { 83 Top[x]=Tp; id[x]=++id_cnt; 84 if(!son[x]) return; 85 dfs2(son[x],Tp); 86 int j; 87 for (int i=g1.firs[x];i;i=g1.g[i].nex) 88 { 89 j=g1.g[i].too; 90 if(j==f[x]||j==son[x]) continue; 91 dfs2(j,j); 92 } 93 } 94 95 int lca (int x,int y) 96 { 97 while(Top[x]!=Top[y]) 98 { 99 if(dep[ Top[x] ]>dep[ Top[y] ]) swap(x,y); 100 y=f[ Top[y] ]; 101 } 102 return (dep[x]<dep[y])?x:y; 103 } 104 105 int d (int x,int y) { return dep[x]+dep[y]-2*dep[lca(x,y)]; } 106 107 void DP (int x) 108 { 109 int j,vl,t1=0,t2=0,len; dp[x]=0; va[x]=0; 110 vv[x]=vl=vis[x]?v[x]:0; 111 for (int i=g.firs[x];i;i=g.g[i].nex) 112 { 113 j=g.g[i].too; 114 if(dep[j]<dep[x]) continue; 115 DP(j); 116 len=g.g[i].c; 117 dp[x]=(dp[x]+dp[j])%mod; 118 dp[x]=(dp[x]+1LL*vl*(vv[j]*len%mod+va[j]))%mod; 119 dp[x]=(dp[x]+1LL*t2*vv[j]+1LL*t1*(vv[j]*len%mod+va[j]))%mod; 120 vv[x]=(vv[x]+vv[j])%mod; 121 va[x]=(va[x]+va[j]+vv[j]*len%mod)%mod; 122 t1=(t1+vv[j])%mod; t2=(t2+va[j]+vv[j]*len%mod)%mod; 123 } 124 } 125 126 int ask (int x) 127 { 128 int m=0,Tp=0; 129 for (int i=1;i<=n/x;i++) 130 a[++m]=p[i*x]; 131 g.h=0; 132 sort(a+1,a+1+m,cmp); 133 g.firs[1]=0; 134 for (int i=1;i<=m;++i) g.firs[ a[i] ]=0,vis[ a[i] ]=1; 135 for (int i=2;i<=m;++i) g.firs[ lca(a[i],a[i-1]) ]=0; 136 sta[++Tp]=1; 137 int st=1; if(a[1]==1) st=2; 138 for (int i=st;i<=m;++i) 139 { 140 int LCA=lca(sta[Tp],a[i]); 141 while(Tp>1&&LCA!=sta[Tp]) 142 { 143 if(sta[Tp-1]==LCA) { g.add(sta[Tp],sta[Tp-1],dep[ sta[Tp] ]-dep[ sta[Tp-1] ]); Tp--; break;} 144 if(id[ sta[Tp-1] ]<id[ LCA ]) { g.add(LCA,sta[Tp],dep[ sta[Tp] ]-dep[ LCA ]); sta[Tp]=LCA; break; } 145 g.add(sta[Tp],sta[Tp-1],dep[ sta[Tp] ]-dep[ sta[Tp-1] ]); Tp--; 146 } 147 sta[++Tp]=a[i]; 148 } 149 while(Tp>1) g.add(sta[Tp],sta[Tp-1],dep[ sta[Tp] ]-dep[ sta[Tp-1] ]),Tp--; 150 DP(1); 151 for (int i=1;i<=m;++i) vis[ a[i] ]=0; 152 return dp[1]*2%mod; 153 } 154 155 inline int read() 156 { 157 int x=0; 158 char c=getchar(); 159 while (!isdigit(c)) c=getchar(); 160 while (isdigit(c)) x=(x<<3)+(x<<1)+(c^48),c=getchar(); 161 return x; 162 } 163 164 int main() 165 { 166 scanf("%d",&n); init(); 167 for (int i=1;i<=n;++i) v[i]=read(),p[ v[i] ]=i,v[i]=phi[ v[i] ]; 168 for (int i=1;i<n;++i) 169 { 170 x=read(),y=read(); 171 g1.add(x,y,1); 172 } 173 dep[1]=1; dfs1(1); dfs2(1,1); 174 memset(vis,0,sizeof(vis)); 175 for (int i=1;i<=n;++i) 176 ans=(ans+1LL*F[i]*ask(i)%mod)%mod; 177 ans=1LL*ans*qui(1LL*n*(n-1)%mod,mod-2)%mod; 178 printf("%d ",ans); 179 return 0; 180 }
---shzr