【YBT2023寒假Day2 C】网格与圆（莫比乌斯反演）（斐蜀定理）

网格与圆

题目链接：YBT2023寒假Day2 C

题目大意

一个长宽都是 n 的网格，每个格子长宽是 1，除了最左下角，每个格子里有一个半径为 R 的圆，圆心在格子正中心。
然后问你站在最左下角格子正中心，问你能看到多少个圆。
（1e-6<=R<=0.5）

思路

首先如果圆看做点，那不就是仪仗队，所以我们思考可能是在哪个算法上改进一下。

考虑找更多的性质，与圆的半径有关的。

---

首先自己尝试一下会发现可能一个圆被挡住的话与它的圆心被挡住是充要的。
首先如果圆心没有被挡住一定能看到，那我们要证明的就是被挡住的话整个圆都会被挡住。

考虑挡住 \((u,v)\) 的是 \((i,j)\)，设 \((i,j)\) 在直线 \(uy=vx\) 的下面。
然后因为 \((i,j)\) 比 \((u,v)\) 距离观察点近，那它下切线跟观察点连的线跟 \(uy=vx\) 的夹角就更大。
那也就是 \((u,v)\) 直线下面的部分都会被挡住。
同理折过去，\((u-i,v-j)\) 也会挡住它，而且会在上面，那上面的部分也会被挡住，那整个圆就被挡住了。

---

那我们就只用看圆心，通过数学的公式可以做到 \(O(n^2)\)，考虑快速判断一个点是否被挡住。
考虑观察点与直线的距离公式：
\(\dfrac{|Ax+By+C|}{\sqrt{A^2+B^2}}\)
带入 \((u,v)\) 对于的直线 \(uy=vx\)，以及所有 \(0\leqslant i\leqslant u,0\leqslant j\leqslant v,(i,j)\neq(u,v)\) 的 \((i,j)\)，要有所有的距离都要大于 \(R\)。
那带入进去：
\(\dfrac{|iv-ju|}{\sqrt{u^2+v^2}}>R\)
\(|iv-ju|>R\sqrt{u^2+v^2}\)

然后注意到我们利用上仪仗队的性质，那就是要 \(\gcd(u,v)=1\)。
那左边这个我们发现可以利用上 \(\gcd(u,v)=1\)，通过斐蜀定理，它的最小值就是 \(1\)。
那就是 \(1>R\sqrt{u^2+v^2}\)
\(u^2+v^2<\dfrac{1}{R^2}\)

---

那问题就转化了，变成：
有多少对 \((x,y)\in\{0,1,...,n-1\}^2\) 使得 \(\gcd(x,y)=1,x^2+y^2<\dfrac{1}{R^2}\)。
那会发现这就是一个经典的莫反问题，那前面 \(\gcd\) 的我们可以枚举 \(d\)（因为这里有了 \(R\) 的限制所以其实范围就是 \(\min(n-1,\dfrac{1}{R})\)），数里面 \(x^2+y^2<\dfrac{1}{R^2}\) 的个数然后乘上 \(\mu(d)\) 即可。

然后至于数个数你就枚举一个 \(x\)，然后通过开根求出 \(y\) 的最大范围。
不过要注意的一点是这里不能取等号，所以你必须注意除了是下取整，如果刚好是 \(x.00000\)，你的答案应该是 \(x-1\)。

代码

#include<cmath>
#include<cstdio>
#include<vector>
#define ll long long


using namespace std;

const int N = 1e6 + 100;
ll n, pr[N], mu[N];
bool np[N];
long double r, m;

void Init() {
	mu[1] = 1;
	for (int i = 2; i < N; i++) {
		if (!np[i]) {
			mu[i] = -1;
			pr[++pr[0]] = i;
		}
		for (int j = 1; j <= pr[0] && pr[j] * i < N; j++) {
			np[i * pr[j]] = 1;
			if (i % pr[j] != 0) mu[i * pr[j]] = -mu[i];
				else break;
		}
	}
}

ll slove(ll d) {
	ll re = 0; long double temp = m * m / (long double)d / (long double)d;
	for (ll i = 1; i * d <= n && i * d <= m; i++) {
		long double tmp = sqrt(temp - i * i);
		ll now = floor(tmp);
		if (now == tmp) now--;
		now = min(now, (ll)(n / d));
		re += now;
	}
	return re;
}

int main() {
	freopen("circle.in", "r", stdin);
	freopen("circle.out", "w", stdout); 
	
	Init();
	
	scanf("%lld %Lf", &n, &r); 
	
//	if (n == 48 && 200000 == r) while (1);
	
	n--; r = r / 1000000; m = 1.0 / r;
	
	if (!n) {
		printf("0"); return 0; 
	}
	
	ll ans = 0;
	for (ll d = 1; d <= m && d <= n; d++) {
		ans += mu[d] * slove(d);
	}
	printf("%lld", ans + 2);
	
	return 0;
}