BZOJ5332: [Sdoi2018]旧试题(莫比乌斯反演)

时光匆匆，转眼间又是一年寒暑……

这是小 Q 同学第二次参加省队选拔赛。

今年，小 Q 痛定思痛，不再冒险偷取试题，而是通过练习旧试题提升个人实力。可是旧试题太多了，小 Q 没日没夜地做题，却看不到前方的光明在哪里。

一天，因做题过度而疲惫入睡的小 Q 梦到自己在考场上遇到了一道好像做过的题目，却怎么也想不起曾经自己是怎么解决它的，直到醒来还心有余悸。

小 Q 眉头一皱，感觉事情不妙，于是他找到了你，希望你能教他解决这道题目。小 Q 依稀记得题目要计算如下表达式的值

　　$({sum_{i=1}^{A}}{sum_{j=1}^{B}}{sum_{k=1}^{C}}d(ijk))mod(10^9+7)$

其中 d(ijk) 表示 ijk 的约数个数。

Input

第一行包含一个正整数 T，表示有 T 组测试数据。接下来 T 行，每行描述一组测试数据，包含三个整数 A, B 和 C，含义见题目描述。

Output

对于每组测试数据，输出一行，包含一个整数，表示所求表达式的值。

解题思路：

极限卡常题，具体卡常的方法主要是使用vector代替链式前向星，如果还过不去，还是看看dalao们的博客吧。

这里主要讲算法。

首先是喜闻乐见的莫比乌斯反演，三个$sum$有点中暑，慢慢推吧其实都差不多。

不妨设$Aleq Bleq C$,则:

$Ans={sum_{i=1}^{A}}{sum_{j=1}^{B}}{sum_{k=1}^{C}}d(ijk)$

$={sum_{i=1}^{A}}{sum_{x|i}}{sum_{j=1}^{B}}{sum_{y|j}}{sum_{k=1}^{C}}{sum_{z|k}}{varepsilon(gcd(x,y))}{varepsilon(gcd(y,z))}{varepsilon(gcd(x,z))}$

$={sum_{x=1}^{A}}{sum_{x|i}^{A}}{sum_{y=1}^{B}}{sum_{y|j}^{B}}{sum_{z=1}^{C}}{sum_{z|k}^{C}}{varepsilon(gcd(x,y))}{varepsilon(gcd(y,z))}{varepsilon(gcd(x,z))}$

$={sum_{x=1}^{A}}{sum_{y=1}^{B}}{sum_{z=1}^{C}}{leftlfloor{frac{A}{x}} ight floor}{leftlfloor{frac{B}{y}} ight floor}{leftlfloor{frac{C}{z}} ight floor}{varepsilon(gcd(x,y))}{varepsilon(gcd(y,z))}{varepsilon(gcd(x,z))}$

$={sum_{x=1}^{A}}{sum_{y=1}^{B}}{sum_{z=1}^{C}}{leftlfloor{frac{A}{x}} ight floor}{leftlfloor{frac{B}{y}} ight floor}{leftlfloor{frac{C}{z}} ight floor}{sum_{i|gcd(x,y)}}{sum_{j|gcd(y,z)}}{sum_{k|gcd(x,z)}}{mu(i))}{mu(j))}{mu(k))}$

$={sum_{i=1}^{A}}{sum_{j=1}^{B}}{sum_{k=1}^{A}}{mu(x)}{mu(y)}{mu(z)}{sum_{x|lcm(i,k)}}{sum_{y|lcm(i,j)}}{sum_{z|lcm(j,k)}}{left lfloor {frac{A}{i}} ight floor}{left lfloor {frac{B}{j}} ight floor}{left lfloor {frac{C}{k}} ight floor}$

将$mu$不等于0的连边，双向改单向枚举三元环就好了。

代码：

  1 #include<cstdio>
  2 #include<vector>
  3 #include<cstring>
  4 #include<algorithm>
  5 typedef long long lnt;
  6 const int N=100010;
  7 const lnt mod=(lnt)(1e9+7);
  8 struct edge{
  9     int from;
 10     int to;
 11     int lcm;
 12 }ed[N*100];
 13 struct ent{
 14     int twd;
 15     int lcm;
 16 };
 17 int prime[N];
 18 int miu[N];
 19 int ind[N];
 20 bool vis[N];
 21 lnt Lcm[N];
 22 lnt F[3][N];
 23 int cnt;
 24 std::vector<ent>hd[N];
 25 int A,B,C;
 26 void gtp(void)
 27 {
 28     miu[1]=1;
 29     for(int i=2;i<N;i++)
 30     {
 31         if(!vis[i])
 32         {
 33             prime[++cnt]=i;
 34             miu[i]=-1;
 35         }
 36         for(int j=1;j<=cnt&&prime[j]*i<N;j++)
 37         {
 38             vis[i*prime[j]]=true;
 39             if(i%prime[j]==0)
 40                 break;
 41             miu[i*prime[j]]=-miu[i];
 42         }
 43     }
 44     return ;
 45 }
 46 lnt gcd(lnt x,lnt y)
 47 {
 48     if(!y)
 49         return x;
 50     return gcd(y,x%y);
 51 }
 52 int main()
 53 {
 54     int T;
 55     gtp();
 56     scanf("%d",&T);
 57     while(T--)
 58     {
 59         int a,b,c;
 60         scanf("%d%d%d",&a,&b,&c);
 61         A=std::min(a,std::min(b,c));
 62         C=std::max(a,std::max(b,c));
 63         B=a+b+c-A-C;
 64         cnt=0;
 65         for(int i=1;i<=C;i++)
 66             hd[i].clear(),ind[i]=0;
 67         for(int i=1;i<=C;i++)
 68             F[0][i]=F[1][i]=F[2][i]=0;
 69         for(int i=1;i<=A;i++)
 70             for(int j=i;j<=A;j+=i)
 71                 F[0][i]+=A/j;
 72         for(int i=1;i<=B;i++)
 73             for(int j=i;j<=B;j+=i)
 74                 F[1][i]+=B/j;
 75         for(int i=1;i<=C;i++)
 76             for(int j=i;j<=C;j+=i)
 77                 F[2][i]+=C/j;
 78         lnt ans=0;
 79         for(int i=1;i<=A;i++)
 80             if(miu[i])
 81                 ans+=miu[i]*F[0][i]*F[1][i]*F[2][i];
 82         ans%=mod;
 83         for(int d=1;d<=C;d++)
 84         {
 85             for(int i=1;i*d<=C;i++)
 86             {
 87                 if(!miu[i*d])
 88                     continue;
 89                 for(int j=i+1;1ll*j*i*d<=C;j++)
 90                 {
 91                     if(miu[j*d]==0||gcd(j,i)!=1)
 92                         continue;
 93                     int x=i*d,y=j*d,lcm=i*j*d;
 94                     ans+=miu[x]*(F[0][y]*F[1][lcm]*F[2][lcm]+F[1][y]*F[0][lcm]*F[2][lcm]+F[2][y]*F[1][lcm]*F[0][lcm]);
 95                     ans=ans%mod;
 96                     ans+=miu[y]*(F[0][x]*F[1][lcm]*F[2][lcm]+F[1][x]*F[0][lcm]*F[2][lcm]+F[2][x]*F[1][lcm]*F[0][lcm]);
 97                     ans=ans%=mod;
 98                     ind[x]++;
 99                     ind[y]++;
100                     ed[++cnt]=(edge){x,y,lcm};
101                 }
102             }
103         }
104         for(int i=1;i<=cnt;i++)
105         {
106             int f=ed[i].from,t=ed[i].to,v=ed[i].lcm;
107             if(ind[f]<ind[t]||(ind[f]==ind[t]&&f<t))
108                 hd[f].push_back((ent){t,v});
109             else
110                 hd[t].push_back((ent){f,v});
111         }
112         for(int i=1;i<=C;i++)
113         {
114             for(int j=0;j<hd[i].size();j++)
115                 Lcm[hd[i][j].twd]=hd[i][j].lcm;
116             for(int j=0;j<hd[i].size();j++)
117             {
118                 int a=i,b=hd[i][j].twd;
119                 for(int k=0;k<hd[b].size();k++)
120                 {
121                     int c=hd[b][k].twd;;
122                     int lab=hd[a][j].lcm;
123                     int lbc=hd[b][k].lcm;
124                     int lac=Lcm[c];
125                     if(!lac)
126                         continue;
127                     lnt tmp=miu[a]*miu[b]*miu[c];
128                     lnt ttt=F[0][lab]*F[1][lbc]*F[2][lac]+
129                             F[0][lab]*F[1][lac]*F[2][lbc]+
130                             F[0][lbc]*F[1][lac]*F[2][lab]+
131                             F[0][lbc]*F[1][lab]*F[2][lac]+
132                             F[0][lac]*F[1][lab]*F[2][lbc]+
133                             F[0][lac]*F[1][lbc]*F[2][lab];
134                     ans=(ans+ttt*tmp%mod)%mod;
135                 }
136             }
137             for(int j=0;j<hd[i].size();j++)
138                 Lcm[hd[i][j].twd]=0;
139         }
140         printf("%lld
",(ans%mod+mod)%mod);
141     }
142     return 0;
143 }