题意
给定两个长度为 $n$ 的序列 $a,b$。
对于一个 $1$ 到 $n$ 的排列 $p$，记 $c_i=gcd(a_i,b_{p_i})$，$sigma(c)$ 表示序列 $c$ 中所有元素的方差。
求 $$sumlimits_{p}sigma(c)$$ 对 $10^9+7$ 取模。
$2leq n,a_i,b_ileq 10^6$ 。
题解
方差可以写成以下形式
$$egin{aligned}sigma(c)&=dfrac{1}{n}sum_{i=1}^n (c_i-ar{c})^2\&=dfrac{1}{n}sum_{i=1}^n c_i^2-dfrac{1}{n^2}(sum_{i=1}^n c_i)^2end{aligned}$$
那么当前问题是求出 $sum_p sum_{i=1}^n c_i^2$ 与 $sum_p (sum_{i=1}^n c_i)^2$ ，先考虑左式。
显然可以考虑 $(a_i,b_j)$ 点对贡献最后相加，那么答案可以写成
$$(n-1)! imes sum_{i=1}^nsum_{j=1}^ngcd(A_i,B_j)gcd(A_i,B_j)$$
对于右面的可以从大到小容斥得出有多少个点对的 $gcd=k$ 。
故处理当前部分的时间复杂度是调和级数 $mathcal O(mlog m)$ ，其中 $m$ 为数的值域。
同理考虑右式子，由于平方故要考虑两个点对的影响
$$(n-1)! imes sum_{i=1}^nsum_{j=1}^ngcd(A_i,B_j)gcd(A_i,B_j)\(n-2)!sum_{i=1}^nsum_{i'=1}^nsum_{j=1}^nsum_{j'=1}^n [i eq i'][j eq j']gcd(A_i,B_j)gcd(A_{i'},B_{j'})$$
为上面两式之和，第一行在上面已经算出故我们仅用考虑第二行。
由于 $ eq $ 关系非常难受，对 $ eq $ 进行容斥
$$sum_{i=1}^nsum_{i'=1}^nsum_{j=1}^nsum_{j'=1}^ngcd(A_i,B_j)gcd(A_{i'},B_{j'})\sum_{i=1}^nsum_{i'=1}^nsum_{j=1}^nsum_{j'=1}^n [i= i']gcd(A_i,B_j)gcd(A_{i'},B_{j'})\sum_{i=1}^nsum_{i'=1}^nsum_{j=1}^nsum_{j'=1}^n [j= j']gcd(A_i,B_j)gcd(A_{i'},B_{j'})\sum_{i=1}^nsum_{i'=1}^nsum_{j=1}^nsum_{j'=1}^n [i=i'][j= j']gcd(A_i,B_j)gcd(A_{i'},B_{j'})$$
其中第一个和第四个可以直接用算左式的结果，二三进行欧拉反演也可以做到 $mathcal O(mlog m)$ .
总时间复杂度 $mathcal O(mlog m)$ 。
#include<iostream> #include<cstring> #include<cstdio> #include<cstring> #include<vector> #include<queue> #include<algorithm> #include<climits> #define pii pair<int,int> #define pb push_back #define mp make_pair #define fi first #define se second #define int long long #define mod 1000000007 using namespace std; inline int read(){ int f=1,ans=0;char c=getchar(); while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();} while(c>='0'&&c<='9'){ans=ans*10+c-'0';c=getchar();} return f*ans; } const int MAXN=1e6+11; int v[MAXN],pri[MAXN],phi[MAXN],N,A[MAXN],B[MAXN],bucA[MAXN],bucB[MAXN],Lim=1000000,Ans1,Ans2; int BucA[MAXN],BucB[MAXN],f[MAXN],fac[MAXN]; void init(){ phi[1]=1; for(int i=2;i<=Lim;i++){ if(!v[i]){pri[++pri[0]]=i,v[i]=i,phi[i]=i-1;} for(int j=1;j<=pri[0]&&pri[j]*i<=Lim;j++){ v[i*pri[j]]=pri[j]; if(!(i%pri[j])){phi[i*pri[j]]=phi[i]*pri[j];break;} phi[i*pri[j]]=phi[i]*(pri[j]-1); } }return; } int gcd(int a,int b){return !b?a:gcd(b,a%b);} int ksm(int a,int b){int ans=1;while(b){if(b&1) ans*=a,ans%=mod;a*=a,a%=mod;b>>=1;}return ans;} int res1,res2,res3,res4,S1[MAXN],S2[MAXN]; signed main(){ //freopen("B.in","r",stdin); N=read(); for(int i=1;i<=N;i++) A[i]=read(),bucA[A[i]]++; for(int i=1;i<=N;i++) B[i]=read(),bucB[B[i]]++; init(); fac[0]=1; for(int i=1;i<=N;i++) fac[i]=fac[i-1]*i%mod; for(int d=1;d<=Lim;d++) for(int i=d;i<=Lim;i+=d) BucA[d]+=bucA[i],BucB[d]+=bucB[i]; for(int i=1;i<=Lim;i++) f[i]=BucA[i]*BucB[i]%mod; for(int i=Lim;i>=1;i--) for(int j=2*i;j<=Lim;j+=i) f[i]=(f[i]-f[j]+mod)%mod; for(int i=1;i<=Lim;i++) Ans1+=f[i]*i%mod*i%mod,Ans1%=mod; res4=Ans1; Ans1*=fac[N-1],Ans1%=mod; Ans2=Ans1; for(int i=1;i<=Lim;i++) res1+=f[i]*i%mod,res1%=mod; res1=res1*res1%mod; for(int i=1;i<=Lim;i++){ int w1=phi[i]*BucB[i]%mod,w2=phi[i]*BucA[i]%mod; for(int j=i;j<=Lim;j+=i) S1[j]+=w1,S1[j]%=mod,S2[j]+=w2,S2[j]%=mod; } for(int i=1;i<=N;i++) res2+=S1[A[i]]*S1[A[i]]%mod,res2%=mod; for(int i=1;i<=N;i++) res3+=S2[B[i]]*S2[B[i]]%mod,res3%=mod; Ans2+=(res1+res4-res2-res3)*fac[N-2]%mod,Ans2=(Ans2%mod+mod)%mod; //cerr<<Ans1<<" "<<Ans2<<endl; printf("%lld ",(Ans1*ksm(N,mod-2)%mod-Ans2*ksm(N*N%mod,mod-2)%mod+mod)%mod); }/* 3 1 2 3 1 2 3 */