P4619 [SDOI2018]旧试题

P4619 [SDOI2018]旧试题

hoho！！！终于做掉了！！！Orz (color{black}{ exttt{c}}color{red}{ exttt{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}{xoperatorname{lcm}(a,c)}dfrac{B}{yoperatorname{lcm}(a,b)}dfrac{C}{zoperatorname{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}{xoperatorname{lcm}(a,c)}dfrac{B}{yoperatorname{lcm}(a,b)}dfrac{C}{zoperatorname{lcm}(b,c)}) 其实是可以预处理之后 (O(log)) 的。

[sum_{i=1}^{frac{A}{operatorname{lcm}(a,c)}}dfrac{A}{xoperatorname{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(msqrt m)) 级别的。这个可以参考网上一大堆三元环计数的博客，我们是通过枚举来统计答案的。

我们惊奇的发现可以枚举 ((a,b,c)) 这种三元环了，只要把上面统计个数直接改成拉边即可，而且数量不会很大。上面的暴力没有白写，随便改改就是正解。。。

对了有个显然的事情，邻接表拉边快遍历慢，vector 拉边慢遍历快，而且有O2加持。这题的瓶颈显然是遍历吧，拉边不过8e5次，应该用 vector 存图。不过邻接表会慢成什么样子没试过，小概率比 vector 快。

可是 (msqrt m) 上界不是 (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 了。

还有上面说可以预处理 (O(log)) 的那个式子，我们发现他可以不用取模，只不过极限情况有可能炸 long long （甚至可能根本炸不掉）。于是我们可以先拿一个大的数存下来，最后再取模，减少取模次数

发现两个优化很管用，本地直接跑到了 (5.2s)

勇了一发，最大点5.00s，AC了！！！喜提次劣解。woc怎么还有更慢的？？？

upd：更慢的那个不知道哪里去了，最劣解坐稳了。。。

注释代码来测瓶颈。发现瓶颈还是 (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
",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

*/