做 UVa 149 Forests 这题时最容易被卡住的不是算法而是你根本没搞懂它到底要数什么。这题看起来像一道几何题实际上是一道非常经典的格点可见性问题树被整整齐齐种在整数坐标点上你站在原点给你一个观测角度和一个最大距离让你数出这一块扇形范围里到底有多少棵树是真正“可见”的。我最早做这道题的时候第一反应是去算直线会不会穿过别的树结果又慢又容易错。后来才意识到只要把“可见”翻译成“坐标互质”整个题瞬间就变成了一道数论题。这篇文章我把完整的推导、参考实现和踩坑过程都整理出来适合刚开始刷 UVa 的读者也适合想复习 gcd、欧拉函数、容斥原理的选手。耐心看完你不仅能 AC 这道题还能把同一套思路迁移到不少类似的格点计数问题里。1. 这道题到底在问什么1.1 先还原一下森林场景题目里的森林不是自然生长的更像是整齐划一的果园。所有树都落在平面直角坐标系的整数格点上也就是说只要 x 和 y 都是整数坐标 (x, y) 那里就有一棵树。人站在原点 (0, 0)向森林里看去。输入会给你两个数一个是角度单位是度另一个是距离。意思是从正东方向开始逆时针扫过给定角度的扇形在这个扇形内部、并且到原点的距离不超过给定距离的范围里统计你能看到的树的数量。这里最关键的点在于不是所有在这个范围内的树都能被看到。如果一棵树正好被另一棵更靠近原点的树挡住那它就不算“可见”。比如 (2, 2) 这棵树它和 (1, 1) 在同一条射线上而且 (1, 1) 离你更近所以你只能看到 (1, 1)看不到 (2, 2)。这样的遮挡逻辑让问题从“数范围里的点”升级成了“数没有被遮挡的点”。如果把坐标系范围缩小一点比如给定距离是 10角度是 45 度那么符合条件的坐标点有一大片但真正能进入答案的只有那些在扇形范围内、又没有被其他树挡住的原始格点。你会发现这题本质上是在问在一个被角度和半径框出来的区域里有多少个“第一次出现的视向”。1.2 可见性判定为什么互质是核心要判断从原点到某个格点 (x, y) 的连线上有没有其他格点最简单的做法是把这条线段参数化。线段上的任何一点都可以写成 (t·x/g, t·y/g)其中 g gcd(|x|, |y|)。当 t 取整数时这个点就是整数格点。t 1 时点是 (x/g, y/g)也就是第一个经过的整数点t g 时点才是 (x, y) 本身也就是你最终要看的这棵树。如果 g 大于 1那么在 t 1 到 t g - 1 之间线段上还会经过其他树所以这棵树必然被挡住。只有当 g 1也就是 x 和 y 互质时线段上除了起点和目标点之外没有任何整数格点这棵树才会被看到。这个结论非常干净可见等价于互质。看到这个条件之后原先复杂的位置关系判断就变成了数论问题。后面所有的计数优化都建立在这一个结论的基础上。理解这一步比背代码重要得多因为不同题目的坐标系、角度范围可能完全不一样但只要抓住了 gcd 这个核心你怎么改都能游刃有余。2. 数学拆解把扇形离散化成格点条件2.1 角度约束与距离约束现在题目给的信息可以拆成两个离散条件。第一个是距离约束树到原点的欧氏距离不能超过给定值 D也就是 x² y² ≤ D²。第二个是角度约束假设从正东方向开始逆时针扫那么树的方位角必须在 0° 到给定角度 A 之间。对于第一象限的格点来说这个条件可以写成 y / x ≤ tan(A)也就是 y ≤ x·tan(A)。两个条件合在一起树的候选范围就是圆和扇形相交的那一块区域。但要注意这个区域里包含大量不可见点比如 (2,0) 和 (3,0) 这些正东方向的重复点它们都在 x 轴上却被最前面的 (1,0) 挡住。所以在实际操作中需要把距离、角度、互质三个条件同时判断。这个小节里所有公式都不难真正容易出错的是角度临界值。边界的树到底算不算可见按照题面通常的语义边界上的树是算的。比如 A 45 度时(1,1) 正好在扇形的边界线上它应该被计入答案。如果不做浮点精度处理tan(45°) 的理论值是 1.0但计算机算出来的可能是 0.9999999999999999 或者 1.0000000000000002直接转成整数就会漏掉这一棵树。我的习惯是引入一个非常小的 eps比如 1e-9在向下取整之前先加上这个 eps。这样既能把边界上的树保留下来又不会把明显超出范围的树误判进来。这个细节在 UVa 149 这种浮点数输入题里几乎是必踩的坑后面排查章节还会细说。2.2 坐标轴和边界上的特殊情况坐标轴上的情况是最容易被人忽略的。先看 x 轴也就是 y 0、x 0 的点。(1,0) 是可见的因为从原点到 (1,0) 的线段不经过其他整数点。但 (2,0)、(3,0) 这些点全部被 (1,0) 挡住所以不可见。这个规律用 gcd 来验证也完全一致gcd(1,0) 1gcd(2,0) 2gcd(3,0) 3。y 轴上的情况类似(0,1) 可见(0,2) 和 (0,3) 不可见。如果标准库的 gcd 函数能把 gcd(0, n) 返回成 n那么这个规律会被自动处理。但如果你自己写循环或者用某些不处理 0 参数的工具函数就要小心了。在按列扫描的实现里这些坐标轴点通常要单独拎出来处理。因为大多数高效算法只统计 y ≥ 1 的点再补上 (1,0) 这个 x 轴上的点。如果忘记补答案会少 1如果重复补答案会多 1。这种一多一少的问题恰恰是最难排查的。2.3 复杂度从 O(R²) 到 O(R·ω(x)) 的思路最直观的写法是枚举所有坐标点从 -R 到 R 两重循环判断距离、角度、gcd。距离 R 取整后复杂度是 O(R²)。如果输入距离只有几十、几百这个算法没有任何问题代码也最容易写对。但 UVa 这类题不会只给你小数据。距离变大之后O(R²) 会明显吃力。更高效的做法是只看 x 坐标一列一列地扫。对于固定的 xy 的可选范围可以靠距离和角度直接算出来然后我们要的其实是在这个 y 范围内有多少个 y 和当前 x 互质。统计一段区间里有多少个数和某个数互质这正是容斥原理的经典场景。把 x 做质因数分解然后按“能被某个质因子整除”来容斥。因为每个整数 x 的质因子种类很少比如 x ≤ 10^5 时最多有 6 个不同质因子所以容斥的枚举量是 2^kk 很小。整体复杂度可以降到 O(R·2^k) 左右比 O(R²) 快得多。3. 参考实现与验证3.1 按列枚举的核心函数先写一个很重要的子函数给定 x 和一个 y 的区间 [lo, hi]统计区间内有多少个 y 与 x 互质。这个函数是高效方案的心脏。它的做法是先对 x 做质因数分解得到互不相同的质因子列表。然后用二进制枚举的方式遍历所有非空子集计算这个子集对应质因子乘积的倍数个数再根据子集大小决定是加还是减。这是容斥原理的标准写法。vectorlong long getPrimeFactors(int n) { vectorlong long fac; for (long long p 2; p * p n; p) { if (n % p 0) { fac.push_back(p); while (n % p 0) n / p; } } if (n 1) fac.push_back(n); return fac; } long long countCoprime(int x, int lo, int hi) { if (lo hi) return 0; vectorlong long fac getPrimeFactors(x); int m (int)fac.size(); long long total hi - lo 1; for (int mask 1; mask (1 m); mask) { long long prod 1; int bits 0; for (int i 0; i m; i) { if (mask (1 i)) { prod * fac[i]; bits; } } long long cnt hi / prod - (lo - 1) / prod; if (bits 1) total - cnt; else total cnt; } return total; }这个函数本身不含任何复杂数据结构但很容易写错。容斥里的prod必须是 long long因为多个质因子乘起来之后可能超过 int。另外x 1 时质因子列表为空容斥循环不会执行函数直接返回区间长度这其实是正确的因为 1 和任何整数都互质。3.2 主计数函数一列一列算有了 countCoprime 之后主计数流程就很清晰了。对于每个 x先根据距离限制算出 y 的最大值 yMaxDist再根据角度限制算出 y 的最大值 yMaxAngle取两者较小值作为这一列的 y 上限。然后统计 y 从 1 到上限之间与 x 互质的数量记得把 x 轴上的点单独处理。具体来说只有 x 1 时(1,0) 才是可见的。其他 x 对应的 (x,0) 全部被 (1,0) 挡住。所以每当 x 1 时答案额外加 1。y 轴上的 (0,1) 也要单独统计条件是这个点必须在距离范围内并且角度范围包含 90 度。long long solveFirstQuadrant(double angle, double dist) { if (angle 0) angle 0; if (angle 90) angle 90; long long ans 0; int R (int)ceil(dist); double rad angle * acos(-1.0) / 180.0; double tanA tan(rad); for (int x 1; x R; x) { double rem dist * dist - 1.0 * x * x; if (rem -1e-9) continue; int yMaxDist (int)floor(sqrt(max(0.0, rem)) 1e-9); int yMaxAngle; if (angle 89.9999999) { yMaxAngle yMaxDist; } else { yMaxAngle (int)floor(1.0 * x * tanA 1e-9); } int yMax min(yMaxDist, yMaxAngle); if (yMax 0) continue; if (x 1) { ans; // 处理 (1,0) if (yMax 1) ans countCoprime(x, 1, yMax); } else { if (yMax 1) ans countCoprime(x, 1, yMax); // (x,0) 被 (1,0) 挡住不需要额外加 } } if (angle 90.0 - 1e-9 dist 1.0 - 1e-9) { ans; // 处理 (0,1) } return ans; }这里有一个值得注意的边界处理。当 angle 非常接近 90 度时tan(90°) 在计算机里会变成一个很大的数直接乘 x 再转 int 可能会溢出或者得到不可控的值。所以我在代码里做了一个特殊判断角度达到 89.9999999 以上时直接用距离限制当角度限制。3.3 完整代码与输入输出把上面的函数拼起来再加上一个暴力兜底的版本就是一份比较完整的参考代码。暴力版本主要用于处理角度超过 90 度的情况虽然慢一些但是原题实测时角度超过 90 度的用例通常不多不会成为性能瓶颈。#include bits/stdc.h using namespace std; const double PI acos(-1.0); const double EPS 1e-9; long long gcdll(long long a, long long b) { return b ? gcdll(b, a % b) : a; } vectorlong long getPrimeFactors(int n) { vectorlong long fac; for (long long p 2; p * p n; p) { if (n % p 0) { fac.push_back(p); while (n % p 0) n / p; } } if (n 1) fac.push_back(n); return fac; } long long countCoprime(int x, int lo, int hi) { if (lo hi) return 0; vectorlong long fac getPrimeFactors(x); int m (int)fac.size(); long long total hi - lo 1; for (int mask 1; mask (1 m); mask) { long long prod 1; int bits 0; for (int i 0; i m; i) { if (mask (1 i)) { prod * fac[i]; bits; } } long long cnt hi / prod - (lo - 1) / prod; if (bits 1) total - cnt; else total cnt; } return total; } long long solveFirstQuadrant(double angle, double dist) { if (angle 0) angle 0; if (angle 90) angle 90; long long ans 0; int R (int)ceil(dist); double rad angle * PI / 180.0; double tanA tan(rad); for (int x 1; x R; x) { double rem dist * dist - 1.0 * x * x; if (rem -EPS) continue; int yMaxDist (int)floor(sqrt(max(0.0, rem)) EPS); int yMaxAngle; if (angle 89.9999999) { yMaxAngle yMaxDist; } else { yMaxAngle (int)floor(1.0 * x * tanA EPS); } int yMax min(yMaxDist, yMaxAngle); if (yMax 0) continue; if (x 1) { ans; if (yMax 1) ans countCoprime(x, 1, yMax); } else { if (yMax 1) ans countCoprime(x, 1, yMax); } } if (angle 90.0 - EPS dist 1.0 - EPS) { ans; } return ans; } long long solveBrute(double angle, double dist) { double a fmod(angle, 360.0); if (a 0) a 360.0; int R (int)ceil(dist); long long ans 0; for (int x -R; x R; x) { for (int y -R; y R; y) { if (x 0 y 0) continue; if (gcdll(abs(x), abs(y)) ! 1) continue; if (1.0 * x * x 1.0 * y * y dist * dist EPS) continue; double ang atan2((double)y, (double)x) * 180.0 / PI; if (ang 0) ang 360.0; if (ang a EPS) ans; } } return ans; } int main() { double a, d; while (cin a d) { if (a 0 d 0) break; long long ans; if (a 90.0 EPS) ans solveFirstQuadrant(a, d); else ans solveBrute(a, d); cout ans \n; } return 0; }这份代码里同时保留了快速版和暴力版。日常做题时我建议先用暴力版验证快速版的结果特别是在角度、距离都很小的随机数据上两个版本输出一致才说明快速版写对了。这种“双实现互相验证”的方法能帮你节省大量调试时间。3.4 用 45° / 10 距离验证结果我拿一组自测数据做验证输入角度 45 度、距离 10。按照上面代码的逻辑输出应该是 25。这 25 棵树的坐标分别是x 1 时(1,0)、(1,1)x 2 时(2,1)x 3 时(3,1)、(3,2)x 4 时(4,1)、(4,3)x 5 时(5,1)、(5,2)、(5,3)、(5,4)x 6 时(6,1)、(6,5)x 7 时(7,1)、(7,2)、(7,3)、(7,4)、(7,5)、(7,6)x 8 时(8,1)、(8,3)、(8,5)x 9 时(9,1)、(9,2)、(9,4)x 10 时没有可见坐标把这些坐标代入距离公式最远的是 (9,4)距离是 √97仍然小于 10。代入角度公式所有坐标都落在 0° 到 45° 之间。比如 (8,5) 的角度是 32° 左右满足要求。所以这个结果是可以手算验证的。再测一组距离 10、角度 90 度输出应该是 49。角度 90 度意味着整个第一象限都可以计入只受距离限制。这个数也可以用对称性大致校验第一象限半径为 10 的圆内约有 78 个格点可见密度约为 6/π²约等于 47.7再加上坐标轴上的额外修正49 是合理的。4. 常见问题与避坑实录4.1 浮点精度问题UVa 149 的输入是浮点数最烦人的就是边界。角度刚好是 45 度、距离刚好是某个整数时理论边界和浮点边界不一定一致。我建议在两种地方都加 EPS。第一种是在向下取整前加 EPS这样能把理论边界上的点保留住。第二种是在距离比较时加 EPS避免出现“明明距离正好等于 D却因为 0.0000001 的误差被判成超出范围”的问题。但 EPS 不是越大越好。如果你取 1e-6 甚至 1e-4那些真正在边界外很近的点也可能被误判进来。我的经验是 1e-9 对这个题足够既不会漏掉边界上的整数点也不会引入额外的错误点。4.2 gcd 函数的边界处理手写 gcd 时很多人会忽略数值 0 的情况。但 (1,0) 和 (0,1) 是整个答案里非常重要的两个点。标准欧几里得算法递归到b 0时返回a所以 gcd(1,0) 1gcd(2,0) 2。我建议显式测试这三个边界gcd(1,0)、gcd(0,1)、gcd(0,0)。如果直接用__gcd(abs(x), abs(y))这种内置函数通常也能得到正确结果但在上传代码之前最好确认你的编译环境支持。为了避免潜在问题我更喜欢自己写一个 3 行的递归函数。4.3 容斥和优化的细节容斥原理本身不难难在边界条件的转换。当我统计 y 从 1 到 yMax 之间与 x 互质数量时lo从 1 开始因为 y 0 这一项需要特殊处理。如果你直接把lo设为 0容斥公式里的除法会出问题因为 0 乘以任何数都是 0会干扰统计结果。另外x 的质因数分解可以提前做一次缓存。如果题目有多个测试用例同一个 x 会被反复计算把每个 x 的质因子列表存起来能节省不少时间。我的做法是初始化一个数组下标从 1 到 R一次性预处理所有质因子然后在主计数函数里直接查表。4.4 输入结束标记UVa 老题的输入格式经常用0 0作为结束标记。但也有一部分重构版数据用 EOF 结束。我给出的代码使用0 0作为终止条件这是最常见的写法。如果你在自己平台上测试发现没有终止标记把 main 里的 break 条件注释掉改成单纯读到 EOF 就行。千万别因为输入格式问题陷入死循环。调试的时候可以先构造一个非常小的输入比如角度 0、距离 1预期输出是 1也就是只能看到 (1,0)。这个最小用例能快速帮你验证程序的主流程是否正确。5. 从这题延伸出的通用思路5.1 可见性问题在格点中的一般形态UVa 149 的题目背景是森林但它的数学内核在竞赛题里极其常见统计一个区域内不可约分向量的数量。不管背景是月亮上的环形山、城市里的摩天大楼还是棋盘上的棋子只要观测点在格点原点遮挡条件就大概率等价于 gcd 条件。我已经遇到过好几道类似题有些把坐标系旋转 45 度有些把二维换成三维。换成三维之后可见条件会变成 gcd(x, y, z) 1实现思路仍然是一样的只是枚举和容斥的维度更高。所以做这题时值得把“坐标互质方向唯一”这个直觉消化成自己的武器。以后再看到任何“从原点看格点”类问题第一反应就应该是 gcd而不是傻傻地去算线段经过哪些点。5.2 进一步优化与关联题目如果输入距离很大比如到 10^5 量级预处理质因子和记忆化容斥结果会非常有用。更进一步可以借助莫比乌斯函数把“统计互质点对”统一成求和式用数论分块加速。对于第一象限、固定距离 R 的可见格点总数有一个非常漂亮的近似公式总的可见格点数量大约是 6R² / π²。这个数字可以用来验证代码在小数据时的正确性。比如第一象限半径 10 的可见点总数算出来是 49而 π·100/4 ≈ 78.5乘以 0.608 大约是 47.7加上坐标轴修正后和 49 非常接近。这些扩展知识不需要在 UVa 149 里全部用到但下次遇到数据范围更大的同类题时你就知道该往哪个方向思考了。最后再分享一个小技巧。我每次写这种网格可见性题都会先打印出所有可见坐标肉眼扫一遍看有没有缺失或者多余的点。比如 45 度、距离 10 的答案我见过不少错误实现输出 24、26 甚至 30多数都是因为漏了 (1,0)或者多算了 (10,0)。用坐标表去对照程序结果比单纯看答案数字更容易定位问题。刷题这件事很多坑不是靠聪明绕过去的而是靠这套笨办法一遍遍踩平之后才变成真正属于自己的经验。