U149858
作者：互联网
题意
给定两个长度为 $n$ 的序列 $a,b$。
对于一个 $1$ 到 $n$ 的排列 $p$，记 $c_i=\gcd(a_i,b_{p_i})$，$\sigma(c)$ 表示序列 $c$ 中所有元素的方差。
求 $$\sum\limits_{p}\sigma(c)$$ 对 $10^9+7$ 取模。
$2\leq n,a_i,b_i\leq 10^6$ 。
题解
方差可以写成以下形式
$$\begin{aligned}\sigma(c)&=\dfrac{1}{n}\sum_{i=1}^n (c_i-\bar{c})^2\\&=\dfrac{1}{n}\sum_{i=1}^n c_i^2-\dfrac{1}{n^2}(\sum_{i=1}^n c_i)^2\end{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)!\times \sum_{i=1}^n\sum_{j=1}^n\gcd(A_i,B_j)\gcd(A_i,B_j)$$
对于右面的可以从大到小容斥得出有多少个点对的 $\gcd=k$ 。
故处理当前部分的时间复杂度是调和级数 $\mathcal O(m\log m)$ ，其中 $m$ 为数的值域。
同理考虑右式子，由于平方故要考虑两个点对的影响
$$(n-1)!\times \sum_{i=1}^n\sum_{j=1}^n\gcd(A_i,B_j)\gcd(A_i,B_j)\\(n-2)!\sum_{i=1}^n\sum_{i'=1}^n\sum_{j=1}^n\sum_{j'=1}^n [i\neq i'][j\neq j']\gcd(A_i,B_j)\gcd(A_{i'},B_{j'})$$
为上面两式之和，第一行在上面已经算出故我们仅用考虑第二行。
由于 $\neq $ 关系非常难受，对 $\neq $ 进行容斥
$$\sum_{i=1}^n\sum_{i'=1}^n\sum_{j=1}^n\sum_{j'=1}^n\gcd(A_i,B_j)\gcd(A_{i'},B_{j'})\\\sum_{i=1}^n\sum_{i'=1}^n\sum_{j=1}^n\sum_{j'=1}^n [i= i']\gcd(A_i,B_j)\gcd(A_{i'},B_{j'})\\\sum_{i=1}^n\sum_{i'=1}^n\sum_{j=1}^n\sum_{j'=1}^n [j= j']\gcd(A_i,B_j)\gcd(A_{i'},B_{j'})\\\sum_{i=1}^n\sum_{i'=1}^n\sum_{j=1}^n\sum_{j'=1}^n [i=i'][j= j']\gcd(A_i,B_j)\gcd(A_{i'},B_{j'})$$
其中第一个和第四个可以直接用算左式的结果，二三进行欧拉反演也可以做到 $\mathcal O(m\log m)$ .
总时间复杂度 $\mathcal O(m\log 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\n",(Ans1*ksm(N,mod-2)%mod-Ans2*ksm(N*N%mod,mod-2)%mod+mod)%mod);
}/*
3
1 2 3
1 2 3
*/
View Code
标签：gcd,int,sum,U149858,neq,include,define 来源： https://www.cnblogs.com/si-rui-yang/p/14615246.html