ARTICLE · INTELLIGENCE

战地情报 · 详情页

来自尧图项目组的一线实战观察与深度解析

灯泡开关问题:从暴力模拟到常数级解的高性能计算实践

灯泡开关问题:从暴力模拟到常数级解的高性能计算实践 先说结论n个灯泡从关闭开始第i轮翻转所有编号为i的倍数的开关最终亮着的灯泡数量是floor(√n)。这不是抖机灵而是整道题从暴力模拟到数学优化自然收敛的结果。作为一份高性能计算笔记我想把这条路完整走一遍——包括为什么暴力解在大规模场景下必死、约数奇偶性如何完成降维、常数级代码里浮点精度那些坑以及这些问题在真实性能优化里对应的思维方式。LeetCode 319 的解法说一句就完了但背后那套“先归约再优化”的 HPC 方法论值得展开。这篇文章适合做性能优化的人、刷题刷到这道题但对“常数级”三个字感兴趣的开发者也适合所有习惯把“能不能更省”挂在嘴边的人。灯泡开关问题规模一大人人都会意识到硬算是走不通的——但怎么一步步推导出数学解怎么把数学解落成无错的工程代码这才是真正有价值的干货。1. 先别急着开超算暴力解法在灯泡问题里到底要消耗多少资源1.1 问题是如何从“开关棋盘”变成“抽象序列”的先把问题形式化有n个灯泡编号1..n初始全部关闭。第i轮操作是翻转所有能被i整除的编号对应的灯泡状态。也就是说第i轮翻转的灯泡集合是{i, 2i, 3i, ...}。一直做到第n轮问最终有多少灯泡亮着。最直接的模拟思路是这样的uint64_t brute_force(uint64_t n) { std::vectorchar bulb(n 1, 0); for (uint64_t i 1; i n; i) for (uint64_t j i; j n; j i) bulb[j] ^ 1; uint64_t ans 0; for (uint64_t j 1; j n; j) ans bulb[j]; return ans; }这个代码的麻烦不是“慢一点”这么简单。外层跑n轮第i轮的内层循环次数是n / i所以总翻转次数是调和级数T(n) n/1 n/2 n/3 ... n/n n * H_n ≈ n (ln n γ)这里γ≈0.5772是欧拉常数。当n10^9时T ≈ 2.1×10^10也就是大约两百一十亿次翻转。两百一十亿次字节异或单看指令数量已经不小但真正的杀手是后面那笔访存账。1.2 用带宽和延迟估算暴力解法的“真实成本”我们在 HPC 里评估一个循环能不能算第一件事不是数指令而是算“内存流量”。灯泡模拟里每个字节状态至少要经历一次读、一次写读旧状态、异或、写回新状态。于是理想情况下每次翻转产生2字节访存流量。把几个量级摊开看n状态数组体积总翻转次数理想内存流量单机可行性1e91 GB约 210 亿约 420 GB一台大内存机器勉强能放但跑起来很慢1e121 TB约 2.8×10^13约 56 TB需要分布式存储跨节点通信成本爆炸1e151 PB约 3.5×10^16约 70 PB数据中心级贵到离谱新一代 DDR5 内存带宽大约在 50~100 GB/s 这个量级。照这个估算n1e9时最小访存时间也要四五秒真实运行因为 cache miss、TLB 失效、stride 访问导致缓存完全发挥不出作用实际时间大概率再翻几倍。请注意这只是单轮扫描的乐观下限还没算分支、循环开销和调度抖动。更有意思的是访问模式。内层循环j i当i很大的时候相邻两次访问的地址相隔很远比如i接近n/2时每次翻转都要从内存加载一个全新的 cache line之前缓存里存的东西全浪费了。这种“跨越式访问”在 HPC 里叫 stride 访问是最伤性能的模式之一。要解决它完全可以用#pragma omp parallel for把外层循环并行化但线程能摊掉一部分工作量摊不掉内存流量本身。这是典型的“内存墙”问题——算力再大数据搬不动就是搬不动。换句话说这道题如果走暴力模拟你需要的根本不是更强的计算能力而是更聪明的问题转换。2. 从翻转次数到约数奇偶性数学上的一次降维2.1 每个灯泡被翻转的次数是一个数论函数回到第 1 章的分析一个灯泡最终是亮还是灭只看它被翻转了多少次翻转奇数次从关变开翻转偶数次状态不变。那么灯泡j在整个过程中到底被翻转了几次它在第i轮被翻转当且仅当i能整除j。所以灯泡j的翻转次数就是j的所有正整数约数的个数标准记号是τ(j)。于是问题第一次从“过程模拟”变成“性质统计”统计1..n中哪些数满足τ(j)为奇数。这一步其实就是 HPC 里最常用的一招——不要逐步演化系统而是直接找系统中的“守恒量”或“解析特征”。比如粒子模拟里我们有时不需要每帧追踪所有相互作用而是用统计分布代替个体行为。2.2 约数配对为什么只有完全平方数能亮现在看τ(j)的奇偶性。拿12举例它的约数是1, 2, 3, 4, 6, 12。除了1和122和6是一对乘积等于123和4也是一对乘积还是12。一般情况下如果d是j的约数那么j/d也一定是约数因为d * (j/d) j天然成立。所以除了满足d j/d的特殊情况约数永远成对出现。成对意味着偶数个。唯一打破偶数配对的条件是d² j也就是说j是一个完全平方数。此时约数配对的“中间项”只有一个τ(j)会变成奇数。结论非常干净完全平方数编号的灯泡被翻转奇数次最终亮着非完全平方数编号的灯泡被翻转偶数次最终灭着。所以这道题从“模拟 n 轮开关”变成了“数一数 1..n 里有多少个完全平方数”。2.3 从模拟到计数思维转变的代价一旦换成计数视角答案就成了求最大的整数k满足k² ≤ n也就是floor(√n)。从O(n log n)的翻转量直接掉到O(1)的开方。这个降幅让我想起一个经典比喻你原本要手动清点一排抽屉里的硬币数学优化告诉你“这一排抽屉只有两类一类会亮一类不会”那你还点什么呢直接数类就行了。这种“过程型 → 计数型”的归约在真实高性能计算里天天出现。比如流体模拟中的某些算子原始定义是逐个网格点做 NYU 求和但通过格林函数解析展开可以直接得到远场近似再比如稀疏矩阵乘法如果知道矩阵结构满足某种低秩性质就压根不用建全矩阵。灯泡开关问题只是把同一个思路浓缩成一道简答题。3. O(1)算法的工程细节整数平方根没那么天真3.1 从数学答案到可运行代码为什么必须做边界修正数学结论是floor(√n)落到代码小学生都会写return (int)sqrt(n)。可这道题一旦进入“高性能计算”语境就不能这么草率。std::sqrt接收和返回的是double而double的有效数字只有 53 位。整数在2^53以内时double可以精确表示每一个整数超过2^53就有两个相邻整数被折叠成同一个浮点数。你给std::sqrt传一个快接近2^63的uint64_t它先要把这个数转成double这一步就可能发生舍入。算出来的平方根和真实数学平方根在整数边界上可能相差 1。更隐蔽的问题出现在“完全平方数恰好紧邻边界”时。比如某个n的数学平方根是k但浮点运算结果落在k - 1那边直接截断成uint64_t就错了。所以任何依赖浮点做精确整数运算的代码都必须用整数乘法回验。3.2 一个带边界修正的整数平方根实现我最常用的版本是这样的uint64_t isqrt_u64(uint64_t n) { uint64_t r static_castuint64_t(std::sqrt(static_castdouble(n))); while ((__uint128_t)(r 1) * (r 1) n) r; while ((__uint128_t)r * r n) --r; return r; }两个while循环各自最多跑一次或两次修正成本基本可以忽略不计。关键在于比较时用了__uint128_t防止(r1)*(r1)在 64 位无符号整数里溢出。r最大约2^32平方后约2^64刚好碰到uint64_t的上限不额外扩展类型就有边界风险。为什么不用sqrtl一劳永逸long double在 x86 上通常是 80 位扩展精度64 位尾数对 64 位整数平方根确实足够。但在某些平台上long double跟double没有区别跨平台代码不能赌这个。__uint128_t是 GCC/Clang 的扩展MSVC 不支持如果你在 Windows 上工作可以用两个uint64_t手动比较乘积或者干脆用整数二分兜底。3.3 常数级不是免费午餐乘法验证也要算成本有读者会说都 O(1) 了还纠结什么但在高性能计算里O(1) 和 O(1) 的差别可能就是工程成败。比如你要对一亿个n分别求灯泡开关答案每个都做一次std::sqrt加两次乘法修正累积起来就值得算。一次uint32_t * uint32_t的 64 位乘法在 x86 上大概是几条指令的事std::sqrt走硬件指令通常延迟十几个周期、吞吐数个周期一条。如果你在循环里调一百万次一次开方多花 10 个周期总量就是一千万个周期折算到 3 GHz 的机器上大约几毫秒。这个量级普通应用无所谓但放在批量微基准测试或者 GPU kernel 里就值得用更激进的策略。这也正是 HPC 的真实面貌算法复杂度降到位只是获得了参赛资格接下来常数级的小数点后优化才是一流性能团队拉开差距的地方。4. 把常数级再压一压sqrt指令、精度损失与边界检查4.1 当 n 超过双精度有效位数时会发生什么double的 53 位有效数字把它能“无损分辨”的整数上限钉死在2^53 ≈ 9.007×10^15。在那个数字之上double能表示的最小间隔变成 2、4、8……这意味着当你执行static_castdouble(n)时n可能已经被舍入到附近的偶数。std::sqrt在这个被打磨过的输入上做运算结果自然可能偏离真实平方根。一个实际例子假设n接近2^63真实平方根在3037000499附近。浮点开方的相对误差通常在几个 ULP 以内但ULP在这个量级对应的绝对误差是几十到几百。如果没有修正循环返回结果可能差出好几个数。好在我们有第 3.2 节的while修正往下减一个或者往上加一个把结果拉回数学真值。这里要澄清一个常见误解**修正循环不是“处理异常情况”它本身就是算法的必要组成部分。**只要你的输入可能超过2^53浮点结果就不能直接信任。哪怕概率只有百万分之一数值代码的失败概率也必须为 0。4.2 从 sqrt(double) 到整数二分、牛顿迭代到底怎么选“常数级解法”内部也有多条路线。我按工程场景整理成一个表方法精确保真范围大致成本特征适用场景std::sqrt(double) 乘法修正≤2^53 放心更大也能靠修正兜底极快硬件指令通用首选std::sqrtl(long double)x86 上 64 位整数基本安全比double慢几倍可读性优先的跨平台谨慎版整数二分完全精确无浮点依赖O(log n)几十次乘法循环嵌入式、无硬件 sqrt、需要可证明正确时牛顿迭代 查表初值可以达到 64 位整数全范围正确2~3 次迭代非常适合向量化GPU/DSP 上的批量求根低精度近似指令如 rsqrt通常只有几到十几位精度延迟低、吞吐高图形、近似搜索、对误差不敏感的 HPC 算法在 x86 上sqrtsd指令的延迟大约是十几周期吞吐则取决于微架构。相比一次 L1 命中约 4 周期和一次整数乘法约 3 周期它不算贵但也不算免费。如果你的热循环里真的每一轮都要开方可以考虑用 SIMD 同时对多个double求平方根。AVX2 的vsqrtpd、AVX512 的vrsqrt14ps这类指令把“常数级”的常数再除以向量宽度。4.3 批量场景下的向量化思路玩具问题里的工程感有人会问灯泡开关问题本身只需要一次开方谈向量化是不是小题大做我觉得不是。换个输入模式就合理了如果传给函数的是一组很长的n例如要统计考虑多个规模的灯泡开关数量你需要对一整个数组的n求isqrt这时候 SIMD 开方加修正就是正经的高性能优化。伪代码长这样void bulbs_batch(const uint64_t* in, uint64_t* out, size_t len) { #pragma omp simd for (size_t i 0; i len; i) { uint64_t r isqrt_u64(in[i]); out[i] r; } }编译器会尝试把循环里的sqrt向量化。但如果内部有分支修正向量化效果会打折。更激进的写法是先用近似指令计算一批初值再统一做两轮牛顿迭代最后用整数乘法把所有通道的结果拉回正确范围。细节极其吃平台这里先放一个方向。我强调一遍这个案例里正统路线永远是“数学降维优先指令级优化随后”。如果你试图用 AVX512 优化一个本来就要跑十几小时的暴力模拟那是本末倒置反过来先证明算法必须做一次开方再去榨干这次开方的每个周期顺序就对了。5. 灯泡问题教给我的高性能计算思考方式5.1 性能上界由什么决定复杂度、访存模式和向量宽度在高性能计算项目里我判断一段代码有没有救通常看三件事算法复杂度决定“总工作量”数据布局决定“访存模式”向量化宽度决定“理论峰值”。灯泡问题把这三件事全演示了一遍暴力翻转的复杂度是n log n工作量巨大这是第一层原罪。每个字节状态被反复跨越式访问访存模式极其低效这是第二层原罪。即便你把它并行化、向量化前面两层损失也补不回来。数学优化直接把第一层从n log n降到 1后面两层的优化空间才变得有意义。这也是 Amdahl 定律的朴素版本如果你的工作总量里有 99% 是无效翻转把剩下 1% 的开方提速一万倍也毫无帮助。5.2 先降算法复杂度再谈并行化我踩过的一个真实坑有一次优化一个数值积分模块最开始我把物理网格加密后直接上 OpenMP开了 32 个线程跑起来仍然慢得离谱。用 perf 一看内存带宽直接被写满线程数再翻倍只会更糟。后来我把被积函数做了符号化展开消掉一个嵌套求和也就是把每步计算量从 O(N) 降到 O(1)再上并行墙钟时间瞬间降了两个数量级。灯泡开关问题就是这件事的微缩模型。它反复提醒我同一个教训并行化是最后的优化手段不是第一反应。如果数学结构允许你消掉大量中间态那么任何调度器、缓存优化、线程库都接不住这个红利。5.3 简化与近似的边界这个案例里“近似”必须被精确回验第 4 章那个while修正本质上是“有界近似 精确回验”的组合拳。先用快的浮点指令得到一个近似根再用整数乘法判断它是否落在正确答案的邻域内最终把它拉回精确值。这种模式在 HPC 里常见到让人习惯性忽略低秩近似算法会用残差检查决定是否该增秩稀疏迭代求解器会用收敛判据决定何时停止蒙特卡洛会用方差估计决定样本量。所有近似都必须配一个“怎么知道自己错了”的机制。灯泡问题里的回验机制就是两个while循环几行字却是整套方案可靠性的地基。我自己后来在写大规模粒子平衡树代码时也把这种“近似计算 整数/统计回验”的套路带进去先用浮点快速排序算分箱位置再对边界元素做精确比较修正。性能和正确性一次性都要。6. 避坑清单从浮点修正到代码可读性如果只记住一条经验那就是“答案就是 sqrt(n)”是一个数学陈述不是代码规范。写成生产级代码时必须考虑溢出、舍入和平台差异。我把这几条容易踩的坑集中列出来直接(int)sqrt(n)返回。在n 2^53或接近完全平方数时可能差 1。最经典错误。修正时用uint64_t做乘法。(r1)*(r1)可能溢出 64 位。用__uint128_t或等价的手动高半区乘法。迷信sqrtl。x86 上没问题迁移到某些 ARM64 或 MSVC 环境就退化。循环热路径里不加分支修正。当批量输入时分支预测和 SIMD 退化会让理论上的 O(1) 变成实际上的大常数。把“算法复杂度低”当成“性能必然好”。复杂度是上限常量因子、缓存、向量化等工程因素依旧决定最终落地效果。还有一个容易被忽视的点代码可读性。灯泡开关问题的答案只有一行但你如果在团队代码库里写return std::sqrt(n);下一个人看见会一脸茫然。我通常会在上面留两行注释说明为什么答案等于完全平方数个数为什么需要修正。高性能计算的代码最难维护的不是寄存器分配而是未来有人要理解你的假设。现在回到最初那个问题为什么一个玩具题值得用“高性能计算笔记”来记录因为它的优化链条太典型了——暴力模拟、资源估算、数论归约、浮点近似、边界修正、向量化每一步都能在真实 HPC 项目里找到对应。灯泡开关只是恰好让这条链在一道题里全部跑了一遍。我的习惯是拿到一个计算密集任务先问三个问题这个问题有没有解析解没有的话能不能用统计/近似绕过最重的循环如果都必须算我能不能把数据布局改成顺序访问灯泡开关问题对第一个问题的回答是“有”所以后面两个压根不用启动。下次遇到看起来需要超多算力的问题不妨先停下来想想你要的到底是“算得快”还是“不用算”。
RELATED READING

延伸阅读

更多一线实战笔记与深度复盘,助您持续精进