P4619 [SDOI2018]旧试题

题目描述

\(T\) 组数据，计算

\[\sum_{i=1}^a\sum_{j=1}^b\sum_{k=1}^cd(ijk) \]

对 \(10^9+7\) 取模的结果。

数据范围

• \(1\le T\le 10,1\le a,b,c\le 10^5,\sum\max(a,b,c)\le 2\cdot 10^5\) 。

时间限制 \(\texttt{5s}\) ，空间限制 \(\texttt{512MB}\) 。

分析

题目名字中的 "旧试题" 指的是P3327 [SDOI2015] 约数个数和。

先考虑 \(d(ijk)\) 如何处理，熟知结论：

\[d(ijk)=\sum_{x=1}^i\sum_{y=1}^j\sum_{z=1}^k[\gcd(x,y)=1\and\gcd(y,z)=1\and\gcd(z,x)=1]\\ \]

证明很简单，每个素因子对左右两侧的贡献都是 \(\alpha+\beta+\gamma+1\) 。

带回原式，枚举 \(x,y,z\) ：

\[\sum_{x=1}^a\sum_{y=1}^b\sum_{z=1}^c[\gcd(x,y)=1\and\gcd(y,z)=1\and\gcd(z,x)=1]\lfloor\frac ax\rfloor\lfloor\frac by\rfloor\lfloor\frac cz\rfloor \]

莫比乌斯反演：

\[\sum_{x=1}^a\sum_{y=1}^b\sum_{z=1}^c\sum_{u|\gcd(x,y)}\mu(u)\sum_{v|\gcd(y,z)}\mu(v)\sum_{w|\gcd(z,x)}\mu(w)\lfloor\frac ax\rfloor\lfloor\frac by\rfloor\lfloor\frac cz\rfloor \]

枚举 \(u,v,w\) ：

\[\sum_{u=1}^{\min(a,b)}\mu(u)\sum_{v=1}^{\min(b,c)}\mu(v)\sum_{w=1}^{\min(c,a)}\mu(w)\sum_{\text{lcm}(u,w)|x}\lfloor\frac ax\rfloor\sum_{\text{lcm}(v,u)|y}\lfloor\frac by\rfloor\sum_{\text{lcm}(w,v)|z}\lfloor\frac cz\rfloor \]

记 \(f(n)=\sum_{i=1}^n\lfloor\frac ni\rfloor\) ，差分前缀和可以做到 \(\mathcal O(n\log n)\) 预处理。

于是原式变为：

\[\sum_{u=1}^{\min(a,b)}\sum_{v=1}^{\min(b,c)}\sum_{w=1}^{\min(c,a)}\mu(u)\mu(v)\mu(w)f(\lfloor\frac a{\text{lcm}(u,v)}\rfloor)f(\lfloor\frac b{\text{lcm(v,w)}}\rfloor)f(\lfloor\frac c{\text{lcm(w,u)}}\rfloor) \]

考虑哪些三元组 \((u,v,w)\) 会产生贡献。

• \(\mu(u)\neq 0,\mu(v)\neq 0,\mu(w)\neq 0\) 。
• \(\text{lcm}(u,v),\text{lcm}(v,w),\text{lcm}(w,u)\le\max(a,b,c)\) 。

将满足上述限制的 \((u,v)\) 连边，那么产生贡献的三元组就是三元环！

先考虑如何建图。枚举 \(d=\gcd(u,v)\) ，记 \(n=\max(a,b,c),u'=\frac ud,v'=\frac vd\) 。

限制为 \(\gcd(u',v')=1,d\cdot u'\cdot v'\le n\) ，可以在调和级数 \(\mathcal O(\frac nd\cdot\log\frac nd)\) 的时间内枚举出所有合法的 \((u',v')\) 。

因此完整的建图复杂度为 \(\sum_{d=1}^n\mathcal O(\frac nd\cdot\log\frac nd)=\mathcal O(n\log ^2n)\) 。

这张图其实非常稀疏，输出边数发现 \(m=760741\) ！

熟知结论三元环个数上限为 \(\mathcal O(m\sqrt m)\) ，直接枚举即可。

将计算 \(\text{lcm}\) 的代价视为 \(\mathcal O(1)\) ，则时间复杂度 \(\mathcal O(m\sqrt m)\) 。注意枚举全排列时需要精细实现，减少 \(\text{lcm}\) 的计算次数。

#include<bits/stdc++.h>
#define fi first
#define se second
#define mp make_pair
#define pii pair<int,int>
using namespace std;
const int maxn=1e5+5,mod=1e9+7;
int a,b,c,m,n,t,res;
pii e[800000];
int f[maxn],deg[maxn],vis[maxn];
vector<int> g[maxn];
#define b B
int b[maxn],p[maxn],mu[maxn];
inline void add(int &x,int y)
{
    if((x+=y)>=mod) x-=mod;
}
inline void dec(int &x,int y)
{
    if((x-=y)<0) x+=mod;
}
inline void init(int n)
{
    mu[1]=1;
    for(int i=2,cnt=0;i<=n;i++)
    {
        if(!b[i]) p[++cnt]=i,mu[i]=-1;
        for(int j=1;j<=cnt&&i*p[j]<=n;j++)
        {
            b[i*p[j]]=1;
            if(i%p[j]==0) break;
            mu[i*p[j]]=mu[i]*mu[p[j]];
        }
    }
    for(int i=1;i<=n;i++)
        for(int j=1;i*j<=n;j++)
            add(f[i*j],i),dec(f[min((i+1)*j,n+1)],i);
    for(int i=1;i<=n;i++) add(f[i],f[i-1]);
}
#undef b
inline int gcd(int a,int b)
{
    if(b==0) return a;
    return gcd(b,a%b);
}
inline int lcm(int a,int b)
{
    return a/gcd(a,b)*b;
}
inline int calc(int x,int y,int z)
{
    return 1ll*f[a/x]*f[b/y]%mod*f[c/z]%mod;
}
int main()
{
    scanf("%d",&t),init(maxn-5);
    while(t--)
    {
        scanf("%d%d%d",&a,&b,&c),m=0,n=max({a,b,c}),res=0;
        for(int i=1;i<=n;i++) deg[i]=vis[i]=0,g[i].clear();
        for(int d=1;d<=n;d++)
            for(int i=1;i<=n/d;i++)
                for(int j=i+1;1ll*i*j<=n/d;j++)
                    if(mu[i*d]&&mu[j*d]&&gcd(i,j)==1)
                        e[++m]=mp(i*d,j*d),deg[i*d]++,deg[j*d]++;
        for(int i=1;i<=m;i++)
        {
            int u=e[i].fi,v=e[i].se;
            if(mp(deg[u],u)<mp(deg[v],v)) g[u].push_back(v);
            else g[v].push_back(u);
        }
        for(int u=1;u<=n;u++)
        {
            if(!mu[u]) continue;
            res=(res+mu[u]*calc(u,u,u))%mod;
            for(auto v:g[u])
            {
                int x=lcm(u,v);
                vis[v]=u;
                res=(res+mu[v]*(0ll+calc(u,x,x)+calc(x,u,x)+calc(x,x,u)))%mod;
                res=(res+mu[u]*(0ll+calc(v,x,x)+calc(x,v,x)+calc(x,x,v)))%mod;
            }
            for(auto v:g[u])
                for(auto w:g[v])
                    if(vis[w]==u)
                    {
                        int x=lcm(u,v),y=lcm(v,w),z=lcm(w,u);
                        res=(res+mu[u]*mu[v]*mu[w]*(0ll+calc(x,y,z)+calc(x,z,y)
                        +calc(y,x,z)+calc(y,z,x)+calc(z,x,y)+calc(z,y,x)))%mod;
                    }
        }
        printf("%d\n",(res+mod)%mod);
    }
    return 0;
}