SDOI2018-旧试题

题意

题意：给定\(A,B,C\)，求

\[\sum_{i=1}^A \sum_{j=1}^B \sum_{k=1}^C d(i\times j\times k) \]

其中\(d(n)\)表示\(n\)的约数个数，即\(d(n)=\sum_{k|n}1\)，\(1\leq \max(A,B,C)\leq 10^5\)，多组询问，\(\sum \max(A,B,C)\leq 2\times 10^5\)，时限5s.

题解

题解：

首先这题是有一些莫反以外的做法的，可以直接考虑容斥，但我还没看懂，比如博客： https://www.cnblogs.com/abigjiong/articles/17446718.html 就有说，但这题即使单看莫反+三元环计数的做法也算是比较巧妙了，还是值得一说。

前置知识

有一个经典的问题是（弱化版本）：\(\sum_{i=1}^n \sum_{j=1}^m d(i\times j)\) 如何计算？如果去考虑 \(i\times j\) 的唯一分解，会得到

\[d(i\times j)=\sum_{i'|i}\sum_{j'|j}[\gcd(i',j')=1] \]

证明起来大概是要考虑一种一一对应关系，这里就略去了。

类似的可以得到多个乘积的情形：

\[d(\prod_{i=1}^n x_i)=\sum_{\forall i:y_i|x_i}[y_i两两互质] \]

反演

所以到这个题里来，就有：

\[\begin{aligned} ans&=\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 \lfloor\frac{A}{x}\rfloor \lfloor\frac{B}{y}\rfloor \lfloor\frac{C}{z}\rfloor [\gcd(x,y)=1]\times [\gcd(y,z)=1]\times [\gcd(z,x)=1]\\ &=\sum_{x=1}^A\sum_{y=1}^B\sum_{z=1}^C \lfloor\frac{A}{x}\rfloor \lfloor\frac{B}{y}\rfloor \lfloor\frac{C}{z}\rfloor \sum_{p|x,p|y}\mu(p)\times \sum_{q|y,q|z}\mu(q) \times \sum_{r|z,r|x}\mu(r)\\ &=\sum_{p} \sum_q \sum_r \mu(p)\mu(q)\mu(r) \sum_{x=1}^A \lfloor\frac{A}{lcm(p,r)\times x}\rfloor \sum_{y=1}^B \lfloor\frac{B}{lcm(p,q)\times y}\rfloor \sum_{z=1}^C \lfloor\frac{C}{lcm(q,r)\times z}\rfloor \end{aligned} \]

注意后面的式子的形式，形如：

\[F(A,L)=\sum_{x=1}^{A/L} \lfloor \frac{A}{xL}\rfloor \]

至多只有\(\frac{A}{L}\)项，\(A,B,C\) 至多三种取值，对每个\(L\)是可以\(O(A\log A)\)计算的。

所以答案形如：

\[\begin{aligned} ans&=\sum_{p} \sum_q \sum_r \mu(p)\mu(q)\mu(r) \times F(A,lcm(p,r))\times F(B,lcm(p,q))\times F(C,lcm(q,r)) \end{aligned} \]

很好…这时候一个只学过数论的选手停止了思考…

上界估计

注意到\(F(A,L)\) 非零，当且仅当\(L\leq A\) 的，放在这里的情况，就是需要形如 \(lcm(p,r)\leq A\) 的条件，这是不是还挺苛刻的？

估计一个 \(\sum_{i=1}^n \sum_{j=1}^n [lcm(i,j)\leq n]\)的上界：

\[\begin{aligned} S&=\sum_{i=1}^n \sum_{j=1}^n [ij\leq n\times \gcd(i,j)]=\sum_{i=1}^n \sum_{j=1}^n \sum_{d=1}^n[\gcd(i,j)=d][i\times j\leq n\times d]\\ &=\sum_{k=1}^n \mu(k)\sum_{d=1}^n \ \sum_{i=1}^n \sum_{j=1}^n \ [kd|i]\times [kd|j]\times[ij\leq nd]\\ &\leq \sum_{k=1}^n \sum_{d=1}^n \sum_{i=1}^{n/kd}\sum_{j=1}^{n/kd} [(kdi)\times (kdj)\leq nd] \end{aligned} \]

中间的条件是\(k^2 d^2 \times i\times j\leq nd\)，即\(j\leq \frac{n}{k^2 d \times i}\)，再求和：

\[\sum_{i=1}^{n/kd}\frac{n}{k^2 d\times i}=O(\frac{n}{k^2 d}\log \frac{n}{kd}) \]

所以

\[S= O(\sum_{k=1}^n \sum_{d=1}^n \frac{n}{k^2 d}\log n)=O(n\log n\sum_{k=1}^n \frac{1}{k^2 }\sum_{d=1}^n \frac{1}{d})=O(n\log^2 n) \]

中间有一步 \(\sum_{k=1}^n\frac{1}{k^2}=O(1)\)，如果对高数还有印象的话可能会记得\(\sum_{k=1}^\infty \frac{1}{k^2}=\pi^2 /6\)，所以有限项求和自然是\(O(1)\)的。

到此为止我们会知道，符合\(lcm(i,j)\leq n\)且\(1\leq i,j\leq n\)的序对\((i,j)\)大约是\(O(n\log ^2 n)\)级别的。

注：事实上我们上面估计的时候用到了一个放缩：\(\sum_{k=1}^n \mu(k)/k^2\leq \sum_{k=1}^n 1/k^2\)，实际上这个放缩是非常“松”的，如果对Dirichlet生成函数熟悉的话，设 \(M(z)=\sum_{n=1}^\infty \frac{\mu(z)}{n^z}\)，右边一项趋于无穷时恰是\(M(2)\)，而我们有\(M(z)\times \zeta(z)=1\)，所以\(M(2)=\frac{1}{\zeta(2)}=6/\pi^2\)，大约是 \(0.61\) 的一个常数。

三元环

但是如果把每个数\(p,q,r\)看成一张图上的点，两个点 \(p,q\) 有边当且仅当\(lcm(p,q)\leq n\)，那么这张图的边数是\(O(n\log ^2 n)\)的，我们完全可以暴力建图出来！

（实际实现的时候应该枚举LCM而不是枚举 \(p,q\) ）

原问题变成枚举图上的每个三元环，两个点\(p,q\)的边权就是 \(F(\dots,lcm(p,q))\) 这种东西（对于固定的边，自然可以\(O(1)\) 知道）。

而\(m\)个点的无向图的三元环计数是可以\(O(m\sqrt m)\)做到的！所以这道题的时间上界就是\(O(n\sqrt n\log ^3 n)\)的！！！

而且其实一方面边数卡不满\(O(n\log ^2 n)\)这个上界（常数只有\(0.6\)左右，而且我们只关心双向边中的一条边，实际实现的时候大约只有几十万条边），另一方面是三元环计数的 \(O(m\sqrt m)\) 如果要卡满也是依赖于图的性质，而我们这张图是基于数论性质建出来的，loj上LCA有给出这个问题一个更紧的上界，会发现远远卡不满。

完结撒花

到此为止我就基本完成了这道神奇的题目…顺带分析了某个式子的上界，看了眼Dirichlet生成函数，学了一下三元环计数的trick（说实话在做这题之前还没遇到过三元环计数），当然最关键的还是这个问题转化的思路啦…非常奇妙

存个代码：

#include<bits/stdc++.h>
#define rep(i,a,b) for(int i=(a);i<=(b);i++)
#define fastio ios_base::sync_with_stdio(false);cin.tie(0);cout.tie(0)
#define mp make_pair
#define pb push_back
using namespace std;
typedef long long ll;
const int N=1e5+5;
const int MOD=1e9+7;
const int M=1e6;
int cnt,pri[N],mu[N];
bool not_pri[N];
ll lcm(int x,int y){return (ll)x*y/__gcd(x,y);}
void init(){
	mu[1]=1;
	for(int i=2;i<=N-5;i++){
		if(!not_pri[i]){
			pri[++cnt]=i;
			mu[i]=-1;
		}
		for(int j=1;j<=cnt&&i*pri[j]<=N-5;j++){
			not_pri[i*pri[j]]=1;
			if(i%pri[j]==0){
				mu[i*pri[j]]=0;
				break;
			}else mu[i*pri[j]]=-mu[i];
		}
	}
}
int MX,a[3];
ll f[N][3],ans;
int u[M],v[M],w[M],cnt_edge;
vector<vector<pair<int,int>>> G;
vector<int> deg;
void prework(){
	rep(i,0,2)cin>>a[i];
	sort(a,a+3);
	MX=a[2];
	G=vector<vector<pair<int,int>>>(MX+1);
	deg=vector<int>(MX+1);

	rep(i,1,a[2])rep(x,0,2)f[i][x]=0;
	rep(x,0,2)rep(i,1,a[x])for(int k=a[x]/i;k>=1;k--)f[k][x]+=a[x]/(k*i);

	cnt_edge=0;
	auto addEdge=[](int i,int j,int W){
		if(i==j)return;
		cnt_edge++;
		u[cnt_edge]=i;
		v[cnt_edge]=j;
		w[cnt_edge]=W;
	};
	//i*j/gcd()
	rep(i,1,MX)if(mu[i])for(int L=i;L<=MX;L+=i)for(int j=1;j*j<=L;j++)if(L%j==0){
		int t=L/j;
		if(j<=i&&mu[j]&&lcm(i,j)==L)addEdge(i,j,L);
		if(t<=i&&mu[t]&&t!=j&&lcm(i,t)==L)addEdge(i,t,L);
	}
	rep(i,1,cnt_edge){
		if(deg[u[i]]<deg[v[i]])swap(u[i],v[i]);
		if(deg[u[i]]==deg[v[i]]){
			if(u[i]>v[i])G[u[i]].pb(mp(v[i],w[i]));
			else G[v[i]].pb(mp(u[i],w[i]));
		}else G[u[i]].pb(mp(v[i],w[i]));
	}
}
void add(ll &x,ll y){x+=y;if(x>=MOD)x-=MOD;}
void upd(int sign,ll v){
	if(sign>0)add(ans,v);
	if(sign<0)add(ans,MOD-v);
}
void work(){
	ans=0;
	rep(i,1,a[0])upd(mu[i],f[i][0]*f[i][1]*f[i][2]%MOD);
	rep(i,1,MX)for(auto [j,v1]:G[i]){
		//case: i,j,k different each other
		for(auto [k,v2]:G[j]){
			ll v0=lcm(i,k);
			if(v0>MX)continue;
			int s=mu[i]*mu[j]*mu[k];
			upd(s,f[v0][0]*(f[v1][1]*f[v2][2]+f[v2][1]*f[v1][2])%MOD);
			upd(s,f[v1][0]*(f[v0][1]*f[v2][2]+f[v2][1]*f[v0][2])%MOD);
			upd(s,f[v2][0]*(f[v0][1]*f[v1][2]+f[v1][1]*f[v0][2])%MOD);
		}
		//存在两个数相同
		//i-i-j,i-j-i,j-i-i
		ll sq;

		sq=f[v1][1]*f[v1][2];
		upd(mu[j],f[i][0]*sq%MOD);
		upd(mu[i],f[j][0]*sq%MOD);

		sq=f[v1][0]*f[v1][2];
		upd(mu[j],f[i][1]*sq%MOD);
		upd(mu[i],f[j][1]*sq%MOD);

		sq=f[v1][0]*f[v1][1];
		upd(mu[j],f[i][2]*sq%MOD);
		upd(mu[i],f[j][2]*sq%MOD);
	}
	cout<<ans%MOD<<'\n';
}
int main(){
	fastio;
	init();
	int tc;cin>>tc;
	while(tc--){
		prework();
		work();
	}
	return 0;
}