「SDOI2018」旧试题
题目
点这里看题目。
给定 \(A,B,C\)，求：
\[\sum_{i=1}^A\sum_{j=1}^B\sum_{k=1}^C\sigma_0(ijk)
\]
单个测试点内有 \(T\) 组测试数据。
所有测试点满足 \(1\le T\le 10,1\le A,B,C\le 10^5,1\le \sum \max(A,B,C)\le 2\times 10^5\)。
分析
处理乘积因子个数求和还是挺经典的，手法比较通用：考虑 \(ijk\) 的一个因子 \(d\)，我们令 \(i\) 贡献的部分为 \(\gcd(i,d)\)，\(j\) 贡献的部分为 \(\gcd\left(\frac{d}{\gcd(i,d)},j\right)\)，\(k\) 贡献的就是剩下的。
反过来，就相当于枚举 \(i,j,k\) 的因子（分别为 \(d_1,d_2,d_3\)），然后要求 \(\gcd\left(\frac i {d_1},d_2d_3\right)=\gcd\left(\frac j{d_2},d_3\right)=1\)。有了这个就可以开始推式子了：
\[\newcommand{\lcm}{\operatorname{lcm}}
\newcommand{\divide}[2]{\left\lfloor\frac{#1}{#2}\right\rfloor}
\begin{aligned}
&\sum_{i=1}^A\sum_{j=1}^B\sum_{k=1}^C\sigma_0(ijk)\\
=&\sum_{i=1}^A\sum_{j=1}^B\sum_{k=1}^C\sum_{d_1|i}\sum_{d_2|j}\sum_{d_3|k}[\gcd(i/d_1,d_2)=1][\gcd(i/d_1,d_3)=1][\gcd(j/d_2,d_3)=1]\\
=&\sum_{d_1,d_2,d_3}\sum_{i=1}^{\divide A {d_1} }\sum_{j=1}^{\divide{B}{d_2}}\sum_{k=1}^{\divide{C}{d_3}}[\gcd(i,d_2)=1][\gcd(i,d_3)=1][\gcd(j,d_3)=1]\\
=&\sum_{d_2,d_3}\divide{C}{d_3}\sum_{i=1}^{A}\divide{A}{i}\sum_{j=1}^{\divide{B}{d_2}}\sum_{p|d_2,p|i}\mu(p)\sum_{q|d_3,q|i}\mu(q)\sum_{r|d_3,r|j}\mu(r)\\
=&\sum_{d_2,d_3}\sum_{p|d_2}\sum_{q|d_3}\sum_{r|d_3}\divide{C}{d_3}\mu(p)\mu(q)\mu(r)\sum_{i'=1}^{\divide{A}{\lcm(p,q)}}\divide{A}{i'\lcm(p,q)}\divide{B}{d_2r}\\
=&\sum_{p,q,r}\mu(p)\mu(q)\mu(r)\sum_{i'=1}^{\divide{A}{\lcm(p,q)}}\divide{A}{i'\lcm(p,q)}\sum_{j'=1}^{\divide{B}{pr}}\divide{B}{j'pr}\sum_{k'=1}^{\divide{C}{\lcm(q,r)}}\divide{C}{k'\lcm(q,r)}
\end{aligned}
\]
考虑记 \(f(n)=\sum_{i=1}^n\lfloor\frac{n}{i}\rfloor\)，最后的答案就是：
\[\newcommand{\lcm}{\operatorname{lcm}}
\newcommand{\divide}[2]{\left\lfloor\frac{#1}{#2}\right\rfloor}
\sum_{p,q,r}\mu(p)\mu(q)\mu(r)f\left(\divide{A}{\lcm(q,p)}\right)f\left(\divide{B}{pr}\right)f\left(\divide{C}{\lcm(q,r)}\right)
\]
实际上没有用的观察：
\(p,q,r\) 都是 square-free 的，这样的数比较少。
\((p,r)\) 对是 \(O(B\log B)\) 的。
\(\min\{p,r\}\le \sqrt B\)。
这个时候单纯的数论技巧已经不能拯救我们了，得另辟蹊径。
注意到任一项只受两个数影响，而不存在三个数都参与的项，从这一点出发容易想到“连边”表示二元关系。考虑将 \(\mu\neq 0\) 的正整数作为结点放到图上，并在 \(\operatorname{lcm}(u,v)\le \max\{A,B,C\}\) 的点对 \((u,v)\) 之间连上边，则可以产生贡献的 \((p,q,r)\) 对必须是图上的一个三元环。
三元环计数受限于边数。这里的边数是多少呢？因为不能简单地给出较小的上界，经验判断很容易让人放弃这个想法——然而，在 \(\mu\neq 0\) 和 \(\operatorname{lcm}\le \max\{A,B,C\}\) 两条限制下，有效的边数大约在 \(8\times 10^5\) 左右，小得惊人！
Remark. 经验判断常常有效，但不总是值得信赖。此处忽略实际条件而轻信经验，就会误入歧途。
所以，把边建出来后跑三元环计数即可。有几个小点值得一记：
-
用
std :: vector存边，规避邻接表夸张的 cache miss 和 random access。 -
压缩边数，把自环的情况拿出来单独算。三元环计数只跑平凡三元环。
-
此处边之间不等价，所以可能需要将同一个三元环计数多次——复杂度不变，但是常数大增。
所以得转化成边“等价”三元环计数，但这就必须要手动枚举 \(6\) 种贡献情况。
嘛......卡常有效就行。
代码
#include <bits/stdc++.h>
#define rep( i, a, b ) for( int i = (a) ; i <= (b) ; i ++ )
#define per( i, a, b ) for( int i = (a) ; i >= (b) ; i -- )
const int mod = 1e9 + 7;
const int MAXN = 1e5 + 5, MAXE = 8e5 + 5;
template<typename _T>
inline void Read( _T &x ) {
x = 0; char s = getchar(); bool f = false;
while( s < '0' || '9' < s ) { f = s == '-', s = getchar(); }
while( '0' <= s && s <= '9' ) { x = ( x << 3 ) + ( x << 1 ) + ( s - '0' ), s = getchar(); }
if( f ) x = -x;
}
template<typename _T>
inline void Write( _T x ) {
if( x < 0 ) putchar( '-' ), x = -x;
if( 9 < x ) Write( x / 10 );
putchar( x % 10 + '0' );
}
typedef std :: pair<int, int> Edge;
int su[MAXN][3];
int edgFr[MAXE], edgTo[MAXE], edgWei[MAXE];
int tot = 0;
std :: vector<Edge> grph[MAXN];
int deg[MAXN];
int f[MAXN], mu[MAXN];
int A, B, C;
inline int Mul( int x, const int &v ) { return 1ll * x * v % mod; }
inline int Sub( int x, const int &v ) { return ( x -= v ) < 0 ? x + mod : x; }
inline int Add( int x, const int &v ) { return ( x += v ) >= mod ? x - mod : x; }
inline int& MulEq( int &x, const int &v ) { return x = 1ll * x * v % mod; }
inline int& SubEq( int &x, const int &v ) { return ( x -= v ) < 0 ? ( x += mod ) : x; }
inline int& AddEq( int &x, const int &v ) { return ( x += v ) >= mod ? ( x -= mod ) : x; }
inline int Gcd( int x, int y ) { for( int z ; y ; z = x, x = y, y = z % y ); return x; }
inline int Lcm( int x, int y ) { return x / Gcd( x, y ) * y; }
inline void Init( const int &n ) {
for( int i = 1 ; i <= n ; i ++ )
for( int j = i ; j <= n ; j += i )
AddEq( f[j], 1 );
rep( i, 1, n ) AddEq( f[i], f[i - 1] );
mu[1] = 1;
for( int d = 1 ; d <= n ; d ++ )
for( int k = 2 ; k * d <= n ; k ++ )
SubEq( mu[k * d], mu[d] );
}
int main() {
Init( 1e5 );
int T; Read( T );
while( T -- ) {
Read( A ), Read( B ), Read( C );
int lim = std :: max( A, std :: max( B, C ) );
rep( i, 1, lim ) deg[i] = 0, grph[i].clear();
tot = 0;
rep( p, 1, lim ) if( mu[p] ) {
int limQ = lim / p;
rep( q, 1, limQ ) if( mu[p * q] ) {
int limR = lim / p / q;
rep( r, q + 1, limR ) if( mu[p * r] && Gcd( q, r ) == 1 ) {
edgWei[++ tot] = p * q * r;
deg[edgFr[tot] = p * q] ++;
deg[edgTo[tot] = p * r] ++;
}
}
}
rep( i, 1, lim ) grph[i].reserve( deg[i] );
rep( i, 1, tot ) {
int u = edgFr[i], v = edgTo[i];
if( deg[u] == deg[v] ? u > v : deg[u] > deg[v] )
grph[u].emplace_back( v, edgWei[i] );
else grph[v].emplace_back( u, edgWei[i] );
}
int ans = 0;
for( int p = 1 ; p <= A && p <= C && 1ll * p * p <= B ; p ++ )
AddEq( ans, Mul( mu[p], Mul( Mul( f[A / p], f[C / p] ), f[B / p / p] ) ) );
for( int i = 1, p, q, r, w ; i <= tot ; i ++ ) {
w = edgWei[i];
p = q = edgFr[i], r = edgTo[i];
AddEq( ans, Mul( mu[r], Mul( Mul( f[A / p], f[B / p / r] ), f[C / w] ) ) );
p = q = edgTo[i], r = edgFr[i];
AddEq( ans, Mul( mu[r], Mul( Mul( f[A / p], f[B / p / r] ), f[C / w] ) ) );
q = r = edgFr[i], p = edgTo[i];
AddEq( ans, Mul( mu[p], Mul( Mul( f[C / q], f[B / p / r] ), f[A / w] ) ) );
q = r = edgTo[i], p = edgFr[i];
AddEq( ans, Mul( mu[p], Mul( Mul( f[C / q], f[B / p / r] ), f[A / w] ) ) );
p = r = edgFr[i], q = edgTo[i];
AddEq( ans, Mul( mu[q], Mul( Mul( f[A / w], f[C / w] ), f[B / p / r] ) ) );
p = r = edgTo[i], q = edgFr[i];
AddEq( ans, Mul( mu[q], Mul( Mul( f[A / w], f[C / w] ), f[B / p / r] ) ) );
}
for( int u = 1, v, w ; u <= lim ; u ++ ) {
int n = grph[u].size(), m;
for( int i = 0 ; i < n ; i ++ )
memset( su[grph[u][i].first], 0, 12 );
for( int i = 0 ; i < n ; i ++ ) {
v = grph[u][i].first, m = grph[v].size();
if( 1ll * u * v <= B ) {
int coe = Mul( mu[v], f[B / u / v] );
for( int j = 0, val ; j < m ; j ++ ) {
w = grph[v][j].first, val = grph[v][j].second;
if( val <= A ) AddEq( su[w][2], Mul( coe, f[A / val] ) );
if( val <= C ) AddEq( su[w][0], Mul( coe, f[C / val] ) );
}
}
if( grph[u][i].second <= A ) {
int coe = Mul( mu[v], f[A / grph[u][i].second] );
for( int j = 0, val ; j < m ; j ++ ) {
w = grph[v][j].first, val = grph[v][j].second;
if( val <= C ) AddEq( su[w][1], Mul( coe, f[C / val] ) );
if( 1ll * v * w <= B ) AddEq( su[w][2], Mul( coe, f[B / v / w] ) );
}
}
if( grph[u][i].second <= C ) {
int coe = Mul( mu[v], f[C / grph[u][i].second] );
for( int j = 0, val ; j < m ; j ++ ) {
w = grph[v][j].first, val = grph[v][j].second;
if( val <= A ) AddEq( su[w][1], Mul( coe, f[A / val] ) );
if( 1ll * v * w <= B ) AddEq( su[w][0], Mul( coe, f[B / v / w] ) );
}
}
}
for( int i = 0, val ; i < n ; i ++ ) {
v = grph[u][i].first, val = grph[u][i].second;
int coe = Mul( mu[u], mu[v] );
if( 1ll * u * v <= B ) AddEq( ans, Mul( Mul( coe, su[v][1] ), f[B / u / v] ) );
if( val <= A ) AddEq( ans, Mul( Mul( coe, su[v][0] ), f[A / val] ) );
if( val <= C ) AddEq( ans, Mul( Mul( coe, su[v][2] ), f[C / val] ) );
}
}
Write( ans ), putchar( '\n' );
}
return 0;
}
