U149858

题意

给定两个长度为 $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
  */