(mathcal{Description})
Link.
给定 (n,m,k)，求 (xin [1,n]capmathbb N,yin [1,m]cap mathbb N)，且最简分数 (frac{x}{y}) 在 (k) 进制下是纯循环小数（包括整数）的 ((x,y)) 数量。
(n,mle10^9)，(kle2 imes10^3)。
(mathcal{Solution})
当你举几个十进制的纯循环小数就不难发现规律了。（
考虑一个已有 (xperp y) 的 (frac{x}y)，假设它是 (k) 进制下的纯循环小数，且循环节长度为 (l)。记 ({x}) 表示 (x) 的小数部分值，那么有
[left{frac{xk^l}{y}
ight}=left{frac{x}{y}
ight}\Leftrightarrow~~~~frac{xk^l}{y}-leftlfloorfrac{xk^l}{y}
ight
floor=frac{x}{y}-leftlfloorfrac{x}{y}
ight
floor\Leftrightarrow~~~~xk^l-yleftlfloorfrac{xk^l}{y}
ight
floor=x-yleftlfloorfrac{x}{y}
ight
floor\Leftrightarrow~~~~xk^lequiv xpmod y\Leftrightarrow~~~~kperp y
]
所以题目就是要求
[sum_{i=1}^msum_{j=1}^n[iperp j][iperp k]
]
略微推一下式子嘛：
[egin{aligned}sum_{i=1}^msum_{j=1}^n[iperp j][iperp k]&=sum_{i=1}^m[iperp k]sum_{j=1}^nsum_{dmid i,dmid j}mu(d)\&=sum_{d=1}^{min{n,m}}[dperp k]mu(d)sum_{i=1}^{lfloorfrac{m}{d}
floor}[iperp k]lfloorfrac{n}{d}
floor\&=sum_{d=1}^{min{n,m}}[dperp k]mu(d)lfloorfrac{n}{d}
floorsum_{i=1}^{lfloorfrac{m}{d}
floor}[iperp k]end{aligned}
]
套上整除分块，分别研究两个求和，令
[f(n)=sum_{i=1}^n[iperp k]\g(n,k)=sum_{i=1}^n[iperp k]mu(i)
]
快速解决它们，就能整除分块啦。
先考虑 (f)，显然的事实是 ([iperp k]=[(imod k)perp k])，继而有
[f(n)=lfloorfrac{n}{k}
floor f(n)+f(nmod k)
]
注意到 (k) 很小，(mathcal O(k)) 预处理之后就能 (mathcal O(1)) 求 (f) 了。
对于 (g) 而言，([iperp k]) 还能继续莫反——
[egin{aligned}g(n,k)&=sum_{i=1}^n[iperp k]mu(i)\&=sum_{i=1}^nmu(i)sum_{dmid i,dmid k}mu(k)\&=sum_{dmid k}mu(d)sum_{i=1}^{lfloorfrac{n}{d}
floor}mu(id)\&=sum_{dmid k}mu(d)sum_{i=1}^{lfloorfrac{n}{d}
floor}[iperp d]mu(id)~~~~*\&=sum_{dmid k}(mu(d))^2g(lfloorfrac{n}{d}
floor,d)end{aligned}
]
其中，标注 (*) 的步骤同时利用 (mu) 自身和积性函数普遍的性质进行“无用”转化，巧妙地完成了递推式。直接记忆化计算上式 就可以在可观的复杂度内求出 (g) 了，特别地，当 (k=1)，需要用杜教筛求 (mu) 的前缀和。
复杂度据说是 (mathcal O(sigma_0(k)n^{frac{1}2}+n^{frac{2}3}))。
(mathcal{Code})
/* Clearink */
#include <cmath>
#include <cstdio>
#include <unordered_map>
#define rep( i, l, r ) for ( int i = l, rep##i = r; i <= rep##i; ++i )
#define per( i, r, l ) for ( int i = r, per##i = l; i >= per##i; --i )
typedef long long LL;
const int MAXK = 2e3, MAXS = 1e7;
int n, m, K, pn, pr[MAXS + 5];
bool vis[MAXS + 5];
int f[MAXK + 5], mu[MAXS + 5], mus[MAXS + 5];
inline int imin( const int a, const int b ) { return a < b ? a : b; }
inline int gcd( const int a, const int b ) { return b ? gcd( b, a % b ) : a; }
inline LL calcF( const int n ) {
return ( n / K ) * f[K] + f[n % K];
}
inline void sieve() {
mu[1] = mus[1] = 1;
rep ( i, 2, MAXS ) {
if ( !vis[i] ) mu[pr[++pn] = i] = -1;
for ( int j = 1, t; j <= pn && ( t = i * pr[j] ) <= MAXS; ++j ) {
vis[t] = true;
if ( !( i % pr[j] ) ) break;
mu[t] = -mu[i];
}
mus[i] = mu[i] + mus[i - 1];
}
}
inline int calcM( const int n ) {
static std::unordered_map<int, int> mem;
if ( n <= MAXS ) return mus[n];
if ( mem.count( n ) ) return mem[n];
int ret = 1;
for ( int l = 2, r; l <= n; l = r + 1 ) {
r = n / ( n / l );
ret -= ( r - l + 1 ) * calcM( n / l );
}
return mem[n] = ret;
}
inline LL calcS( const int n, const int k ) {
static std::unordered_map<LL, LL> mem;
if ( !n ) return 0;
if ( k == 1 ) return calcM( n );
LL h = n * 2012ll + k;
if ( mem.count( h ) ) return mem[h];
LL ret = 0;
rep ( i, 1, sqrt( 1. * k ) ) if ( !( k % i ) ) {
ret += mu[i] * mu[i] * calcS( n / i, i );
if ( i * i != k ) {
ret += mu[k / i] * mu[k / i] * calcS( n / ( k / i ), k / i );
}
}
return mem[h] = ret;
}
int main() {
// freopen( "cyclic.in", "r", stdin );
// freopen( "cyclic.out", "w", stdout );
scanf( "%d %d %d", &n, &m, &K );
sieve();
rep ( i, 1, K ) f[i] = f[i - 1] + ( gcd( i, K ) == 1 );
LL ans = 0;
for ( int l = 1, r, t = imin( n, m ); l <= t; l = r + 1 ) {
r = imin( n / ( n / l ), m / ( m / l ) );
ans += ( calcS( r, K ) - calcS( l - 1, K ) )
* ( n / l ) * calcF( m / l );
}
printf( "%lld
", ans );
return 0;
}