[SDOI2018] 旧试题

没见过数论题还卡常的，幸好我常数小。

一、题目

点此看题

二、解法

你看了题目名和给出的柿子，心想这题不是考过的吗？[SDOI2015]约数个数和，不得不说这道题还是很有创意的，如果你做过以前那道题就会发现这个 ( t tirck) 还是可以用：

[d(ijk)=sum_{x|i}sum_{y|j}sum_{z|k}epsilon((x,y))epsilon((x,z))epsilon((y,z)) ]

(epsilon) 表示单位元函数，((x,y)) 表示他们的最大公约数，然后整个表达式可以改写成这样：

[sum_{i=1}^Asum_{j=1}^Bsum_{k=1}^Csum_{x|i}sum_{y|j}sum_{z|k}epsilon((x,y))epsilon((x,z))epsilon((y,z)) ]

然后把 (x,y,z) 放到最前面去枚举：

[sum_{x=1}^Asum_{y=1}^Bsum_{z=1}^Cepsilon((x,y))epsilon((x,z))epsilon((y,z))frac{A}{x}frac{B}{y}frac{C}{z} ]

这时候莫比乌斯反演闪亮登场，把 ( t gcd) 给反演了：

[sum_{x=1}^Asum_{y=1}^Bsum_{z=1}^C(sum_{u|(x,y)}mu(u))(sum_{v|(x,z)}mu(v))(sum_{w|(y,z)}mu(w))frac{A}{x}frac{B}{y}frac{C}{z} ]

然后 (u,v,w) 放在最前面去枚举，这时候 (u,v) 都是 (x) 的因数，所以 (x) 要是 ( t lcm(u,v)) 的倍数，(y,z) 都是这样子的：

[sum_{u=1}^{min(A,B)}sum_{v=1}^{min(A,C)}sum_{w=1}^{min(B,C)}mu(u)mu(v)mu(w)f(lcm(u,v),A)f(lcm(u,w),B)f(lcm(v,w),C) ]

其中 (f(a,b)=sum_{a|x}frac{b}{x}) ，这个东西可以很轻松的 (O(nlog n)) 预处理。但是你发现我们推了这么久的柿子没有什么卵用，时间复杂度还是 (O(n^3)) 的，这时候不得不讲一点前置知识了。

---

感谢这位大佬的博客：Dance of Faith

传统的三元环计数是 (O(nm)) ，也就是我们枚举一条边，再枚举一个点看能否构成环。

我们把原来的无向图改成有向的，每个点定义一个双关键字 ((deg,id)) ，分别表示其度数和编号，每一条边 ((u,v)) 就从值大的连向值小的，现在所有的环都长成这个样子：

这时候统计方法也就呼之欲出了，我们枚举一条有向边，枚举终点所有连出去的边，再看这个点和起点是否构成了三元环，你可能觉得这种做法是辣鸡，对复杂度没有本质的优化，但我们来仔细分析一下复杂度：

我们称出度 (out_ugeqsqrt m) 的点为大度点，否则为小度点，考虑一条边 ((u,v)) 的情况：

• 如果 (v) 是小度点，那没事了，复杂度直接 (O(msqrt m))
• 如果 (v) 是大度点，必然有 (deg_ugeq deg_vgeqsqrt m) ，所以这样的 (u) 最多只有 (O(sqrt m)) 个，那么考虑每个 (v) 的贡献至多是 (O(sqrt mout_v))，因为每个 (v) 至多被 (sqrt m) 个这样的 (u) 相连，因为 (sum out_vleq m) ，所以复杂度 (O(msqrt m))

优美的暴力。

---

回到这道题，要相信我们的反演一定是有用的，观察这个柿子，发现 (u,v,z) 的 (mu) 都不为 (0) 的时候才有贡献，所以 ((u,v,z)) 三元组的数量应该没有那么大，而由于 (lcm<A) 之类的限制三元组的数量就更少了。

从图论的角度来思考，我们把 (lcm(u,v)leqmax(A,B,C)) 的点连边，那么问题就变成了三元环计数，首先我们要知道边数，打一个建图程序发现 (A=B=C=100000) 时，边数是 (760741) 的。

直接套我们讲的东西，时间复杂度变成 (O(msqrt m)) ，建边的时候先枚举 (lcm) ，求出它的质数集合然后做子集枚举（先枚举一个集合再枚举他的子集），写的好一点的话复杂度是 (O(n3^6)) ，因为质数集合不会每个都达到 (6) 所以实际上会快很多。

还有一些细节，比如自环是上述算法难以处理的所以要单独拿出来讨论，建边的时候要避免重边和自环。

你以为这就完了？你发现这个复杂度还是有点悬，所以要卡常，比较重要的一点是跑三元环计数的时候用 ( t vector) 存边会快一些，还有千万不能图方便开 ( t #definespace intspace longlong)

#include <cstdio>
#include <vector>
#include <iostream>
using namespace std;
const int M = 100005;
const int N = 800000;
const int MOD = 1e9+7;
#define ll long long
int read()
{
	int x=0,f=1;char c;
	while((c=getchar())<'0' || c>'9') {if(c=='-') f=-1;}
	while(c>='0' && c<='9') {x=(x<<3)+(x<<1)+(c^48);c=getchar();}
	return x*f;
}
int T,A,B,C,cnt,lim,vis[M],p[M],mu[M],f[M][10],num[M];
int m,ans,a[N],b[N],c[N],fa[M],fb[M],fc[M],d[M];
struct node
{
    int v,c;
};vector<node> g[M];
void init(int n)
{
    mu[1]=1;
    for(int i=2;i<=n;i++)
    {
        if(!vis[i])
        {
            p[++cnt]=i;
            mu[i]=-1;
        }
        for(int j=1;j<=cnt && i*p[j]<=n;j++)
        {
            vis[i*p[j]]=1;
            if(i%p[j]==0) break;
            mu[i*p[j]]=-mu[i];
        }
    }
    for(int i=1;i<=n;i++)
    {
        int x=i;
        for(int j=1;j<=cnt && p[j]*p[j]<=x;j++)
        {
            if(i%p[j]) continue;
            f[i][++f[i][0]]=p[j];
            while(x%p[j]==0) x/=p[j];
        }
        if(x>1) f[i][++f[i][0]]=x;
    }
}
int gcd(int a,int b)
{
    return !b?a:gcd(b,a%b);
}
void work()
{
    lim=max(A,max(B,C));m=ans=0;
    for(int i=1;i<=lim;i++)
    {
        g[i].clear();
        fa[i]=fb[i]=fc[i]=d[i]=0;
        for(int j=i;j<=lim;j+=i)
        {
            fa[i]=(fa[i]+A/j)%MOD;
            fb[i]=(fb[i]+B/j)%MOD;
            fc[i]=(fc[i]+C/j)%MOD;
        }
    }
    for(int i=1;i<=lim;i++)//建图
    {
        if(mu[i]==0) continue;
        int s=1<<f[i][0];num[0]=1;
        for(int j=1;j<=f[i][0];j++)
            num[1<<j-1]=f[i][j];
        for(int j=1;j<s;j++)//处理出质数集合对应的数
        {
            int t=j&(-j);
            if(j^t) num[j]=num[j^t]*num[t];
        }
        for(int j=0;j<s;j++)
            for(int k=j;;k=(k-1)&j)
            {
                int u=num[j],v=i/num[k];
                if(u<v) a[++m]=u,b[m]=v,c[m]=i,d[u]++,d[v]++;
                if(k==0) break;
            }
    }
    for(int i=1;i<=m;i++)
    {
        if(d[a[i]]<=d[b[i]]) swap(a[i],b[i]);
        g[a[i]].push_back(node{b[i],c[i]});//度数大的连向度数小的
    }
    //三个元素都不同
    for(int i=1;i<=m;i++)
    {
        int x=a[i],y=b[i];
        for(int j=0;j<g[y].size();j++)
        {
            int z=g[y][j].v,X=c[i],Y=g[y][j].c,t=gcd(x,z);
            if(1ll*x*z/t>lim)//判断(z,x)有没有边
                continue;
            int Z=x/t*z;
            int tmp=((ll)fa[X]*fb[Y]*fc[Z]+(ll)fa[X]*fb[Z]*fc[Y]
            +(ll)fa[Y]*fb[X]*fc[Z]+(ll)fa[Y]*fb[Z]*fc[X]
            +(ll)fa[Z]*fb[X]*fc[Y]+(ll)fa[Z]*fb[Y]*fc[X])%MOD;
            ans=(ans+mu[x]*mu[y]*mu[z]*tmp)%MOD;
        }
    }
    //其中两个相同
    for(int i=1;i<=m;i++)
    {
        int x=a[i],y=b[i],t=c[i];
        int tmp=((ll)fa[x]*fb[t]*fc[t]+(ll)fa[t]*fb[x]*fc[t]+(ll)fa[t]*fb[t]*fc[x])%MOD;
        ans=(ans+mu[y]*tmp)%MOD;
        tmp=((ll)fa[y]*fb[t]*fc[t]+(ll)fa[t]*fb[y]*fc[t]+(ll)fa[t]*fb[t]*fc[y])%MOD;
        ans=(ans+mu[x]*tmp)%MOD;
    }
    //三个都相同
    for(int i=1;i<=lim;i++)
        ans=(ans+(ll)mu[i]*fa[i]*fb[i]*fc[i])%MOD;
    printf("%d
",(ans+MOD)%MOD);
}
signed main()
{
    init(100000);
    T=read();
    while(T--)
    {
        A=read();B=read();C=read();
        work();
    }
}