P4619 [SDOI2018]旧试题
作者：互联网
hoho！！！终于做掉了！！！Orz \(\color{black}{\texttt{c}}\color{red}{\texttt{yn2006}}\)
这个 \(d(ijk)\) 一看就是要拆的。不然那根本没法下手。
从 \(d(S)\) 的本质入手。\(S\) 的一部分来自 \(i\) ，一部分来自 \(j\) ，一部分来自 \(k\) 。
考虑直接直接 \(d(i)d(j)d(k)\) 算重的原因，发现是来自 \(i,j,k\) 的部分有些不互质，互质就能保证 \(ijk\) 的每一个因数只算一次。
所以我们有 \(d(ijk)=\sum_{x|i}\sum_{y|j}\sum_{z|k}[\gcd(x,y)=1][\gcd(y,z)=1][\gcd(z,x)=1]\)
带回去推式子
\[\sum_{i=1}^{A}\sum_{j=1}^{B}\sum_{k=1}^{C}d(ijk)\\ =\sum_{i=1}^{A}\sum_{j=1}^{B}\sum_{k=1}^{C}\sum_{x|i}\sum_{y|j}\sum_{z|k}[\gcd(x,y)=1][\gcd(y,z)=1][\gcd(z,x)=1]\\ =\sum_{x=1}^{A}\sum_{y=1}^{B}\sum_{z=1}^{C}[\gcd(x,y)=1][\gcd(y,z)=1][\gcd(z,x)=1]\dfrac{A}{x}\dfrac{B}{y}\dfrac{C}{z}\\ =\sum_{x=1}^{A}\sum_{y=1}^{B}\sum_{z=1}^{C}\sum_{a|\gcd(x,y)}\mu(a)\sum_{b|\gcd(y,z)}\mu(b)\sum_{c|\gcd(z,x)}\mu(c)\dfrac{A}{x}\dfrac{B}{y}\dfrac{C}{z}\\ =\sum_{a=1}^{A}\sum_{b=1}^{B}\sum_{c=1}^{C}\mu(a)\mu(b)\mu(c)\sum_{x=1}^{A}\sum_{y=1}^{B}\sum_{z=1}^{C}\sum_{a|x,a|y}\sum_{b|y,b|z}\sum_{c|x,c|z}\dfrac{A}{x}\dfrac{B}{y}\dfrac{C}{z}\\ =\sum_{a=1}^{A}\sum_{b=1}^{B}\sum_{c=1}^{C}\mu(a)\mu(b)\mu(c)\sum_{x=1}^{\frac{A}{\operatorname{lcm}(a,c)}}\sum_{y=1}^{\frac{B}{\operatorname{lcm}(a,b)}}\sum_{z=1}^{\frac{C}{\operatorname{lcm}(b,c)}}\dfrac{A}{x\operatorname{lcm}(a,c)}\dfrac{B}{y\operatorname{lcm}(a,b)}\dfrac{C}{z\operatorname{lcm}(b,c)}\\ \]最后面那个 \(\sum_{x=1}^{\frac{A}{\operatorname{lcm}(a,c)}}\sum_{y=1}^{\frac{B}{\operatorname{lcm}(a,b)}}\sum_{z=1}^{\frac{C}{\operatorname{lcm}(b,c)}}\dfrac{A}{x\operatorname{lcm}(a,c)}\dfrac{B}{y\operatorname{lcm}(a,b)}\dfrac{C}{z\operatorname{lcm}(b,c)}\) 其实是可以预处理之后 \(O(\log)\) 的。
\[\sum_{i=1}^{\frac{A}{\operatorname{lcm}(a,c)}}\dfrac{A}{x\operatorname{lcm}(a,c)}=\sum_{i=1}^{\frac{A}{\operatorname{lcm}(a,c)}}\dfrac{\dfrac{A}{\operatorname{lcm}(a,c)}}{x} \]我们预处理 \(f_x=\sum_{i=1}^{x}\dfrac{x}{i}\) 就可以 \(O(\log)\) 搞上面那个东西。这个东西水爆了，不会建议回炉重造。
但是千万注意这里有个 \(\log\) ，不要像我一样傻乎乎的忽略了 \(\operatorname{lcm}\) 的复杂度当成 \(O(1)\) ，不然后面会死的很惨。
发现完全没法往下做了。。。
这时候需要一点数学（数论）的直觉。
首先 \(\mu(a)\) 不能为 \(0\) ，其次 \(\operatorname{lcm}(a,b) \le \max(A,B,C)\) ，感觉这个东西特别少的样子。
写个暴力输出一下 \(10^5\) 以内 \((a,b)\) 的个数，你发现纯暴力直接T没。
于是我们开始优化这个暴力求个数的过程。
首先套路的枚举一下 \(\gcd\) 。我们在对于每一个数 \(x\) 开vector存它倍数，如果 \(\mu(x)\) 不为 \(0\) 那么我们平方遍历它所有的倍数， \(\gcd(d_j,d_k)=x\) 的时候计数器++。
你发现这个东西跑的飞快，而且，计数器只有 821535 ！！！
有一个特别牛逼的结论，那就是三元环的数量上限是 \(O(m\sqrt m)\) 级别的。这个可以参考网上一大堆三元环计数的博客，我们是通过枚举来统计答案的。
我们惊奇的发现可以枚举 \((a,b,c)\) 这种三元环了，只要把上面统计个数直接改成拉边即可，而且数量不会很大。
可是上界不是 \(10^9\) 吗？
别忘了自然数是很神奇的，这又不是出题人手造恶意卡的图。。。
事实上上限大概 \(10^7\) 。
吼啊，直接枚举！就完结了？？？
极限数据应该是
2
100000 100000 100000
100000 100000 100000
因为这个时候那个 \(\sqrt m\) 最大。
之前提过了 \(\operatorname{lcm}\) 是有 \(\log\) 的。。。于是直接T飞，但是只用8s，没有特别离谱
首先一波乱搞发现 \(\operatorname{lcm}\) 不用开 long long ，就是把那些 \(\operatorname{lcm}\) 大于 \(10^5\) 的直接舍掉，然后 先除以 \(\gcd\) 再乘，即 lcm(x,y)=x/gcd(x,y)*y ，这样就不用 long long 了。
你发现这个优化很管用，本地直接跑到了 \(5.2s\)
勇了一发，最大点5.00s，AC了！！！喜提次劣解。woc怎么还有更慢的？？？
发现瓶颈还是 \(\operatorname{lcm}\) 。我们发现这些边的 \(\operatorname{lcm}\) 是可以直接预处理的。于是大大减少了 \(\operatorname{lcm}\) 的次数。
于是成功跑进了 \(3s\) ，稳了。
//Orz cyn2006
#include<bits/stdc++.h>
using namespace std;
#define fi first
#define se second
#define mkp(x,y) make_pair(x,y)
#define pb(x) push_back(x)
#define sz(v) (int)v.size()
typedef long long LL;
typedef double db;
template<class T>bool ckmax(T&x,T y){return x<y?x=y,1:0;}
template<class T>bool ckmin(T&x,T y){return x>y?x=y,1:0;}
#define rep(i,x,y) for(int i=x,i##end=y;i<=i##end;++i)
#define per(i,x,y) for(int i=x,i##end=y;i>=i##end;--i)
inline int read(){
int x=0,f=1;char ch=getchar();
while(!isdigit(ch)){if(ch=='-')f=0;ch=getchar();}
while(isdigit(ch))x=x*10+ch-'0',ch=getchar();
return f?x:-x;
}
const int N=100005;
const int M=821535+5;
#define mod 1000000007
void fmod(int&x){x+=x>>31&mod,x-=mod,x+=x>>31&mod;}
int A,B,C,mx,ans;
int mu[N],pri[N],pct,sum[N],tot,tag[N];
bool vis[N];
vector<int>d[N];
int U[M],V[M],ord[N],deg[N],tu[M],tv[M],cnt,LCM[M],id[M];
vector<pair<int,int> >e[N];
pair<int,int>p[N];
void init(const int&N){
mu[1]=1;
for(int i=2;i<=N;++i){
if(!vis[i])pri[++pct]=i,mu[i]=-1;
for(int j=1;j<=pct&&i*pri[j]<=N;++j){
vis[i*pri[j]]=1;
if(i%pri[j]==0){mu[i*pri[j]]=0;break;}
mu[i*pri[j]]=-mu[i];
}
}
for(int i=1;i<=N;++i){
if(!mu[i])continue;
for(int j=1;i*j<=N;++j)if(mu[i*j])d[i].pb(i*j);
}
for(int i=1;i<=N;++i)
for(int j=0,up=sz(d[i]);j<up;++j)
for(int k=j;k<up&&1ll*d[i][j]*d[i][k]/i<=N;++k)
if(__gcd(d[i][j],d[i][k])==i)U[++tot]=d[i][j],V[tot]=d[i][k],LCM[tot]=d[i][j]/i*d[i][k];
for(int i=1;i<=N;++i){
int res=0;
for(int l=1,r;l<=i;l=r+1)
r=i/(i/l),res+=(i/l)*(r-l+1);
sum[i]=res;
}
}
int lcm(int x,int y){return x/__gcd(x,y)*y;}
LL f(int x,int y,int z){return 1ll*sum[x]*sum[y]*sum[z];}
void Main(){
A=read(),B=read(),C=read(),mx=max(A,max(B,C)),fill(deg+1,deg+mx+1,0),fill(tag+1,tag+mx+1,0),ans=0,cnt=0;
rep(i,1,tot)if(LCM[i]<=mx)++deg[U[i]],++deg[V[i]],++cnt,tu[cnt]=U[i],tv[cnt]=V[i],id[cnt]=i;
rep(i,1,mx)p[i]=mkp(deg[i],i);
sort(p+1,p+mx+1);
rep(i,1,mx)ord[p[i].se]=i,e[i].clear();
rep(i,1,cnt)
if(ord[tu[i]]>ord[tv[i]])e[tu[i]].pb(mkp(tv[i],LCM[id[i]]));
else e[tv[i]].pb(mkp(tu[i],LCM[id[i]]));
rep(u,1,mx){
for(int i=0,up=sz(e[u]);i<up;++i)tag[e[u][i].fi]=u;
for(int i=0,up1=sz(e[u]);i<up1;++i){
int v=e[u][i].fi;
for(int j=0,up2=sz(e[v]);j<up2;++j){
int w=e[v][j].fi;
if(tag[w]==u){
int x=e[u][i].se,y=e[v][j].se,z=lcm(w,u);unsigned long long res=0;
if(u!=v&&u!=w&&v!=w)
res=f(A/x,B/y,C/z)+f(A/x,B/z,C/y)+f(A/y,B/x,C/z)+f(A/y,B/z,C/x)+f(A/z,B/x,C/y)+f(A/z,B/y,C/x);
else if(u==v&&u==w&&v==w)res=f(A/x,B/y,C/z);
else res=f(A/x,B/y,C/z)+f(A/y,B/z,C/x)+f(A/z,B/x,C/y);
int t=res%mod;
fmod(ans+=1ll*mu[u]*mu[v]*mu[w]*t);
}
}
}
}
printf("%d\n",ans);
}
signed main(){
init(100000);
for(int T=read();T;--T)Main();
}
/*
5
10 10 10
100 100 100
1000 1000 1000
10000 10000 10000
100000 100000 100000
2
100000 100000 100000
100000 100000 100000
*/
标签：P4619,试题,dfrac,sum,mu,SDOI2018,lcm,operatorname,gcd 来源： https://www.cnblogs.com/zzctommy/p/14050558.html