自动搬运

查看原文

来自洛谷，原作者为

	shadowice1984 これが勘違いでも、これが勘違いでも

搬运于2025-08-24 22:02:01，当前版本为作者最后更新于2018-05-22 18:18:43，作者可能在搬运后再次修改，您可在原文处查看最新版

自动搬运只会搬运当前题目点赞数最高的题解，您可前往洛谷题解查看更多

以下是正文

---

神奇反演题……

考验人类智慧的卡常数……

同样复杂度的算法一个跑30s一个跑2s……

然后相信你的信仰，用力卡常数吧

---

本题题解

前置芝士:莫比乌斯反演

请确保你熟练掌握了交换Σ那一套本题解不再讲反演的基础了

---

题目要求

$\sum_{i=1}^{A}\sum_{j=1}^{B}\sum_{k=1}^{C}\sigma_{0}(ijk)$

其中$\sigma_{0}(n)=\sum_{d|n}1$

然后这是一道反演除数函数$\sigma$的题，那么这里其实对于除数函数里边的乘积是有一个套路的，枚举来自于各个因子的互质约数对，然后直接计算

换句话说，对于在函数里边的ijk乘积来讲，我们考虑这个ijk的约数的构成，显然是有来自i的，j的，k的部分，那么我们直接枚举约数三元组求lcm显然是不行的，因为会重复一大堆，但是我们观察重复的部分每一个lcm相同的三元组集当中总是有一个三个约数两两互质的部分，因此可以只枚举两两互质的约数三元组……

说了这么多，只是为了说明刚才的式子可以被重写为下面的表达式

$\sum_{i=1}^{A}\sum_{j=1}^{B}\sum_{k=1}^{C}\sum_{d|i}\sum_{t|j}\sum_{p|k}\epsilon((d,t))\epsilon((t,p))\epsilon((d,p))$

中$\epsilon(x)$是元函数，当且仅当$x=1\epsilon(x)=1$，其他情况为0，$(i,j)$是$gcd(i,j)$的简记形式

下面是喜闻乐见的交换Σ环节！

先把这个式子重写成这样:

$\sum_{i=1}^{A}\sum_{d|i}\sum_{j=1}^{B}\sum_{t|j}\sum_{k=1}^{C}\sum_{p|k}\epsilon((d,t))\epsilon((t,p))\epsilon((d,p))$

然后又是经典的前边是$1~n$顺次枚举后边是$d|i$约数枚举的形式，我们令$i/=d$则$i=di$，交换求和号，得到

$\sum_{d=1}^{A}\sum_{i=1}^{\lfloor\frac{A}{d}\rfloor}\sum_{j=1}^{B}\sum_{t|j}\sum_{k=1}^{C}\sum_{p|k}\epsilon((d,t))\epsilon((t,p))\epsilon((d,p))$

同理我们也可以用相同的手法处理j与t,k与p，得到

$\sum_{d=1}^{A}\sum_{i=1}^{\lfloor\frac{A}{d}\rfloor}\sum_{t=1}^{B}\sum_{j=1}^{\lfloor\frac{B}{t}\rfloor}\sum_{p=1}^{C}\sum_{k=1}^{\lfloor\frac{C}{p}\rfloor}\epsilon((d,t))\epsilon((t,p))\epsilon((d,p))$

简单整理下得到

$\sum_{d=1}^{A}\sum_{t=1}^{B}\sum_{p=1}^{C}\lfloor\frac{A}{d}\rfloor\lfloor\frac{B}{t}\rfloor\lfloor\frac{C}{p}\rfloor\epsilon((d,t))\epsilon((t,p))\epsilon((d,p))$

然后我们运用莫比乌斯反演4个推论中的一个，使用$\mu$来反演$\epsilon$得到

$\sum_{d=1}^{A}\sum_{t=1}^{B}\sum_{p=1}^{C}\lfloor\frac{A}{d}\rfloor\lfloor\frac{B}{t}\rfloor\lfloor\frac{C}{p}\rfloor\sum_{u|t\&u|d}\mu(u)\sum_{v|t\&v|p}\mu(v)\sum_{w|d\&w|p}\mu(w)$

好了又是大家喜闻乐见的交换求和号的环节~

然后这一步交换求和号和其他的交换求和号不同的是，我们必须同时交换3对求和号，此时原来常用的令某一个变量做一个变换的手法已经失效了

此时交换求和号的准则就是，保证每一个单独的变量被枚举了相同的次数

观察限制条件，发现$u|d\&w|d$,$u|t\&v|t$,$v|p\&w|p$

然后直接按照这个变量的限制条件交换求和号会得到

$\sum_{u=1}^{min(A,B)}\sum_{v=1}^{min(B,C)}\sum_{w=1}^{min(A,C)}\mu(u)\mu(v)\mu(w)\sum_{u|d\&w|d}\lfloor\frac{A}{d}\rfloor\sum_{u|t\&v|t}\lfloor\frac{B}{t}\rfloor\sum_{v|p\&w|p}\lfloor\frac{C}{p}\rfloor$

由于$u|d\&w|d$其实就是$lcm(u,w)|d$，所以重新写一下表达式就是

$\sum_{u=1}^{min(A,B)}\sum_{v=1}^{min(B,C)}\sum_{w=1}^{min(A,C)}\mu(u)\mu(v)\mu(w)\sum_{lcm(u,w)|d}\lfloor\frac{A}{d}\rfloor\sum_{lcm(u,v)|t}\lfloor\frac{B}{t}\rfloor\sum_{lcm(v,w)|p}\lfloor\frac{C}{p}\rfloor$

那么函数$f(n,x)=\sum_{n|d}\lfloor\frac{x}{d}\rfloor$其实是可以$O(nlogn)$的枚举约数对于某一个参数x处理出来的，因此刚才的式子可以写成

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

到此为止反演的工作已经做完了，下面请让我们暂时忘记我们会分析算法复杂度这个事实，因为接下来你会觉得这个算法的算法复杂度十分魔幻……

暴力1

我会$O(n^3)$枚举u,v,w有序三元组！

恭喜你白反演了……

暴力2

我们考虑给刚才的暴力剪一剪枝，我们发现一个非常有趣的事实，就是如果n>x那么$f(n,x)=0$那么我们可以根据这个性质建出一张图来,在这张图中当中两个点$u,v$有边当且仅当$\mu(u)!=0,\mu(v)!=0,lcm(u,v)<max(A,B,C)$,然后连一条边权为lcm的边

这样的话我们会发现一个惊人的事实，暴力1中所枚举的三元组可能有效当且仅当这个三元组在这张图当中构成了一个三元环……于是我们开始测试这个剪枝的效率

但是首先我们该如何$O(m)$而不是$O(n^3)$的建图呢?

答案是枚举边权lcm，我们一次性连出边权为i的所有边，那么首先我们可以$O(nloglogn)$的筛出每一个数的互异质因子

下面我们从1到max(A,B,C)枚举边权i，首先因为这个边权是两个$\mu$值不为0 的数的lcm，因此$\mu(i)!=0$，下面我们直接$2^k$的枚举这个i的质因数分解集合的一个子集，由这个子集我们可以得到其中的一个点的编号u，之后我们要确定另外一个点v但是我们发现v不是非常的好枚举了，因此我们可以枚$gcd(u,v)$显然这个gcd一定是u对应的集合的子集，因此我们再次枚举u的子集 确定出gcd，那么对应的v就是$i×gcd/u$了之后我们就在u和v间连一条边

好了听起来刚才的算法复杂度相当不靠谱，如果记$w(d)$为d的互异质因子个数的话，我们的算法复杂度就是$O(\sum3^{w(d)})$或者$O(m)$的，然后我们尝试着用这个程序跑了一遍$A=B=C=10^5$并且记录一下我们建了多少边，结果可能有些令人吃惊，边数是

$\text{760，741}$

这意味着这张图实际上相当稀疏……，因此我们可以考虑枚举这张图的三元环了

那么这里我们采用暴力进行枚举，暴力的遍历每个的出边，遍历另一个端点的出边，采用哈希表判断第3个点是否和第一个点有连边

显然这个算法的复杂度是不对的，随随便便就卡成$O(nm)$了，然后我们考虑给这个算法加个小剪枝，将边从无向边改为有向边，方向为度数大的点指向度数小的点，令人惊讶的是，算法复杂度骤减至$O(m\sqrt{m})$

下面给出神奇的证明

我们称度数超过$\sqrt{m}$的点为大度点，否则就是小度点，然后我们发现两边端点都是大度点的边不超过$O(\sqrt{m})$条，因此这部分复杂度为$O(m\sqrt{m})$对于其余的边每条边贡献的复杂度都不超过$\sqrt{m}$所以总复杂度为$O(m\sqrt{m})$

那么我们信心满满的写了这个算法复杂度为$O(m\sqrt{m})$的暴力交上去一看T成40分，然后我们觉得没啥大不了的，毕竟是暴力嘛

---

正解

于是你去翻了这道题的题解，题解告诉你恭喜你刚才的暴力就是正解，只是你的三元环枚举可能有点不优美于是就华丽的T飞了

你又去刚才的程序跑了一遍$10^5$的点，大概需要12s……

可是我们一个点最多跑2.5s……

于是我们开始了挑战人类智慧的卡常数之旅……

首先我们发现事实上哈希表是不必要的，我们可以直接手动计算两个点之间的lcm从而判断两个点之间是否有连边，于是你信心满满的加上了这个优化，发现没什么卵用的快了2s

接下来我们发现枚举三元组其实是非常消耗计算资源的，于是你去掉了原来的图中所有的自环，拉出来单独枚举，然后你发现这个优化的非常给力，你的程序又跑快了1s

然后我们发现我们事实上是枚举了所有的有序3元组，于是你决定不在遍历到一条边的时候临时决定这个边的方向，而是一开始就把这个边的方向确定好

但是这样会有一个非常诡异的问题就是我们的枚举现在只能枚举无序三元组了…… 因此我们求一个3元组的贡献必须计算6种不同的贡献方式，但是这个优化的效果还是有的，你的程序跑到了5s

于是我们信心满满的修改了一些奥妙重重的运算方式，从而减少了一些取膜和加法运算，很好你的程序(并无卵用的)跑到了4.2s

你对着这份代码又卡了半天常数，发现并无卵用，正当你感叹到底啥才是优美，为什么这么卡常数，按下$Ctrl+A$准备重构代码的时候，你突然看到了你的遍历出边的代码，那是再平常不过的一段代码，大家常用的

for(int i=alist[u];i;i=nxt[i]){}

你眉头一皱发现事情并不对劲，因为这份代码在遍历出边的时候访问的内存并不连续而且跨度还不是一般的大……，这意味着你的程序将会发生$O(m\sqrt{m})$次cache miss，这部分的常数就要占去2s……

于是你快速的把邻接链表换成大家常用的vector，这样你在遍历一个点的出边的时候将会快速访问一段连续的内存，cache miss大大减少，现在你的程序可以在2.2s内完成任务了~

当然代价就是你的代码基本不能看了，所以代码仅供参考

上代码~

#pragma GCC optimize(3)
#include<cstdio>
#include<algorithm>
#include<ctime> 
#include<vector>
using namespace std;const int N=1e5+10;typedef long long ll;const ll mod=1e9+7;
const int E=1e6+10;const int U=1e5;int zhi[N];bool book[N];int miu[N];int hd;int T;int d[N];
int a;int b;int c;int l1;int l2;int l3;int lim;int fj[N][15];int tp[N];int st;int edt;
ll ret=0;int fa[N];int fc[N];int fb[N];struct data{int v;int val;};vector <data> ed[N];
inline int gcd(int a,int b){int c;while(b){c=a%b;a=b;b=c;}return a;}
int eu[E];int ev[E];int et[E];int fr;
inline void solve()
{
    scanf("%d%d%d",&a,&b,&c);lim=max(max(a,b),c);st=clock();
    if(a>b)swap(a,b);if(a>c)swap(a,c);if(b>c)swap(b,c);
    for(int i=1;i<=a;i++)for(int j=1;j*i<=a;j++)(fa[i]+=a/(i*j))%=mod;
    for(int i=1;i<=b;i++)for(int j=1;j*i<=b;j++)(fb[i]+=b/(i*j))%=mod;
    for(int i=1;i<=c;i++)for(int j=1;j*i<=c;j++)(fc[i]+=c/(i*j))%=mod;
    for(int i=1;i<=lim;i++)//建图
    {
        if(miu[i]==0)continue;
        for(int j=0;j<=(1<<tp[i])-1;j++)
        {
            int ret=1;for(int k=0,p=j;p;p>>=1,k++)if(p&1)ret*=fj[i][k];
            for(int k=j;;k=(k-1)&j)
            {
                int g=1;for(int t=0,p=k;p;p>>=1,t++)if(p&1)g*=fj[i][t];
                int v=(ll)i*g/ret;d[ret]++;d[v]++;
                if(ret<v){eu[++fr]=ret;ev[fr]=v;et[fr]=i;}if(k==0)break;//去掉重边
            }
        }
    }
    for(int i=1;i<=fr;i++){if(d[eu[i]]<d[ev[i]])swap(eu[i],ev[i]);ed[eu[i]].push_back((data){ev[i],et[i]});}//变为有向边
    vector <data> ::iterator it1,it2;
    for(int u=1;u<=lim;u++)//vector存边
        for(it1=ed[u].begin();it1!=ed[u].end();++it1)
            for(it2=ed[it1->v].begin();it2!=ed[it1->v].end();++it2)
            {
                int v3;if((v3=(ll)it2->v*u/gcd(it2->v,u))>lim)continue;
                int v1=it1->val;int v2=it2->val;
                if(miu[u]*miu[it1->v]*miu[it2->v]==1)
                {
                    (ret+=((ll)fb[v2]*fc[v3]+(ll)fb[v3]*fc[v2])*fa[v1]+
                          ((ll)fb[v1]*fc[v3]+(ll)fb[v3]*fc[v1])*fa[v2]+
                          ((ll)fb[v2]*fc[v1]+(ll)fb[v1]*fc[v2])*fa[v3])%=mod;//处理三元环贡献，基本不能看了
                }
                else 
                {
                    (ret-=((ll)fb[v2]*fc[v3]+(ll)fb[v3]*fc[v2])*fa[v1]+
                          ((ll)fb[v1]*fc[v3]+(ll)fb[v3]*fc[v1])*fa[v2]+
                          ((ll)fb[v2]*fc[v1]+(ll)fb[v1]*fc[v2])*fa[v3])%=mod;
                    ret>0?ret:ret+mod;
                }
            }	
    for(int i=1;i<=fr;i++)//有一个重复点的三元环
    {
        if(miu[eu[i]]==1)
            (ret+=(ll)fa[et[i]]*fb[et[i]]%mod*fc[ev[i]]+(ll)fa[et[i]]*fb[ev[i]]%mod*fc[et[i]]+(ll)fa[ev[i]]*fb[et[i]]%mod*fc[et[i]])%=mod;
        else 
            (ret-=(ll)fa[et[i]]*fb[et[i]]%mod*fc[ev[i]]+(ll)fa[et[i]]*fb[ev[i]]%mod*fc[et[i]]+(ll)fa[ev[i]]*fb[et[i]]%mod*fc[et[i]])%=mod;
        ret=ret>0?ret:ret+mod;
        if(miu[ev[i]]==1)
            (ret+=(ll)fa[et[i]]*fb[et[i]]%mod*fc[eu[i]]+(ll)fa[et[i]]*fb[eu[i]]%mod*fc[et[i]]+(ll)fa[eu[i]]*fb[et[i]]%mod*fc[et[i]])%=mod;
        else 
            (ret-=(ll)fa[et[i]]*fb[et[i]]%mod*fc[eu[i]]+(ll)fa[et[i]]*fb[eu[i]]%mod*fc[et[i]]+(ll)fa[eu[i]]*fb[et[i]]%mod*fc[et[i]])%=mod;
        ret=ret>0?ret:ret+mod;
    }
    for(int i=1;i<=min(min(a,b),c);i++)//三个重复点的三元环
    {
        if(miu[i]==1)(ret+=(ll)fa[i]*fb[i]%mod*fc[i])%=mod;
        else if(miu[i]==-1)(ret-=(ll)fa[i]*fb[i]%mod*fc[i])%=mod;ret=ret>0?ret:ret+mod;
    }printf("%lld\n",ret);
}
inline void clear()
{
    for(int i=1;i<=a;i++)fa[i]=0;
    for(int i=1;i<=b;i++)fb[i]=0;
    for(int i=1;i<=c;i++)fc[i]=0;
    for(int i=1;i<=lim;i++){vector <data> emp;swap(emp,ed[i]);}
    ret=0;fr=0;
}
int main()
{
    miu[1]=1;for(int i=2;i<=U;i++)//线性筛莫比乌斯函数
    {
        if(book[i]==false){zhi[++hd]=i;miu[i]=-1;}
        for(int j=1;j<=hd&&zhi[j]*i<=U;j++)
        {book[i*zhi[j]]=true;if(i%zhi[j]==0){break;}miu[i*zhi[j]]=miu[i]*-1;}
    }
    for(int i=1;i<=hd;i++)//筛质因数分解
        for(int j=1;j*zhi[i]<=U;j++)
            {int nw=j*zhi[i];fj[nw][tp[nw]]=zhi[i];++tp[nw];}
    scanf("%d",&T);for(int z=1;z<=T;z++){solve();clear();}return 0;//拜拜程序~
}