[NOI2016] 循环之美

题目

点这里看题目。

分析

也许是相对来说不那么难骗分的一道 NOI T3。

首先我们需要知道分数 (frac{x}{y}) 合法的充要条件。这一步在考场上可以通过如下步骤得到：

1. 根据小学奥数可以知道，十进制下所有有限小数都满足分母为 (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])。

2. 打开计算器，或者写一些程序验证 10 进制下的猜想。

3. 推广到 (k) 进制，应该有 ([(x,y)=1][(y,k)=1])。

4. 打开计算器，在 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}))，甚至可以直接递推。

小结：

1. 这里关于 (f) 的处理比较重要，因为如果直接将 (f) 用 Mobius 反演展开反而会得到一个复杂度较高的算法，导致 GG；而此处不仅简化了运算，还揭示了数列 ({a_n=(n,k)}) 的循环特性！
2. 处理 (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. 抓住了问题的相似结构，接着就变成了子问题。这样的思路其实更为常用。

   废话，这不就是分治吗，只不过是对数进行分治

代码

方法 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

你见过鸽子吗