[NOI2016]循环之美

题目

传送门 to luogu

思路

纯循环小数 $a$，设循环节长度为 $L$，那么 $k^{L}a$ 与 $a$ 的小数部分相同，作差有 $(k^L-1)a\in \mathbb{Z}$，所以 $a=\frac{v}{k^L-1}(v\in \mathbb{Z})$ 。另一方面，若 $a$ 满足 $(k^L-1)a\in \mathbb{Z}$，那么 $a$ 一定是纯循环小数，因为小数部分必然需要有长度为 $L$ 的循环节（虽然可能不是最短）。

所以 $a$ 是纯循环小数，当且仅当 $a=\frac{v}{k^L-1}(v\in Z)$ 。哪怕 $\gcd (v,k^L-1)\ne 1$ 也并不影响这一点。于是答案就是
$$\begin{array}{rl}ans& =\sum _{i=1}^n\sum _{j=1}^m\left[\exists L,x,\text{s.t.}\frac{i}{j}=\frac{x}{k^L-1}\right][\gcd (i,j)=1]\\ & =\sum _{i=1}^n\sum _{j=1}^m\left[\exists L,j\mid (k^L-1)\right][\gcd (i,j)=1]\end{array}$$

显然 $i$ 的作用不太大。于是我记 $f(j)={\sum }_{i=1}^n[i\perp j],g(j)=[\exists L,j\mid (k^L-1)]$，那么
$$ans=\sum _{j=1}^{m}f(j)g(j)$$

以前的我是这样做的：由于学的是 $min25$ 筛（只要是积性函数就能筛），所以直接考察函数是否为积性函数。比如 $g(x)$，它是积性函数吗？

我们不难发现 $d\mid x$ 则 $(k^d-1)\mid (k^x-1)$ 这个结论。证明可以利用因式定理，或者用神奇的
$$\gcd (k^a-1,k^b-1)=k^{\gcd (a,b)}-1$$

来很容易地说明。然后，以前的我就看出来了 $g(x)$ 是积性函数。证明很简单，若 $a\mid (k^u-1)$ 且 $b\mid (k^v-1)$，则 $ab\mid (k^{uv}-1)(a\perp b)$ 。而 $a\nmid (k^L-1)$ 也就不可能有 $ab\mid (k^{\nu }-1)$ 了。

于是考虑 $g(p)(p\in \text{Prime})$ 的值。只要 $k$ 不是 $p$ 的倍数，就一定有 $k^{\phi (p)}-1\equiv 0(\mathrm{mod}p)$，故 $g(p)=[\gcd (k,p)=1]$ 。

进一步，$k^{\phi (p^2)}-1\equiv 0(\mathrm{mod}{p}^2)$，只要 $k$ 不是 $p$ 的倍数。那么我们写出

$$g(p^x)=[k\perp p]$$

再利用 $g$ 是积性的，我们得到
$$g(x)=[k\perp x]$$

现在的我回看该式子：简直是显然啊！原式等价于 $k^L\equiv 1(\mathrm{mod}x)$，在 $k\perp x$ 时有解 $L=\phi (x)$，在 $\gcd (k,x)\ne 1$ 时必然无解。更可笑的是，现在的我连原来这种繁杂的证法也不会……

顺便说一句，这里的 $g(x)$ 不仅是积性函数了，而且是 完全积性 函数。

那么 $f(j)$ 怎么算捏？试一试：
$$\begin{array}{rl}f(j)& =\sum _{i=1}^n[\gcd (i,j)=1]\\ & =\sum _{i=1}^n\sum _{p|i,p|j}\mu (p)\\ & =\sum _{p|j}\mu (p)\lfloor \frac{n}{p}\rfloor \end{array}$$

所以我们重写一下开头的式子
$$\begin{array}{rl}ans& =\sum _{j=1}^m\sum _{p|j}\mu (p)\lfloor \frac{n}{p}\rfloor g(j)\\ & =\sum _{p=1}^m\lfloor \frac{n}{p}\rfloor \mu (p)\sum _{p|j,j⩽m}g(j)\\ & =\sum _{p=1}^m\lfloor \frac{n}{p}\rfloor \mu (p)g(p)\sum _{i=1}^{\lfloor \frac{m}{p}\rfloor }g(i)\end{array}$$

最后一步就用到了 $g(x)$ 是完全积性函数的性质。对 $\lfloor \frac{n}{p}\rfloor ,\lfloor \frac{m}{p}\rfloor$ 同时整除分块，对 $\mu (p)g(p)$ 和 $g(i)$ 用 $min25$ 求前缀和即可。

具体怎么求呢？对于 $\mu (p)g(p)$ 和 $g(i)$ 的处理方法都一样，以 $g(i)$ 为例：在 $min25$ 的 $\text{step one}$ 中，用 $f(x)=1$ 去拟合质数处取值，然后暴力枚举 $\mathcal{O}(\sqrt{n})$ 个取值点中的 $\mathcal{O}(\frac{k}{\ln k})$ 个质数，找到不与 $k$ 互质的，减去之。$\text{step two}$ 就挺简单了。

然后看了看题解，发现 $g(x)$ 的前缀和可以更好求，因为 $g(i)$ 的前缀和显然有长度为 $k$ 的 “循环节”。我们可以 $\mathcal{O}(k\log k)$ 暴力求出前 $k$ 项，而后可以直接算。

代码

#include <cstdio>
#include <iostream>
#include <cstring>
#include <vector>
#include <algorithm>
using namespace std;
typedef long long int_;
inline int readint(){
	int a = 0; char c = getchar(), f = 1;
	for(; c<'0'||c>'9'; c=getchar())
		if(c == '-') f = -f;
	for(; '0'<=c&&c<='9'; c=getchar())
		a = (a<<3)+(a<<1)+(c^48);
	return a*f;
}

const int MaxN = 80005;
bool isPrime[MaxN];
vector< int > primes;
int smu[MaxN]; // prefix sum of mu
vector< int > sumf; // real F = mu*g
void sieve(int k){
	memset(isPrime+2,1,MaxN-2);
	smu[1] = 1; // just don't forget it
	sumf.push_back(0); // to be bottom
	for(int i=2,len=0; i<MaxN; ++i){
		if(isPrime[i]){
			primes.push_back(i), ++ len;
			if(k%i == 0) // g(i) = 0
				sumf.push_back(0);
			else sumf.push_back(-1);
		}
		for(int j=0; j<len&&primes[j]<=(MaxN-1)/i; ++j){
			isPrime[i*primes[j]] = false;
			if(i%primes[j] == 0) break;
		}
	}
	for(int i=2,len=primes.size(); i<=len; ++i)
		sumf[i] += sumf[i-1];
}

struct MIN_25{
	int haxi[2][MaxN], w[MaxN];
# define index_(x) (((x) >= MaxN) ? haxi[1][n/(x)] : haxi[0][x])
	int g[MaxN]; // F = mu*g begin with F = 1
	void min_25(int n,int k){
		int id = 0; // allocate index
		for(int i=1; i<=n; i=n/(n/i)+1){
			index_(n/i) = ++ id;
			g[id] = (w[id] = n/i)-1;
		}
		/* Sieve number of primes */ ;
		for(int j=0; primes[j]<=n/primes[j]; ++j)
			for(int i=1; i<=id; ++i){
				if(primes[j] > w[i]/primes[j])
					break; // p_j <= sqrt(w)
				g[i] -= (g[index_(w[i]/primes[j])]-j);
			}
		/* Remove special primes */ ;
		for(int i=id,j=0,x=0; i; --i){
			for(; primes[j]<=k; ++j){
				if(primes[j] > w[i]) break;
				x += !(k%primes[j]); // calc
			}
			g[i] -= x; // remove them
		}
		# define FF(x) ((k%(x)) ? -1 : 0)
		/* Get all ans */ ;
		for(int i=1; i<=id; ++i)
			g[i] = -g[i]; // mu(p)*g(p) = -1
		int j = primes.size();
		for(--j; ~j; --j){ // from big to small
			if(primes[j] > n/primes[j])
				continue; // p_j <= sqrt(n)
			for(int i=1; i<=id; ++i){
				if(primes[j] > w[i]/primes[j])
					break; // good cut branch
				g[i] += FF(primes[j])*(
					g[index_(w[i]/primes[j])]
					- sumf[j+1]); // JiXing
			}
		}
	}
};
MIN_25 fn, fm;

# define sumG(x) (((x)/k)*sxy[k]+sxy[(x)%k])
int sxy[MaxN];
int main(){
	int n = readint(), m = readint();
	int k = readint();
	for(int i=1; i<=k; ++i){
		sxy[i] = sxy[i-1]; // prefix sum
		if(__gcd(i,k) == 1) // coprime
			++ sxy[i]; // add one
	}
	sieve(k); fn.min_25(n,k);
	fm.min_25(m,k); // hc.calc(k);
	int lst = 0, now; int_ ans = 0, fk;
	for(int l=1,r; l<=m&&l<=n; l=r+1){
		r = min(n/(n/l),m/(m/l));
		fk = 1ll*sumG(m/l)*(n/l);
		if(r == n/(n/l)) if(r >= MaxN)
			now = fn.g[fn.haxi[1][n/r]];
			else now = fn.g[fn.haxi[0][r]];
		else if(r >= MaxN)
			now = fm.g[fm.haxi[1][m/r]];
		else now = fm.g[fm.haxi[0][r]];
		ans += fk*(now+1-lst);
		lst = now+1; // mu(1)*g(1) = 1
	}
	printf("%lld\n",ans);
	return 0;
}

后记

去年的 $\mathsf{OneInDark}$ 声称，他花费了 $3$ 天才把这题搞出来。今年的 $\mathsf{OneInDark}$ 却连此题的第一步都推不出来……