题解：P16132 [ICPC 2018 NAIPC] Probe Droids

· · 题解

题意简述

炮台位于网格左下角，初始朝右。若当前方向还有机器人，就继续击毁最近的机器人；否则逆时针旋转到下一个有机器人的方向。

给定若干排名，求对应被击毁机器人的行列坐标。

解题思路

把炮台平移到原点。记 N=n-1、M=m-1，用 (x,y) 表示机器人相对炮台的列偏移和行偏移。

击毁顺序就是按斜率 y/x 从小到大排列；斜率相同时，按离原点的距离从近到远排列。同一射线上的近点消失后，远点会立即进入视线，炮台不会先旋转。

水平轴上的 M 个机器人最先被击毁，竖直轴上的 N 个机器人最后被击毁，可以直接回答。剩下仅考虑 1\le x\le M、1\le y\le N 的 NM 个点。

先解决排名统计。对于正有理数 a/b，记 C(a,b) 为斜率不超过它的内部点数。第 x 列可以选择的行偏移数为 \min(N,\lfloor ax/b\rfloor)。令 t=\min(M,\lfloor Nb/a\rfloor)，则：

C(a,b)=N(M-t)+\sum_{x=1}^t\left\lfloor\frac{ax}{b}\right\rfloor

前 t 列尚未超过行数上限，后面的列均贡献 N。这也正确处理斜率恰好等于某条射线的情况，因为统计使用的是「不超过」。

求和部分使用类欧几里得算法。定义：

F(n,m,a,b)=\sum_{i=0}^{n-1}\left\lfloor\frac{ai+b}{m}\right\rfloor

于是 C(a,b)=N(M-t)+F(t,b,a,a)。每次先取出 a/m 和 b/m 的整商，累加对应的等差数列贡献。余数满足 0\le a,b<m 后，令 Y=an+b；若 Y<m，剩余贡献为 0，否则把参数替换为 (\lfloor Y/m\rfloor,a,m,Y\bmod m)。这相当于交换格点计数的横纵方向，分母按欧几里得过程减小。

有了单调的计数函数，可以二分第 k 个点的斜率。但直接使用浮点数，不容易保证非常接近的两条射线不会混淆。

采用固定分母 K=2^{42}，在整数区间 [0,NK+1) 内，寻找最小的整数 z，使 C(z,K)\ge k。若目标射线的最简斜率为 a/b，则有：

\frac{z-1}{K}<\frac{a}{b}\le\frac{z}{K}

由于 b\le M，任意两条不同的最简斜率 a/b,c/d 都满足：

\begin{aligned} \left|\frac{a}{b}-\frac{c}{d}\right| & =\frac{|ad-bc|}{bd} \\ & \ge\frac{1}{M^2} \end{aligned}

本题中 M<10^6，所以选定的 K 满足 K>4M^2。目标斜率距离 z/K 不超过 1/K，而其他分母不超过 M 的不同分数，距离至少为 1/M^2-1/K>1/K。因此，目标斜率恰好是分母不超过 M、距离 z/K 最近的分数。

这里不需要再限制近似分数的分子。真正的目标分数已经满足分子不超过 N，而上面的距离比较排除了所有其他分母不超过 M 的分数。

用连分数求这个最佳近似。维护相邻两项渐近分数 p/q 和 r/s，初始为 0/1 和 1/0。设当前连分数整商为 h，下一项是 (p+hr)/(q+hs)。若分母不超过 M，就继续欧几里得展开。

当下一项分母将超过 M 时，取最大的允许系数：

h'=\left\lfloor\frac{M-q}{s}\right\rfloor

此时仅需比较 r/s 和 (p+h'r)/(q+h's)。它们位于 z/K 两侧，交叉乘积之差的绝对值为 1，且分母之和超过 M。

两个交叉乘积之差为 1 的分数之间，任意其他最简分数的分母都至少是两者分母之和。可将中间分数的分子、分母表示为这两组分子、分母的正整数线性组合，从而得到这个下界。因此，两者之间没有分母不超过 M 的第三个候选；在两者外侧的分数，也不会比对应一侧的端点更接近 z/K。比较这两个分数与目标的精确距离即可。若连分数已经完整结束且分母满足限制，则直接返回该分数。

代码中的 approx 实现这一过程。所有距离比较都通过交叉相乘完成，不使用浮点数。连分数相邻项的分子和分母互质，所以最终返回的 a/b 已经是最简分数。

恢复斜率后，这条射线上的机器人依次为 (b,a),(2b,2a),\dots，共有：

D=\min\left(\left\lfloor\frac{N}{a}\right\rfloor,\left\lfloor\frac{M}{b}\right\rfloor\right)

严格小于该斜率的内部点数是 C(a,b)-D，所以第 k 个内部点在这条射线上的序号为 k-C(a,b)+D。将最简方向向量乘以这个序号，再把行列偏移各加 1，即可得到原坐标。

代码先处理两条坐标轴，因而后续的 a,b,N,M 都为正，不会出现除以 0。计数及排名使用 long long，类欧几里得中间乘积和近似距离比较使用 __int128。NK+1<2^{62}，二分区间也在有符号整数范围内。

每次询问至多进行 62 次整数二分判断，每次判断和连分数还原均需要对数时间。令 U=\max(n,m)，总时间复杂度为 O(q\log^2 U)，额外空间复杂度为 O(1)。

参考代码

#include <bits/stdc++.h>
using namespace std;

using ll=long long;
using i128=__int128;
const ll K=1LL<<42;
ll n,m;
ll floor_sum(ll n,ll m,ll a,ll b)
{
    ll res=0;
    while(1)
    {
        res+=(i128)(a/m)*n*(n-1)/2+(i128)(b/m)*n;
        a%=m;
        b%=m;
        i128 y=(i128)a*n+b;
        if(y<m)break;
        n=y/m;
        b=y%m;
        swap(a,m);
    }
    return res;
}
ll count(ll a,ll b)
{
    if(!a)return 0;
    ll t=min(m,(ll)((i128)n*b/a));
    return n*(m-t)+floor_sum(t,b,a,a);
}
pair<ll,ll> approx(ll a,ll b,ll m)
{
    ll p=0,q=1,r=1,s=0,x=a,y=b;
    while(y)
    {
        ll k=x/y;
        if(q+k*s>m)break;
        ll u=p+k*r,v=q+k*s;
        p=r;
        q=s;
        r=u;
        s=v;
        ll tmp=x%y;
        x=y;
        y=tmp;
    }
    if(!y)return {r,s};
    ll k=(m-q)/s;
    ll u=p+k*r,v=q+k*s;
    i128 l=(i128)a*v-(i128)b*u,h=(i128)a*s-(i128)b*r;
    if(l<0)l=-l;
    if(h<0)h=-h;
    return l*s<=h*v?pair<ll,ll>{u,v}:pair<ll,ll>{r,s};
}
int main()
{
    ios::sync_with_stdio(false);
    cin.tie(nullptr);
    int q;
    cin>>n>>m>>q;
    n--;
    m--;
    while(q--)
    {
        ll k;
        cin>>k;
        if(k<=m)
        {
            cout<<1<<' '<<k+1<<'\n';
            continue;
        }
        k-=m;
        if(k>n*m)
        {
            cout<<k-n*m+1<<' '<<1<<'\n';
            continue;
        }
        ll l=0,r=n*K+1;
        while(l<r)
        {
            ll mid=l+(r-l)/2;
            if(count(mid,K)>=k)r=mid;
            else l=mid+1;
        }
        auto [a,b]=approx(l,K,m);
        ll pos=k-count(a,b)+min(n/a,m/b);
        cout<<a*pos+1<<' '<<b*pos+1<<'\n';
    }
    return 0;
}