题目
点这里看题目。
分析
也许是相对来说不那么难骗分的一道 NOI T3。
首先我们需要知道分数 (frac{x}{y}) 合法的充要条件。这一步在考场上可以通过如下步骤得到：
根据
小学奥数可以知道，十进制下所有有限小数都满足分母为 (2^a5^b,a,binmathbb Z)；所有混循环小数都有 (frac{t}{10^x-10^y},x,yin mathbb N_+,x>y) 的形式。从而可以猜想 10 进制下条件为 ([(x,y)=1][(y,10)=1])。
打开计算器，或者写一些程序验证 10 进制下的猜想。
推广到 (k) 进制，应该有 ([(x,y)=1][(y,k)=1])。
打开计算器，在 2 进制、8 进制和 16 进制下分别验证猜想，最终敢于下结论。
这个命题的
不太严格证明如下：考虑任意一个 (k) 进制下的纯循环小数，忽略整数位之后，一个纯循环小数可以用一个长度为循环节的序列 ((a_1,a_2,dots,a_n)) 唯一确定：
[frac x y=0.overset{.}{a_1}a_2dots a_{n-1}overset{.}{a_n} ]又有：
[frac{xk^{n}}{y}=a_1a_2dots a_n.overset{.}{a_1}a_2dots a_{n-1}overset{.}{a_n}\ ]注意到此时 (lbrace frac x y brace=lbrace frac{xk^n}{y} brace)，所以有：
[egin{aligned} frac{x}{y}-lfloorfrac x y floor&=frac{xk^n}{y}-lfloorfrac{xk^n}{y} floor\ x-ylfloor frac x y floor&=xk^n-ylfloorfrac {xk^n}{y} floor\ k^nx&equiv xpmod y end{aligned} ]由于 ((x,y)=1)，所以 (x) 存在逆元：
[egin{aligned} k^n&equiv 1pmod yLeftrightarrow (y,k)=1 end{aligned} ]
顺便还可以知道 (operatorname{ord}_y(k)|n)。
因而问题就愉快地变成了：
[sum_{x=1}^nsum_{y=1}^m[(x,y)=1][(y,k)=1]
]
此时我们有多种方法可供选择：
方法 1
该方法目的在于处理掉 ([(x,y)=1])，那么就可以开始推式子了：
[egin{aligned}
&sum_{x=1}^nsum_{y=1}^m[(x,y)=1][(y,k)=1]\
=&sum_{d=1}^{min{n,m}}mu(d)sum_{d|x,xle n}sum_{d|y,yle m}[(y,k)=1]\
=&sum_{d=1}^{min{n,m}}mu(d)lfloorfrac{n}{d}
floorsum_{y'=1}^{lfloorfrac m d
floor}[(dy',k)=1]\
=&sum_{d=1}^{min{n,m}}mu(d)[(d,k)=1]lfloorfrac{n}{d}
floorsum_{y'=1}^{lfloorfrac m d
floor}[(y',k)=1]\
end{aligned}
]
现在我们可以想办法求出下面两个东西，然后就可以整除分块了：
[egin{aligned}
f(n)&=sum_{j=1}^{n}[(j,k)=1]\
S(n,k&)=sum_{j=1}^{n}mu(j)[(j,k)=1]
end{aligned}
]
对于 (f)，注意到对于 (0le r<k)，如果有 ((r,k)=1)，那么对于 (qge 0)，也会有 ((r+qk,k)=1)，因而得到：
[f(n)=lfloorfrac n k
floor f(k)+f(nmod k)
]
然后考虑 (S(n,k))，我们可以继续使用 Mobius 反演：
[egin{aligned}
S(n,k)
&=sum_{j=1}^nmu(j)[(j,k)=1]\
&=sum_{d|k}mu(d)sum_{d|j,jle n}mu(j)\
&=sum_{d|k}mu(d)sum_{j'=1}^{lfloorfrac{n}{d}
floor}mu(dj')\
&=sum_{d|k}mu(d)sum_{j'=1}^{lfloorfrac n d
floor}mu(dj')[(d,j')=1]\
&=sum_{d|k}mu^2(d)S(lfloorfrac n d
floor,d)
end{aligned}
]
边界情况为 (n=0) 或者 (k=1)，后者可以使用杜教筛计算 (mu) 的前缀和。
这里 (S) 的总状态数只有 (O(sigma_0(k)sqrt{n}))，甚至可以直接递推。
小结：
- 这里关于 (f) 的处理比较重要，因为如果直接将 (f) 用 Mobius 反演展开反而会得到一个复杂度较高的算法，导致 GG；而此处不仅简化了运算，还揭示了数列 ({a_n=(n,k)}) 的循环特性！
- 处理 (S(n,k)) 过程中，运用到了 (mu) 取值的特性和积性，值得注意。
方法 2
该方法目的在于处理掉 ([(y,k)]=1)，接着我们又可以开始推式子了：
[egin{aligned}
&sum_{x=1}^nsum_{y=1}^m[(x,y)=1][(y,k)=1]\
=&sum_{x=1}^nsum_{y=1}^m[(x,y)=1]sum_{d|y,d|k}mu(d)\
=&sum_{d|k}mu(d)sum_{x=1}^nsum_{y'=1}^{lfloorfrac m d
floor}[(x,y'd)=1]\
=&sum_{d|k}mu(d)sum_{x=1}^nsum_{y'=1}^{lfloorfrac m d
floor}[(x,y')=1][(x,d)=1]\
end{aligned}
]
此时新的式子里已经出现了相似的结构，我们可以设 (f(n,m,k)=sum_{x=1}^nsum_{y=1}^m[(x,y)=1][(y,k)=1])，那么就有：
[f(n,m,k)=sum_{d|k}mu(d)f(lfloorfrac m d
floor,n,d)
]
边界情况是 (k=1)，直接整除分块配合杜教筛食用。
小结：
-
抓住了问题的相似结构，接着就变成了子问题。这样的思路其实更为常用。
废话，这不就是分治吗，只不过是对数进行分治
代码
方法 1
#include <cmath>
#include <cstdio>
#include <vector>
#include <algorithm>
#include <unordered_map>
using namespace std;
#define rep( i, a, b ) for( int i = (a) ; i <= (b) ; i ++ )
#define per( i, a, b ) for( int i = (a) ; i >= (b) ; i -- )
typedef long long LL;
const int MAXN = 1e6 + 5, MAXK = 2e3 + 5;
//const int MAXN = 10 + 5;
//const int MAXS = 10 + 5, MAXK = 10 + 5;
template<typename _T>
void read( _T &x )
{
x = 0; char s = getchar(); int f = 1;
while( s < '0' || '9' < s ) { f = 1; if( s == '-' ) f = -1; s = getchar(); }
while( '0' <= s && s <= '9' ) { x = ( x << 3 ) + ( x << 1 ) + ( s - '0' ), s = getchar(); }
x *= f;
}
template<typename _T>
void write( _T x )
{
if( x < 0 ) putchar( '-' ), x = -x;
if( 9 < x ) write( x / 10 );
putchar( x % 10 + '0' );
}
template<typename _T>
_T MAX( const _T a, const _T b )
{
return a > b ? a : b;
}
template<typename _T>
_T MIN( const _T a, const _T b )
{
return a < b ? a : b;
}
vector<int> fact[MAXK];
unordered_map<int, LL> mp, memS[MAXK];
LL f[MAXK];
int mu[MAXN], pref[MAXN];
int prime[MAXN], pn;
bool isPrime[MAXN];
int N, M, K, rt, L, U;
inline LL F( const int n ) { return 1ll * ( n / K ) * f[K] + f[n % K]; }
inline int Gcd( int a, int b ) { for( int z ; b ; z = a, a = b, b = z % b ); return a; }
void EulerSieve( const int n )
{
mu[1] = 1;
for( int i = 2 ; i <= n ; i ++ )
{
if( ! isPrime[i] ) prime[++ pn] = i, mu[i] = -1;
for( int j = 1 ; j <= pn && 1ll * i * prime[j] <= n ; j ++ )
{
isPrime[i * prime[j]] = true;
if( ! ( i % prime[j] ) ) break;
mu[i * prime[j]] = - mu[i];
}
}
rep( i, 1, n ) pref[i] = pref[i - 1] + mu[i];
}
LL Smu( const int n )
{
if( n <= 1e6 ) return pref[n];
if( mp.find( n ) != mp.end() ) return mp[n];
LL ret = 1;
for( int l = 2, r ; l <= n ; l = r + 1 )
{
r = n / ( n / l );
ret -= 1ll * ( r - l + 1 ) * Smu( n / l );
}
return mp[n] = ret;
}
LL S( const int n, const int k )
{
if( n <= 0 ) return 0;
if( k == 1 ) return Smu( n );
if( memS[k].find( n ) != memS[k].end() ) return memS[k][n];
LL ret = 0;
rep( i, 0, ( int ) fact[k].size() - 1 )
ret += 1ll * mu[fact[k][i]] * mu[fact[k][i]] * S( n / fact[k][i], fact[k][i] );
return memS[k][n] = ret;
}
int main()
{
// freopen( "cyclic.in", "r", stdin );
// freopen( "cyclic.out", "w", stdout );
read( N ), read( M ), read( K );
L = MIN( N, M ), U = MAX( N, M );
EulerSieve( 1e6 ), rt = sqrt( U );
rep( i, 1, K ) f[i] = f[i - 1] + ( Gcd( i, K ) == 1 );
for( int d = 1 ; d <= K ; d ++ )
for( int j = d ; j <= K ; j += d )
fact[j].push_back( d );
LL ans = 0;
for( int l = 1, r ; l <= L && l <= N && l <= M ; l = r + 1 )
{
r = MIN( MIN( N / ( N / l ), M / ( M / l ) ), L );
ans += 1ll * ( N / l ) * ( S( r, K ) - S( l - 1, K ) ) * F( M / l );
}
write( ans ), putchar( '
' );
return 0;
}
方法 2
你见过鸽子吗