「LOJ6482」LJJ爱数数

「LOJ6482」LJJ爱数数

解题思路 :

打表发现两个数 \(a, b\) 合法的充要条件是(我不管，我就是打表过的)：

\[a + b = \text{gcd}(a, b)^2 \]

设 \(g = \text{gcd(a, b)}\) ，那么相当于是要求：

\[\sum_{g=1}^{\sqrt{2n}}\sum_{i}[\text{gcd}(g^2-ig, ig)=g] \]

化简一波：

\[\sum_{g=1}^{\sqrt{2n}}\sum_{i}[\text{gcd}(g-i, i)=1] \]

根据辗转相除：

\[\sum_{g=1}^{\sqrt{2n}}\sum_{i}[\text{gcd}(g, i)=1] \]

考虑 \(i\) 的上界和下界

\[1 \leq ig \leq n \\ 1 \leq g^2 -ig \leq n \]

解一下这两个不等式：

\[\text{max}_i =\min(\lfloor\frac{n}{g}\rfloor,g - 1) \\ \text{min}_i =\max(g-\lfloor\frac{n}{g}\rfloor,1) \]

原来的式子相当于求：

\[\sum_{g=1}^{\sqrt{2n}}\sum_{i=\text{min}}^{\max}[\text{gcd}(g, i)=1] \]

设 \(f(n)\) 表示 \([1, n]\) 之间与 \(g\) 互质的数的个数，反演一波可以得到：

\[f(n)= \sum_{d|g} \lfloor \frac{n}{d}\rfloor \mu(d) \]

再化简一波式子：

\[\sum_{g=1}^{\sqrt{2n}}f(\max) -f(\min-1) \]

总复杂度 \(O(\sqrt{n}logn)\) 。

code

/*program by mangoyang*/
#pragma GCC optimize("Ofast", "inline")
#include<bits/stdc++.h>
#define inf ((int)(1e9))
#define Max(a, b) ((a) > (b) ? (a) : (b))
#define Min(a, b) ((a) < (b) ? (a) : (b))
typedef long long ll;
using namespace std;
template <class T>
inline void read(T &x){
    int f = 0, ch = 0; x = 0;
    for(; !isdigit(ch); ch = getchar()) if(ch == '-') f = 1;
    for(; isdigit(ch); ch = getchar()) x = x * 10 + ch - 48;
    if(f) x = -x;
}
const int N = 1500005;
int prime[N], d[30000005], tmp[N], mu[N], len[N], b[N], tot;
ll n, ans, m;
inline int solve(int x, int n){
	int ans = 0;
	for(register int i = len[x-1] + 1; i <= len[x]; i++) ans += mu[d[i]] * (n / d[i]);
	return ans;
}
int main(){
	read(n), mu[1] = 1;
	for(int i = 2; i < N; i++){
		if(!b[i]) prime[++tot] = i, mu[i] = -1;
		for(int j = 1; j <= tot && i * prime[j] < N; j++){
			b[i*prime[j]] = 1;
			if(i % prime[j] == 0){ mu[i*prime[j]] = 0; break; }
			mu[i*prime[j]] = -mu[i];
		}
	}
	m = (int) sqrt(2ll * n);
	for(int i = 1; i <= m; i++) if(mu[i])
		for(int j = i; j <= m; j += i) len[j]++;
	for(int i = 1; i <= m; i++) len[i] += len[i-1];
	for(int i = 1; i <= m; i++) if(mu[i])
		for(int j = i; j <= m; j += i) d[(++tmp[j])+len[j-1]] = i;
	for(int g = 1; g <= m; g++)
		ans += solve(g, Min(n / g, g - 1)) - solve(g, Max(1, g - n / g) - 1);
	cout << ans << endl;
	return 0;
}