莫比乌斯反演-题目

$M o b i u s$

　　题意概述：求$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 }

problem b

　　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 }

Zap

　　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 }

Gcd

　　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 }

Crash的数字表格

　　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 }

jzptab

　　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 }

Surprise me

---shzr