约定：
((i,j))表示(gcd(i,j)==1)
([x])表示对(x)向下取整
首先有：
于是题目所求即：
然后我们肯定要先枚举因数，就有：
这个时候再反演一波，得到：
我们肯定要换个枚举顺序
于是有：
那个同时整除两个太丑了，我们换成(u|lcm(d,h))算了
接着有：
我们发现后面那一大坨貌似都之和两个量有关
(A,lcm(d,h))这类变量
我们把其中核心的式子再改写试试：
我们考虑记(f(n,x))表示：
注意到一个引理：
于是好像有：
也就是：
所以我们记(f(x)=sum_{i=1}^x[dfrac{x}{i}])
就有原题所求即：
好像(f(x))可以单次(sqrt x)的求出来
我们考虑更快速的推好了
注意到对于某一个数(ile x)，其对于(f(x))会产生的贡献为(dfrac{x}{i})即(1-x)内有多少个数为其倍数
所以(f)的实际意义即(1-x)内的约数和
好了这样你就可以(O(nlog n)|O(n))把这些值推出来了
当然尽管反演了这么多你的复杂度还是(O(n^3log n))的qwq
我们发现我们要求的实际上是有序对((d,k,h))对答案造成的贡献
于是可以考虑转化问题的模型，建出一张图来
考虑给两个点((u,v))之间连一条边
这是一条无向边，这条边的边权为(lcm(u,v))，我们记边(u,v)的边权为(w_{u,v})
于是好像对于一个有序对((d,k,h))，若其两两不相同，则其可以通过这张图上的三元环来表示
对于一个三元环((d,k,h))其对答案造成贡献当且仅当(dle min(A,B),kle min(B,C),hle min(C,A))
那么它会对答案造成(mu(d)*mu(k)*mu(h)*f(dfrac{A}{w_{d,h}})f(dfrac{B}{w_{d,k}})f(dfrac{C}{w_{h,k}}))
于是我们好像就可以在这张图上做三元环计数然后统计答案了qwq
但好像这样还是很暴力
我们尝试给这张图剪一下枝
发现大量的三元环对答案造成的贡献都是(0)
于是我们可以把(mu(x)=0)或者(lcm(u,v)>max(A,B,C))的点都删去
据shadowice1984聚聚说的这样剪之后原图上的边会很少
当然，暴力建边的复杂度是(O(n^2))的————(Step(1))
我们注意到这个定义下建出的边只能统计两两不相等的三元组对答案造成的贡献
我们就考虑暴力统计，每次枚举一个点并枚举其出点和其出点的出点，这样做复杂度是(O(nm))的————(Step(2))
于是你还需要额外写个特判去判断两两相等的三元组对答案造成的贡献
可以发现算两两相等的三元组的复杂度是(O(m))而计算三个全部相等的复杂度则是(O(n))
我们尝试把(5000)的点的部分的边数数出来，发现它是 (49528)
如果你写了个(map)来判两点之间有没有连边肯定会T飞以至于一分都无法获得
否则的话如果你不带(log)而是单纯的非严格(nm)你应该可以获得(30pts)的高分
于是我们要考虑优化
(1.)优化(Step2)
于是我们把三元环计数中的技巧拿出来用，就可以将复杂度变成(msqrt m)了
当然三元环计数统计的是无向三元环计数
所以你还要枚举其六种不同的排列来确保正确性
这样就可以将复杂度做到(O(n^2+msqrt m))了
如果你尝试写了一发并交了上去你会惊人的发现这个做法好像还是(30)分
(2.)优化(Step1)
我们考虑不以(n^2)的枚举来建边
考虑枚举一个(gcd=x)，那么(i=kx,j=px)
则(gcd(k,p)==1)
且(mu(i) e0,mu(j) e0)并且((kpx)le max(A,B, C))
我们枚举这个点对((k,p))就行了
当然这样的建边复杂度并不是(O(m))
但也比(O(n^2))快，而且快得多
关于近似(O(m))的建边就是考虑枚举(lcm)，同时(mu(lcm) e0)，然后质因数分解其，，枚举可能的(mu(x) e0)的点，对于每个点可以通过一个关于因数的取或不取二进制去标号它，记录一个桶，然后枚举其子集，那么要连的边就是其子集与其补集(&)得到的结果
于是你就十分开心自以为能过的打出了这道题的代码
当然如果你用的是vector你是真的能过的
否则你应该会(T)
关于代码中特判两个相等的三元组对答案的贡献，因为这个菜得要死的笔者是先写完暴力再把两个优化加上的，所以他没有在连边的时候把这个情况处理掉
实际上是可以在连边的过程中处理掉的qwq
小提示(1:) 貌似枚举三元组的时候受到的限制很麻烦的样子(min(A,B),min(B,C),min(C,A))
然而实际上你并不需要管，统一把限制设成(max(A,B,C))就好了
这是因为如果(u>min(A,B))那么(f(min(A,B)/u)=0)
下面是又丑又长的代码
#include<bits/stdc++.h>
using namespace std;
#define rep( i, s, t ) for( register int i = s; i <= t; ++ i )
#define re register
#define int long long
int read() {
char cc = getchar(); int cn = 0, flus = 1;
while(cc < '0' || cc > '9') { if( cc == '-' ) flus = -flus; cc = getchar(); }
while(cc >= '0' && cc <= '9') cn = cn * 10 + cc - '0', cc = getchar();
return cn * flus;
}
const int N = 100000 + 500 ;
const int P = 1e9 + 7 ;
int A, B, C, f[N], u[N], p[N], d[N], r[N], deg[N], vis[N], book[N], cnt, tot, top ;
bool isp[N] ;
struct E { int to, w ; };
struct node { int u, v, w ; } e[N * 10];
vector<E> mp[N] ;
inline void init() {
isp[1] = 1, u[1] = 1, r[++ tot] = 1, book[1] = tot ;
for( re int i = 2; i <= N - 5; ++ i ) {
if( !isp[i] ) p[++ top] = i, u[i] = -1 ;
for( re int j = 1; j <= top && i * p[j] <= N - 5; ++ j ) {
isp[i * p[j]] = 1 ;
if( i % p[j] == 0 ) continue ;
u[i * p[j]] = - u[i] ;
}
if( u[i] != 0 ) r[++ tot] = i, book[i] = tot ;
}
for( re int i = 1; i <= N - 5; ++ i ) {
for( re int j = i; j <= N - 5; j += i ) ++ d[j] ;
f[i] = ( f[i - 1] + d[i] ) % P ;
}
}
inline int gcd( int x, int y ) {
if( x == 0 ) return y ;
return gcd( y % x, x ) ;
}
signed main()
{
int T = read() ; init() ;
while( T-- ) {
memset( mp, 0, sizeof(mp) ), memset( deg, 0, sizeof(deg) ) ;
A = read(), B = read(), C = read(), cnt = 0 ;
int Ans = 0, D = 0 ; D = max( A, max( B, C ) ) ;
//连边
for( re int g = 1; r[g] <= D; ++ g ) {
for( re int i = 1; r[i] * r[g] <= D; ++ i ) {
if( !u[r[i] * r[g]] ) continue ;
for( re int j = i + 1; j <= tot && r[i] * r[j] * r[g] <= D; ++ j ) {
if( !u[r[j] * r[g]] || gcd( r[i], r[j] ) != 1 ) continue ;
int u = book[r[i] * r[g]], v = book[r[j] * r[g]], w = r[i] * r[j] * r[g] ;
++ deg[u], ++ deg[v], e[++ cnt] = (node){ u, v, w } ;
mp[u].push_back((E){v, w}), mp[v].push_back((E){u, w}) ;
}
}
}
//两个相等
int rD = min( A, min( B, C ) ) ;
for( re int i = 1; r[i] <= min( A, B ); ++ i ) {
int siz = mp[i].size() - 1 ;
rep( j, 0, siz ) {
int u1 = mp[i][j].to ;
Ans += u[r[i]] * u[r[u1]] * u[r[u1]] * f[A / mp[i][j].w] * f[B / mp[i][j].w] * f[C / r[u1]] ;
Ans += u[r[i]] * u[r[i]] * u[r[u1]] * f[A / r[i]] * f[B / mp[i][j].w] * f[C / mp[i][j].w] ;
Ans += u[r[i]] * u[r[i]] * u[r[u1]] * f[A / mp[i][j].w] * f[B / r[i]] * f[C / mp[i][j].w] ;
}
}
// 三个相等
for( re int i = 1; r[i] <= rD; ++ i ) Ans += u[r[i]] * u[r[i]] * u[r[i]] * f[A / r[i]] * f[B / r[i]] * f[C / r[i]] ;
//互不相等
memset( mp, 0, sizeof(mp) ) ;
for( re int i = 1; i <= cnt; ++ i )
if( deg[e[i].u] >= deg[e[i].v] ) mp[e[i].u].push_back((E){ e[i].v, e[i].w }) ;
else mp[e[i].v].push_back((E){ e[i].u, e[i].w }) ;
for( re int i = 1; r[i] <= D; ++ i ) {
int siz = mp[i].size() - 1 ;
rep( j, 0, siz ) vis[mp[i][j].to] = mp[i][j].w ;
rep( j, 0, siz ) {
int u1 = mp[i][j].to, siz2 = mp[u1].size() - 1 ;
rep( k, 0, siz2 ) {
int v = mp[u1][k].to ;
if( !vis[v] ) continue ;
int muu = u[r[i]] * u[r[u1]] * u[r[v]], b = mp[i][j].w, c = mp[u1][k].w, a = vis[v] ;
int rc = f[A / a] * ( f[B / b] * f[C / c] + f[B / c] * f[C / b] )
+ f[A / b] * ( f[B / a] * f[C / c] + f[B / c] * f[C / a] )
+ f[A / c] * ( f[B / a] * f[C / b] + f[ B / b ] * f[C / a] ) ;
Ans += muu * rc ;
}
}
rep( j, 0, siz ) vis[mp[i][j].to] = 0 ;
}
printf("%d
", Ans % P ) ;
}
return 0;
}