【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
题解：
考虑总共的直线个数是\(8\sum_{i=1}^n\phi(i)\)。
第一个问题是\(n=8\sum_{i=1}^n\phi(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}\rceil\)
估计复杂度，当\(n<=k\)时，\(x=1\)
当\(n>k\)时，每次要加上\(n/k\)左右，也就是乘上\({k+1\over k}\)。
所以大概是\(O(k+log_{k+1\over k}~n)≈O(k ~ln ~n)\)。
第2个问题是求第x条直线，先全部转成第一象限。
这个比较套路，在Stern-Brocot Tree上二分，然后变成数数问题，
\(cnt=\sum_{i=1}^n\sum_{j=1}^n[(i,j)=1]*[{j\over i}<{b \over a}]\)
\(=\sum_{d=1}^n\mu(d)\sum_{i=1}^{\lfloor n/d \rfloor}min({\lfloor(b*i)/a\rfloor},\lfloor n/d \rfloor)\)
预处理\(\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("\n")
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\n", ax, ay);
}
