再做一遍,算是复健一下数论。
Description
\(T\) 组数据,给出 \(A,B,C\),求
\[\sum_{i=1}^A\sum_{j=1}^B\sum_{k=1}^Cd(ijk) \]\(d\) 表示正因子个数。答案对 \(P=10^9+7\) 取模。
\(1\le T\le 10,1\le A,B,C\le 10^5,1\le \sum\max\{A,B,C\}\le 2\times10^5.\)
Solution
只能说是神仙题。
以下内容会进行一个正解的棒读,但并不保证都是人能想出来的东西。
引理:
\[d(ijk)=\sum_{x\mid i}\sum_{y\mid j}\sum_{z\mid k}[x\perp y][y\perp z][z\perp x] \]证明:
记 \(\nu_p(a)\) 表示 \(a\) 包含因子 \(p\) 的个数。
注意到每个质因子 \(p\) 对两边的式子的贡献都是独立的。从而可以得到
\[LHS=\prod(\nu_p(ijk)+1)=\prod (\nu_p(i)+\nu_p(j)+\nu_p(k)+1)=RHS \]证毕。
其实这段证明按 MO 的标准来看是非常平凡的。
但问题是考场上可能根本想不到这个引理。
不过日常做题还是有机会遇到这个的:
\[d(ij)=\sum_{x\mid i}\sum_{y\mid j}[x\perp y] \]如果以前见过的话,自己推出引理的可能性微存。(大嘘
总之,有了引理之后,就可以直接代换掉原式中的 \(d(ijk)\) 了。
\[ans=\sum_{i=1}^A\sum_{j=1}^B\sum_{k=1}^C\sum_{x\mid i}\sum_{y\mid j}\sum_{z\mid k}[x\perp y][y\perp z][z\perp x] \]交换求和号。
\[\def\div#1#2{\left\lfloor\dfrac{#1}{#2}\right\rfloor} ans=\sum_{x=1}^A\sum_{y=1}^B\sum_{z=1}^C[x\perp y][y\perp z][z\perp x]\div Ax\div By\div Cz\]注意到 \(x\perp y,y\perp z,z\perp x\) 这三个约束条件可以构成一个三元环。
众所周知,一张图上三元环的个数最多为 \(O(m\sqrt m)\),其中 \(m\) 为边数。这也是三元环计数的时间复杂度。
但是我们写个代码算一下边数,发现 \(m\) 可以达到 3e9,所以这个做法寄了。
寄掉的做法为什么还要提呢?因为待会还有用。
调整一下思路,采取数论题的常规做法:大力反演。
\[\def\div#1#2{\left\lfloor\dfrac{#1}{#2}\right\rfloor} ans=\sum_{x=1}^A\sum_{y=1}^B\sum_{z=1}^C\sum_{u|x,u|y}\mu(u)\sum_{v|y,v|z}\mu(v)\sum_{w|z,w|x}\mu(w)\div Ax\div By\div Cz\]交换求和号。
\[\def\div#1#2{\left\lfloor\dfrac{#1}{#2}\right\rfloor} ans=\sum_{u=1}^A\sum_{v=1}^B\sum_{w=1}^C\mu(u)\mu(v)\mu(w) \sum_{{\rm lcm}(w,u)|x}\div Ax\sum_{{\rm lcm}(u,v)|y}\div By\sum_{{\rm lcm}(v,w)|z}\div Cz\]记 \(\displaystyle f(d,A)=\sum_{d\mid x}\left\lfloor\dfrac Ax\right\rfloor\),则有
\[ans=\sum_{u=1}^A\sum_{v=1}^B\sum_{w=1}^C\mu(u)\mu(v)\mu(w) \cdot f({\rm lcm}(w,u),A)\cdot f({\rm lcm}(u,v),B)\cdot f({\rm lcm}(v,w),C)\]到目前为止复杂度仍然是 \(O(n^3)\) 的,而且感觉已经推不下去了,非常寄。
但推到这一步之后又可以用三元环计数做了。
我们将 \(u,v,w\) 看作点,两点 \(u,v\) 之间连一条权值为 \({\rm lcm}(u,v)\) 的边。
不难想到如下优化:
-
如果 \(u\) 满足 \(\mu(u)=0\),则可以直接删掉这个点;
-
如果 \(u,v\) 满足 \({\rm lcm}(u,v)>\max\{A,B,C\}\),那么可以删掉 \(u\) 和 \(v\) 之间的边。
写个程序验证一下,发现极限情况下的边数为 \(760741\)。
而三元环计数的时间复杂度为 \(O(m\sqrt m)\),已经足以通过。
当然,常数太大也不行(
一些细节问题:
建图时,建议用 dfs 依次确定 \(u,v\) 的每个质因子,并且用 \(\rm lcm\) 进行剪枝。
这样的做法应该是近似 \(O(m)\) 的,会比 \(O(n^2)\) 暴力建图快不少。
建议用 vector 存边而不是链式前向星。这样可以使内存访问更加连续,常数更优。
记得要去掉重边。
以及,本题中三元环不一定是三个不同的点,注意考虑特殊情况。
View Code
#include<stdio.h>
#include<vector>
#define repE(i,x) for (auto i:E[x]) if (i.y<=n)
typedef long long ll;
const int N=1e5+5,M=8e5+5,P=1e9+7;
const double eps=1e-14;
int n,m,cnt,T,A,B,C,ans;
int p[20005],mu[N],fA[N],fB[N],fC[N],d[N];
bool vp[N]; double I[N];
struct Node { int x,y,lcm; }a[M];
struct Node2 {int y,lcm; }v[N];
std::vector<Node2>E[N];
inline int max(int x,int y) { return x>y?x:y; }
void dfs(int x,int y,int i,int lcm) {
if (x<y)
++d[x],++d[y],
a[++m]={x,y,lcm};
for (int j=i+1; j<=cnt; ++j) {
if (1ll*lcm*p[j]>n) return ;
int tl=lcm*p[j];
dfs(x*p[j],y,j,tl);
dfs(x,y*p[j],j,tl);
dfs(x*p[j],y*p[j],j,tl);
}
}
inline void inc1(int x) { //1个点的三元环
ans=(1ll*mu[x]*fA[x]*fB[x]*fC[x]%P+ans+P)%P;
}
inline void inc2(int x,int y,int l) { //2个点的三元环
ll s1=1ll*fA[x]*fB[l]*fC[l]
+1ll*fA[l]*fB[x]*fC[l]
+1ll*fA[l]*fB[l]*fC[x];
ll s2=1ll*fA[y]*fB[l]*fC[l]
+1ll*fA[l]*fB[y]*fC[l]
+1ll*fA[l]*fB[l]*fC[y];
ans=((s1*mu[y]+s2*mu[x])%P+ans+P)%P;
}
inline void inc3(int x,int y,int z,int l1,int l2,int l3) { //3个点的三元环
ll s1=(1ll*fA[l1]*fB[l2]+1ll*fA[l2]*fB[l1])*fC[l3];
ll s2=(1ll*fA[l2]*fB[l3]+1ll*fA[l3]*fB[l2])*fC[l1];
ll s3=(1ll*fA[l3]*fB[l1]+1ll*fA[l1]*fB[l3])*fC[l2];
ans=((s1+s2+s3)*mu[x]*mu[y]*mu[z]%P+ans+P)%P;
}
int main() {
n=100000;
mu[1]=1,I[1]=1;
for (int i=2; i<=n; ++i) { //筛mu
if (!vp[i]) mu[i]=-1,p[++cnt]=i;
for (int j=1,t; j<=cnt&&(t=p[j]*i)<=n; ++j) {
vp[t]=1;
if (i%p[j]==0) break;
mu[t]=-mu[i];
}
I[i]=1./i+eps;
}
scanf("%d",&T);
while (T--) {
scanf("%d%d%d",&A,&B,&C);
n=max(max(A,B),C),m=0,ans=0;
dfs(1,1,0,1); //建边
for (int i=1; i<=m; ++i) { //建图
int x=a[i].x,y=a[i].y,lcm=a[i].lcm;
if (d[x]<=d[y]) E[x].push_back({y,lcm});
else E[y].push_back({x,lcm});
}
for (int i=1; i<=n; ++i) //预处理fA,fB,fC
for (int j=i; j<=n; j+=i)
fA[i]+=A*I[j],
fB[i]+=B*I[j],
fC[i]+=C*I[j];
for (int i=1; i<=m; ++i) //三元环为(x,x,y)或(x,y,y)
inc2(a[i].x,a[i].y,a[i].lcm);
for (int i=1; i<=n; ++i) if (mu[i]) {
inc1(i); //三元环为(i,i,i)
repE(j,i)
v[j.y]={i,j.lcm};
repE(j,i) repE(k,j.y) if (v[k.y].y==i)
inc3(i,j.y,k.y,j.lcm,k.lcm,v[k.y].lcm); //三元环为(i,j.y,k.y)
}
printf("%d\n",ans);
for (int i=1; i<=n; ++i) //初始化
d[i]=0,E[i].clear(),
fA[i]=fB[i]=fC[i]=0;
}
return 0;
}