[SDOI2018]旧试题
t组询问,询问\(\sum_{i=1}^A\sum_{j=1}^B\sum_{k=1}^Cd(ijk)mod\ 10^9+7\),\(1 ≤ T ≤ 10, 1 ≤ A, B, C ≤ 10^5, ,\)\(1≤∑max(A,B,C)≤2\times10 ^5\) 。
解
不难由约数拆分公式类比出
\[ans=\sum_{i=1}^A\sum_{j=1}^B\sum_{k=1}^C\sum_{x|i}\sum_{y|j}\sum_{z|k}\epsilon(gcd(x,y))\epsilon(gcd(y,z))\epsilon(gcd(x,z))=
\]
\[\sum_{x=1}^A\sum_{y=1}^B\sum_{z=1}^C\epsilon(gcd(x,y))\epsilon(gcd(y,z))\epsilon(gcd(x,z))[A/x][B/y][C/z]
\]
注意到多个元素的gcd,要保证gcd的相等,而倍数型更适用于只有gcd,约数条件的问题,于是考虑辅助公式,于是我们有
\[ans=\sum_{x=1}^A\sum_{y=1}^B\sum_{z=1}^C\sum_{i|x,i|y}\mu(i)\sum_{j|x,j|z}\mu(j)\sum_{k|y,k|z}\mu(k)[A/x][B/y][C/z]=
\]
\[\sum_{i=1}^{min(A,B)}\sum_{j=1}^{min(A,C)}\sum_{k=1}^{min(B,C)}\mu(i)\mu(j)\mu(k)\sum_{lcm(i,j)|x}\sum_{lcm(i,k)|y}\sum_{lcm(j,k)|z}[A/x][B/y][C/z]
\]
于是我们发现三组数,考虑三元环mobius,注意到一组i,j,k要有效,必然它们的\(\mu\)要有值,且小于A,B,C,而显然后式子可以\(O(nlog(n))\)维护,即两个数\(\mu\)有值,范围在A,B,C内,而lcm也没有超过A,B,C,我们就将其连边,这样我们就得到一张图,而解就是其中的三元环。
建图根据\(lcm(i,j)gcd(i,j)=ij\),枚举gcd,再枚举互质的两个数,边判\(\mu\)的是否有值以及是否超范围即可。
而三元环寻找显然是\(n^3\),考虑优化,这里采取度数大的点与度数小的点连边,这样我们的时间复杂度就只有\(O(m\sqrt{m})\)(n为点,m为边)
证明:
考虑度数大于\(\sqrt{m}\)的点,这样的点不能超过\(\sqrt{m}\)个,这里的枚举不超过\(m\sqrt{m}\),而小于的点,点不超过n所以枚举次数应为\(n\sqrt{m}\),因此当边数较多时可近似认为时间复杂度\(O(m\sqrt{m})\)。
于是我们有预处理后式的方法,又能较快建出了图,又有很快地找出三元环的方式,于是问题,注意判自环即可,至此得到很好解决。
参考代码:
#include <iostream>
#include <cstdio>
#include <cstring>
#define il inline
#define ri register
#define ll long long
#define plim 100000
#define yyb 1000000007
#define swap(x,y) x^=y^=x^=y
using namespace std;
struct point{
int next,to,len;
}ar[plim<<3];int at;
int head[plim+1];
struct edge{
int p1,p2,len;
}e[plim<<3];int et;
bool check[plim+1];
ll fa[plim+1],fb[plim+1],fc[plim+1];
int pt,prime[10000],mu[plim+1],ds[plim+1],
opt[plim+1],mark[plim+1];
il int max(int,int),gcd(int,int);
il void prepare(int),link(int,int,int);
int main(){
ll ans;
int lsy,A,B,C,mx,i,j,k;
scanf("%d",&lsy),prepare(plim);
while(lsy--){
scanf("%d%d%d",&A,&B,&C),mx=max(A,max(B,C)),ans&=0;
for(i=1;i<=A;++i)for(j=i;j<=A;j+=i)fa[i]+=A/j;
for(i=1;i<=B;++i)for(j=i;j<=B;j+=i)fb[i]+=B/j;
for(i=1;i<=C;++i)for(j=i;j<=C;j+=i)fc[i]+=C/j;
for(i=1;i<=mx;++i)ans+=mu[i]*fa[i]*fb[i]*fc[i];
for(i=1;i<=mx;++i){
if(!mu[i])continue;
for(j=1;j*i<=mx;++j){
if(!mu[j*i])continue;
for(k=j+1;(ll)j*i*k<=mx;++k){
if(!mu[i*k]||gcd(j,k)>1)continue;
int u(i*j),v(i*k),l(i*j*k);
ans+=mu[u]*(fa[v]*fb[l]*fc[l]+fa[l]*fb[v]*fc[l]+fa[l]*fb[l]*fc[v]);
ans+=mu[v]*(fa[u]*fb[l]*fc[l]+fa[l]*fb[u]*fc[l]+fa[l]*fb[l]*fc[u]);
e[++et].p1=u,e[et].p2=v,e[et].len=l,++ds[u],++ds[v];
}
}
}
for(i=1;i<=et;++i)
if(ds[e[i].p1]>ds[e[i].p2])link(e[i].p1,e[i].p2,e[i].len);
else link(e[i].p2,e[i].p1,e[i].len);
for(i=1;i<=mx;++i){
for(j=head[i];j;j=ar[j].next)mark[ar[j].to]=i,opt[ar[j].to]=ar[j].len;
for(j=head[i];j;j=ar[j].next)
for(k=head[ar[j].to];k;k=ar[k].next){
if(mark[ar[k].to]!=i)continue;
int d1(ar[j].len),d2(ar[k].len),d3(opt[ar[k].to]);
ans+=mu[i]*mu[ar[j].to]*mu[ar[k].to]*
(fa[d1]*fb[d2]*fc[d3]+fa[d1]*fb[d3]*fc[d2]+
fa[d2]*fb[d1]*fc[d3]+fa[d2]*fb[d3]*fc[d1]+
fa[d3]*fb[d1]*fc[d2]+fa[d3]*fb[d2]*fc[d1]);
}
}printf("%lld\n",ans%yyb);
memset(mark,0,sizeof(mark)),memset(fa,0,sizeof(fa)),memset(fb,0,sizeof(fb));
memset(fc,0,sizeof(fc)),memset(head,0,sizeof(head)),memset(ds,0,sizeof(ds)),
et&=at&=0;
}
return 0;
}
il void link(int x,int y,int len){
ar[++at].len=len,ar[at].to=y;
ar[at].next=head[x],head[x]=at;
}
il int gcd(int a,int b){
while(b)swap(a,b),b%=a;return a;
}
il int max(int a,int b){
return a>b?a:b;
}
il void prepare(int n){
int i,j;mu[1]|=true;
for(i=2;i<=n;++i){
if(!check[i])prime[++pt]=i,mu[i]=-1;
for(j=1;j<=pt&&prime[j]*i<=n;++j){
check[i*prime[j]]|=true;
if(!(i%prime[j]))break;
mu[i*prime[j]]=~mu[i]+1;
}
}
}