Luogu

LOJ

分析

首先要知道
$$
d(ijk)=\sum_{x|i}\sum_{y|j}\sum_{z|k}[\gcd(x,y)=1][\gcd(x,z)=1][\gcd(y,z)=1]
$$
现在可以来推式子了。（假设 $A<B<C$ ）
$$
\begin{aligned}
&\sum_{i=1}^A\sum_{j=1}^B\sum_{k=1}^Cd(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(x,z)=1][\gcd(y,z)=1]\\
=&\sum_{i=1}^A\sum_{j=1}^B\sum_{k=1}^C\Big\lfloor\frac{A}{i}\Big\rfloor\Big\lfloor\frac{B}{j}\Big\rfloor\Big\lfloor\frac{C}{k}\Big\rfloor[\gcd(i,j)=1][\gcd(i,k)=1][\gcd(j,k)=1]\\
=&\sum_{i=1}^A\sum_{j=1}^B\sum_{k=1}^C\Big\lfloor\frac{A}{i}\Big\rfloor\Big\lfloor\frac{B}{j}\Big\rfloor\Big\lfloor\frac{C}{k}\Big\rfloor\sum_{x|\gcd(i,j)}\mu(x)\sum_{y|\gcd(i,k)}\mu(y)\sum_{z|\gcd(j,k)}\mu(z)\\
=&\sum_{x=1}^{A}\mu(x)\sum_{y=1}^B\mu(y)\sum_{z=1}^C\mu(z)\sum_{\operatorname{lcm}(x,y)|i}^A\Big\lfloor\frac{A}{i}\Big\rfloor\sum_{\operatorname{lcm}(x,z)|j}^B\Big\lfloor\frac{B}{j}\Big\rfloor\sum_{\operatorname{lcm}(y,z)|k}^C\Big\lfloor\frac{C}{k}\Big\rfloor
\end{aligned}
$$
设
$$
f_y(x)=\sum_{x|i}^y\Big\lfloor\frac{y}{i}\Big\rfloor
$$
原式变为
$$
\sum_{x=1}^A\mu(x)\sum_{y=1}^B\mu(y)\sum_{z=1}^C\mu(z)f_a(\operatorname{lcm}(x,y))f_b(\operatorname{lcm}(x,z))f_c(\operatorname{lcm}(y,z))
$$
显然我们可以 $O(C\log C)$ 预处理出 $f_a,f_b,f_c$ ，然而计算还是 $O(n^3)$ 的。

但是注意到 $O(n^3)$ 的枚举中会出现很多贡献为 $0$ 的情况：

  • $\mu(x),\mu(y),\mu(z)$ 中有 $0$ 。
  • $\operatorname{lcm}(x,y)>A$ 或 $\operatorname{lcm}(x,z)>B$ 或 $\operatorname{lcm}(y,z)>C$ 。

那么有没有什么方法可以快速找出贡献不为 $0$ 的 $(x,y,z)$ 三元组呢？

可以枚举最大公约数 $d$ ，再枚举 $x=id,y=jd$ ，满足 $\gcd(i,j)=1,\mu(x)\neq 0,\mu(y)\neq 0,\operatorname{lcm}(u,v)=ijd\leq C$ ，然后连一条边 $(x,y)$ 。这样子所有的三元环 $(x,y,z)$ 就代表一个贡献不为 $0$ 的三元组。

感觉不太靠谱，但是复杂度 $O(m)$ 的，而且令人震惊的是只连了 $760741$ 条边。

于是我们可以直接 $O(m\sqrt{m})$ 的统计三元环并计算贡献了。

注意到这里的三元环是有序的，所以我们要把六种情况的贡献都加上。

以及上面的连边是没有算自环的，要把 $x,y,z$ 都相等和 $x,y,z$ 有两个相等的情况拿出来单独处理：第一种直接枚举是多少，第二种可以在连边的时候处理，将 $(x,x,y)$ 和 $(x,y,y)$ 的贡献都加上即可。

然后这题据说非常卡常，注意一下技巧：

  • 用 vector 存边，可以减少 cache miss 。
  • 最后答案不会超过 long long 范围，可以最后再取模。

上面算贡献的过程讲的可能不是很清楚，可以结合代码理解一下。

代码

大概还是能看的

// ===================================
//   author: M_sea
//   website: http://m-sea-blog.com/
// ===================================
#include <algorithm>
#include <iostream>
#include <cstdlib>
#include <cstring>
#include <cstdio>
#include <vector>
#include <cmath>
#define re register
using namespace std;

inline int read() {
    int X=0,w=1; char c=getchar();
    while (c<'0'||c>'9') { if (c=='-') w=-1; c=getchar(); }
    while (c>='0'&&c<='9') X=X*10+c-'0',c=getchar();
    return X*w;
}

const int N=100000+10,M=760741+10;
const int mod=1e9+7;

int p[N],v[N],cnt=0,mu[N];
inline void sieve(int n) {
    mu[1]=1;
    for (re int i=2;i<=n;++i) {
        if (!v[i]) p[++cnt]=i,mu[i]=-1;
        for (re int j=1;j<=cnt&&i*p[j]<=n;++j) {
            v[i*p[j]]=1;
            if (i%p[j]) mu[i*p[j]]=-mu[i];
            else break;
        }
    }
}

long long fa[N],fb[N],fc[N];
struct edge { int u,v,w; } e[M];
int deg[N],ecnt=0;
vector<pair<int,int> > G[N];
int vis[N];
inline int calc(int A,int B,int C) {
    if (A>C) swap(A,C); if (B>C) swap(B,C);
    for (re int i=1;i<=C;++i) deg[i]=fa[i]=fb[i]=fc[i]=0;
    for (re int i=1;i<=A;++i)
        for (re int j=i;j<=A;j+=i) fa[i]+=A/j;
    for (re int i=1;i<=B;++i)
        for (re int j=i;j<=B;j+=i) fb[i]+=B/j;
    for (re int i=1;i<=C;++i)
        for (re int j=i;j<=C;j+=i) fc[i]+=C/j;
    ecnt=0;
    for (re int i=1;i<=C;++i) G[i].clear();

    long long ans=0;
    // 三个相等
    for (re int i=1;i<=C;++i)
        if (mu[i]) ans+=mu[i]*mu[i]*mu[i]*fa[i]*fb[i]*fc[i];
    // 两个相等
    for (re int d=1;d<=C;++d)
        for (re int i=1;i*d<=C;++i) {
            if (!mu[i*d]) continue;
            for (re int j=i+1;1ll*i*j*d<=C;++j) {
                if (!mu[j*d]||__gcd(i,j)!=1) continue;
                int u=i*d,v=j*d,lcm=i*j*d;
                ++deg[u],++deg[v],e[++ecnt]=(edge){u,v,lcm};
                ans+=mu[u]*mu[u]*mu[v]*fa[u]*fb[lcm]*fc[lcm];
                ans+=mu[u]*mu[u]*mu[v]*fa[lcm]*fb[u]*fc[lcm];
                ans+=mu[u]*mu[u]*mu[v]*fa[lcm]*fb[lcm]*fc[u];
                ans+=mu[u]*mu[v]*mu[v]*fa[v]*fb[lcm]*fc[lcm];
                ans+=mu[u]*mu[v]*mu[v]*fa[lcm]*fb[v]*fc[lcm];
                ans+=mu[u]*mu[v]*mu[v]*fa[lcm]*fb[lcm]*fc[v];
            }
        }
    // 都不相等
    for (re int i=1;i<=ecnt;++i) {
        int u=e[i].u,v=e[i].v,w=e[i].w;
        if (deg[u]<deg[v]||(deg[u]==deg[v]&&u>v)) swap(u,v);
        G[u].push_back(make_pair(v,w));
    }
    for (re int u=1;u<=C;++u) { if (!mu[u]) continue;
        for (re auto i:G[u]) vis[i.first]=i.second;
        for (re auto i:G[u]) {
            int v=i.first; if (!mu[v]) continue;
            for (re auto j:G[v]) {
                int w=j.first; if (!mu[w]||!vis[w]) continue;
                int w_uv=i.second,w_uw=vis[w],w_vw=j.second;
                ans+=mu[u]*mu[v]*mu[w]*fa[w_uv]*fb[w_uw]*fc[w_vw];
                ans+=mu[u]*mu[v]*mu[w]*fa[w_uv]*fb[w_vw]*fc[w_uw];
                ans+=mu[u]*mu[v]*mu[w]*fa[w_uw]*fb[w_uv]*fc[w_vw];
                ans+=mu[u]*mu[v]*mu[w]*fa[w_uw]*fb[w_vw]*fc[w_uv];
                ans+=mu[u]*mu[v]*mu[w]*fa[w_vw]*fb[w_uv]*fc[w_uw];
                ans+=mu[u]*mu[v]*mu[w]*fa[w_vw]*fb[w_uw]*fc[w_uv];
            }
        }
        for (re auto i:G[u]) vis[i.first]=0;
    }
    return (ans%mod+mod)%mod;
}

int main() { sieve(1e5);
    for (re int T=read();T;--T) {
        int A=read(),B=read(),C=read();
        printf("%d\n",calc(A,B,C));
    }
    return 0;
}
最后修改：2021 年 03 月 23 日 10 : 13 PM
© 允许规范转载