【NOIP2019模拟2019.10.07】果实摘取 (约瑟夫环、Mobius反演、类欧、Stern-Brocot Tree)

Description:


小 D 的家门口有一片果树林,果树上果实成熟了,小 D 想要摘下它们。
为了便于描述问题,我们假设小 D 的家在二维平面上的 (0, 0) 点,所有坐标范围的绝对值不超过 N 的整点坐标上都种着一棵果树。((0, 0) 这个点没有果树)
小 D 先站在 (0, 0) 处,正对着 (1, 0) 的方向。
每次摘果实时,小 D 会逆时针选择他能看到的第 K 棵还未摘取果实的果树,然后向着这个方向走去,在行走的过程中摘下沿路的所有的果树上的果树果实,直到走到果树林的边缘。
接下来,小 D 回到 (0, 0) 处,正对着上一次摘果实的果树的方向。
小 D 会重复这个过程,直到所有的果实都被摘取,小 D 感兴趣的是,最后一棵被摘下果实的果树是哪一棵?
注意小 D 不能看到被任何其他果树遮挡着的果树。

1 ≤ N, K ≤ 10^5

题解:


考虑总共的直线个数是(8sum_{i=1}^nphi(i))

第一个问题是(n=8sum_{i=1}^nphi(i))的环,每次取当前数的第k个,最后剩哪个的约瑟夫问题。

约瑟夫递推式(编号0-n-1):
(f[1]=0,f[n]=(f[n-1]+k)~mod~n)

推导显然很简单:
假设从0开始,第一次取了k-1,剩下(k,k+1,…,n-1,0,1,…k-2),以这个做一个(n-1)的子问题,得到的答案加k模n就是了。

由于(n很大),这个递推式显然会超时。

当n远大于k时,很久才会模一次,不妨一次跳多步,假设要x步,相当于下面的不等式:
(f[n]+k*x>=n+x)

(x=lceil{n-f[n] over k - 1} ceil)

估计复杂度,当(n<=k)时,(x=1)

(n>k)时,每次要加上(n/k)左右,也就是乘上({k+1over k})

所以大概是(O(k+log_{k+1over k}~n)≈O(k ~ln ~n))

第2个问题是求第x条直线,先全部转成第一象限。

这个比较套路,在Stern-Brocot Tree上二分,然后变成数数问题,

(cnt=sum_{i=1}^nsum_{j=1}^n[(i,j)=1]*[{jover i}<{b over a}])

(=sum_{d=1}^nmu(d)sum_{i=1}^{lfloor n/d floor}min({lfloor(b*i)/a floor},lfloor n/d floor))

预处理(mu)的前缀和,每次查询分块+类欧即可。

总复杂度:

(O(k~ln~n+n(n^{2/3})+sqrt n *log^3n))

Code:


#include<bits/stdc++.h>
#define fo(i, x, y) for(int i = x, B = y; i <= B; i ++)
#define ff(i, x, y) for(int i = x, B = y; i <  B; i ++)
#define fd(i, x, y) for(int i = x, B = y; i >= B; i --)
#define ll long long
#define pp printf
#define hh pp("
")
using namespace std;

const int N = 1e5 + 5;

int n, k;

int bz[N], p[N], p0, mu[N], phi[N], smu[N];

void sieve(int n) {
	phi[1] = mu[1] = 1;
	fo(i, 2, n) {
		if(!bz[i]) p[++ p0] = i, phi[i] = i - 1, mu[i] = -1;
		for(int j = 1; i * p[j] <= n; j ++) {
			int k = i * p[j]; bz[k] = 1;
			if(i % p[j] == 0) {
				phi[k] = phi[i] * p[j];
				mu[k] = 0;
				break;
			}
			phi[k] = phi[i] * phi[p[j]];
			mu[k] = -mu[i];
		}
	}
	fo(i, 1, n) smu[i] = smu[i - 1] + mu[i];
}

ll divs(ll x, ll y) {
	return x / y + (x % y != 0);
}
ll calc(ll n, ll m) {
	if(m == 1) return n;
	ll x = 1, s = 0;
	while(x < n) {
		ll k = divs(x - s, m - 1);
		if(x + k > n) return s + (n - x) * m + 1;
		x += k; s = (s + m * k) % x;
	}
	return s + 1;
}

ll calc(ll n, ll a, ll b, ll c) {
	if(n < 0) return 0;
	if(a >= c || b >= c) return calc(n, a % c, b % c, c) + n * (n + 1) / 2 * (a / c) + (n + 1) * (b / c);
	ll m = (a * n + b) / c;
	return n * m - calc(m - 1, c, c - b - 1, a);
}
ll ca(int n, int a, int b) {
	ll m = (ll) b * n / a;
	if(n <= m) { return calc(n, a, 0, b);}
	return calc(m, a, 0, b) + (ll) (n - m) * n;
}
ll count(int b, int a) {
	ll ans = 0;
	for(int i = 1, j; i <= n; i = j + 1) {
		j = n / (n / i);
		ans += ca(n / i, b, a) * (smu[j] - smu[i - 1]);
	}
	return ans;
}

ll sp;
int ax, ay;

void work(ll t) {
	int x1 = 0, y1 = 1, x2 = 1, y2 = 0;
	while(1) {
		int x3 = x1 + x2, y3 = y1 + y2;
		if(max(x3, y3) > n) break;
		if(count(x3, y3) <= t) {
			int c = 1;
			while(1) {
				int x4 = x3 + c * x2, y4 = y3 + c * y2;
				if(x4 > n || y4 > n || count(x4, y4) > t) break;
				c *= 2;
			}
			for(; c; c /= 2) {
				int x4 = x3 + c * x2, y4 = y3 + c * y2;
				if(x4 <= n && y4 <= n && count(x4, y4) <= t) x3 = x4, y3 = y4;
			}
			x1 = x3, y1 = y3;
			ax = x3, ay = y3;
		} else {
			int c = 1;
			while(1) {
				int x4 = x3 + c * x1, y4 = y3 + c * y1;
				if(x4 > n || y4 > n || count(x4, y4) <= t) break;
				c *= 2;
			}
			for(; c; c /= 2) {
				int x4 = x3 + c * x1, y4 = y3 + c * y1;
				if(x4 <= n && y4 <= n && count(x4, y4) > t) x3 = x4, y3 = y4;
			}
			x2 = x3, y2 = y3;
		}
	}
	int c = min(n / ax, n / ay);
	ax *= c, ay *= c;
	swap(ax, ay);
}

int main() {
	freopen("garden.in", "r", stdin);
	freopen("garden.out", "w", stdout);
	sieve(1e5);
	scanf("%d %d", &n, &k);
	fo(i, 1, n) sp += 2 * phi[i];
	ll t = calc(sp * 4, k);
	if(t <= sp) {
		if(t == 1) ax = n, ay = 0; else
		work(t - 1);
	} else
	if(t <= 2 * sp) {
		if(t == sp + 1) ax = 0, ay = n; else
		work(2 * sp - t + 1), ax = -ax;
	} else
	if(t <= 3 * sp) {
		if(t == 2 * sp + 1) ax = -n, ay = 0; else
		work(t - 2 * sp - 1), ax = -ax, ay = -ay;
	} else {
		if(t == 3 * sp + 1) ax = 0, ay = -n; else
		work(4 * sp - t + 1), ay = -ay;
	}
	pp("%d %d
", ax, ay);
}
原文地址:https://www.cnblogs.com/coldchair/p/11630565.html