仅用三条指令判断闰年：一种 is_leap_year_fast 的技术解析

Falk Hüffner · 2025年5月15日

通过以下代码，我们只需约 3 个 CPU 指令即可判断年份 0 ≤ y ≤ 102499 是否为闰年：

bool is_leap_year_fast(uint32_t y) {
  return ((y * 1073750999) & 3221352463) <= 126976;
}

这是如何实现的呢？答案出乎意料地复杂。本文将解释它，主要目的是为了享受位运算的乐趣；最后，我会简要讨论它的实际用途。

以下是闰年检查的典型实现：

bool is_leap_year(uint32_t y) {
  if ((y % 4) != 0) return false;
  if ((y % 100) != 0) return true;
  if ((y % 400) == 0) return true;
  return false;
}

我们使用的是 proleptic Gregorian calendar，它将 Gregorian calendar 从 1582 年的引入向后扩展，并包括 0 年。 因此，我们不需要以任何不同的方式处理 1582 年之前的年份。 为了简单起见，我们忽略负年份并使用无符号年份。

优化标准方法

让我们首先进行一些简单的加速，以便获得一个良好的基准。 我不确定该归功于谁——这些技巧可能已经被多次独立地重新发明了。

我们可以用 (y % 25) != 0 替换 (y % 100) != 0：我们已经知道 y 是 2 的倍数，所以如果它也是 5 的倍数，那么它就是 2 ⋅ 5 = 100 的倍数。 类似地，我们可以用 (y % 16) == 0 替换 (y % 400) == 0：我们已经知道 y 是 5 的倍数，所以如果它也是 2 的倍数，那么它就是 5 ⋅ 2 = 400 的倍数。

bool is_leap_year1(uint32_t y) {
  if ((y & 3) != 0) return false;
  if ((y % 25) != 0) return true;
  if ((y & 15) == 0) return true;
  return false;
}

这很有用，因为我们现在可以用位掩码替换模 4 和模 16。 编译器实现者还熟知一个技巧，即如何降低模 25。 使用 gcc 编译 (x % 25) != 0 并翻译回 C，我们得到 x * 3264175145 > 171798691。 在典型延迟为 3 个周期的乘法和至少 20 个周期的模运算下，这是一个很大的改进。 我将只给出它是如何工作的直觉； 更多详细信息可以在以下位置找到：

这些魔术数字从何而来？ 我们有 232 ⋅ 19/25 = 3264175144.96（精确）。 因此，通过乘以 3264175145，我们近似计算乘以 (19/25) 的小数部分。 如果这是一个精确的乘法，对于 25 的倍数，我们将获得 0 的小数部分。 然而，我们乘以的数字高了 0.04，因此误差可能高达 0.04 ⋅ (232 - 1) = 171798691.8，这就是第二个魔术数字的来源。

对于 x % 100，此技巧的效果不佳，我们需要一个额外的修复指令，因此将 y % 100 减少到 y % 25 非常有帮助。

我们的闰年检查现在如下所示：

bool is_leap_year2(uint32_t y) {
  if ((y & 3) != 0) return false;
  if (y * 3264175145u > 171798691u) return true;
  if ((y & 15) == 0) return true;
  return false;
}

请注意，像 gcc 或 clang 这样的现代编译器将从 is_leap_year1 生成 something 像 is_leap_year2 这样的代码，因此在 C 源代码中这样做没有多大意义，但在其他编程语言中可能很有用。

此代码通常编译为分支程序集。 在实践中，输入通常是相当可预测的，因此这不一定是一件坏事。 如果我们想避免分支预测错误惩罚，但同时也会减慢最佳情况的速度，我们可以稍微调整一下，并获得 branchless code（以及一个像样的代码高尔夫候选者）：

bool is_leap_year3(uint32_t y) {
  return !(y & ((y % 25) ? 3 : 15));
}

如果您想查看更多日历计算加速，请务必查看 Jacob Pratt 的 [https://jhpratt.dev/blog/optimizing-with-novel-calendrical-algorithms/](https://hueffner.de/falk/blog/)。

寻找一种位运算方法

我们是否可以通过放弃所有输入的正确性来改进闰年计算？ 毕竟，我们通常不关心 3584536493 年是否为闰年； 事实上，Python、C# 和 Go 仅支持 0（或 1）到 9999 年（此时相对于季节的漂移已经超过 4 天）。 我的想法是，如果存在更短的东西，它基本上看起来像一些使用魔术常数的奇怪哈希，所以我想要尝试一些小的形式并通过暴力破解来猜测常量。 (y * f) <= t 这种形式看起来很有用，但还不够强大。 我的候选方案之一是添加一个掩码：((y * f) & m) <= t。 现在我们有 96 位需要猜测，这仅靠暴力破解是不可行的。 让我们使用 z3，这是一个支持位向量约束的求解器，非常适合此：

import z3
BITS = 32
f, m, t, y = z3.BitVecs('f m t y', BITS)
def target(y):
  return z3.And((y & 3) == 0, z3.Or(z3.URem(y, 25) != 0, (y & 15) == 0))
def candidate(x):
  return z3.ULE((x * f) & m, t)
solver = z3.Solver()
solver.add(z3.ForAll(y, z3.Implies(z3.ULE(y, 400),
                  candidate(y) == target(y))))
if solver.check() == z3.sat:
  print(f'found solution: {solver.model()}')
else:
  print('no solution found')

在几秒钟内，这会找到在非平凡年份范围内产生正确结果的常量。 扩展范围，我最终找到了在 0 年到 102499 年范围内产生正确结果的常量，大约花费了半小时的计算时间，并证明了这是 32 位的最佳结果：

bool is_leap_year_fast(uint32_t y) {
  const uint32_t f = 1073750999u;
  const uint32_t m = 3221352463u;
  const uint32_t t = 126976u;
  return ((y * f) & m) <= t;
}

解释

它是如何工作的？ 令人惊讶，几乎是神奇的，我们可以将所有这些计算压缩到三个指令中。 然而，上面的考虑为我们提供了理解它的大部分工具。

以下是我们的二进制常量，其中四个相关的位范围用字母表示：

  A         D
f: 01000000000000000010001111010111
m: 11000000000000011111000000001111
t: 00000000000000011111000000000000
         B      C

让我们首先考虑对于乘积 p := y ⋅ f，用 m 进行掩码并与 t 进行比较会产生什么影响。 在块 A 中，t 中的位为 0，因此，只要 p 中设置了 A 中的任何位，结果就为假。 否则，块 B 变得相关。 这里 t 中的所有位都为 1，因此只要 p 中取消设置了 B 中的任何位，结果就为真。 否则，对于块 C，我们要求 p 中取消设置所有位。 这样，一系列位范围的比较都合并为一个 <=。

因此，我们可以将 is_leap_year_fast 重写为如下形式：

bool is_leap_year_fast2(uint32_t y) {
  uint32_t p = y * 1073750999u;
  const uint32_t A = 0b11000000000000000000000000000000;
  const uint32_t B = 0b00000000000000011111000000000000;
  const uint32_t C = 0b00000000000000000000000000001111;
  if ((p & A) != 0) return false;
  if ((p & B) != B) return true;
  if ((p & C) == 0) return true;
  return false;
}

这看起来很像 is_leap_year2！ 事实上，这三个条件具有完全相同的目的。 我们将展示

  1. 当 (y % 4) != 0 时，(p & A) != 0 触发；
  2. 当 (y % 100) != 0 时，(p & B) != B 触发；
  3. 当 (y % 16) == 0 时（因此 (y % 400) == 0，因为我们已经知道 y 是 25 的倍数），(p & C) == 0 触发。

两个简单的情况：（1）和（3）

(1)：f 中 A 中的 1 位将 y 的两个低位复制到 p 中的 A 中。 这不能被与 D 中的位乘法的结果破坏：我们可以获得的最大值是 102499 ⋅ (f & D) = 940428325，它只有 30 位。 因此，检查 A 在 p 中是否为零等效于检查 y 是否为 0 模 4。

(3)：检查最低 4 位中是否没有在 p 中设置任何位，就是检查 p 是否为 0 模 16。 但是，我们实际上想检查 y。 这没问题：仅查看 f 的最低 4 位就足够了，并且 f 在那里是 11112 = 15。 乘以 15 = 3 ⋅ 5 不会引入额外的因子 2，因此 16 的可除性不会改变。

有趣的情况：（2）

接下来，让我们尝试找出哪些数字 p & B ≠ B。 为此，f & A 中的 1 位不起作用，因此请考虑 f & D 中的位。 它们是 100011110101112 = 9175。 让我们看看测试适用于哪些数字。

>>> B = 0b00000000000000011111000000000000
>>> s = [y for y in range(5000) if ((y * 9175) & B) == B]
>>> for i in range(0, len(s), 16): print(*(f'{n:4d}' for n in s[i:i+16]))
 14  57  71 100 114 157 171 200 214 257 271 300 314 357 371 400
 414 457 471 500 514 557 571 600 614 657 671 700 714 757 771 800
 814 857 871 900 914 957 971 1000 1014 1057 1071 1100 1114 1157 1171 1200
1214 1257 1271 1300 1314 1357 1371 1400 1414 1457 1471 1500 1514 1557 1571 1600
1614 1657 1671 1700 1714 1757 1771 1800 1814 1857 1871 1900 1914 1957 1971 2000
2014 2057 2071 2100 2114 2157 2171 2200 2214 2257 2271 2300 2314 2357 2371 2400
2414 2457 2471 2500 2514 2557 2571 2600 2614 2657 2671 2700 2714 2757 2771 2800
2814 2857 2871 2900 2914 2957 2971 3000 3014 3057 3071 3100 3114 3157 3171 3200
3214 3257 3271 3300 3314 3357 3371 3400 3414 3457 3471 3500 3514 3557 3571 3600
3614 3657 3671 3700 3714 3757 3771 3800 3814 3857 3871 3900 3914 3957 3971 4000
4014 4057 4071 4100 4114 4157 4200 4214 4257 4300 4314 4357 4400 4414 4457 4500
4514 4557 4600 4614 4657 4700 4714 4757 4800 4814 4857 4900 4914 4957

正如所期望的那样，100 的倍数都在这里，但还有很多其他数字。 只要它们都不是 4 的倍数，这就不是问题，因为这些数字已经在上一步中被过滤掉了。 此外，0 丢失了，但这也不是问题，因为 0 也是 400 的倍数。

让我们尝试理解这个模式。 乍一看，它看起来非常简单：我们有 *14、*57、*71 和 *00。 然而，从 4171 开始，*71 消失了（你注意到了吗？）。 我们稍后也会出现新的模式。 让我们分析一下。 以下 Python 程序

def test(y):
  B = 126976
  return ((y * 9175) & B) == B
active = set()
for y in range(120000):
  r = y % 100
  if test(y):
    if r not in active:
      print(f'{y:6}: started *{r:02}')
      active.add(r)
  else:
    if r in active:
      print(f'{y:6}: stopped *{r:02}')
      active.remove(r)

输出 ``` 14: started *14 57: started *57 71: started *71 100: started *00 4171: stopped *71 32843: started *43 36914: stopped *14 65586: started *86 69657: stopped *57 98329: started *29 102500: stopped *00


因此，从 102500 开始，我们不再捕获 100 的倍数，这解释了为什么 102499 是 `is_leap_year_fast` 返回正确结果的最后一个数字。 我们还看到，在此之下，除了 100 的倍数外，没有一个数字是 4 的倍数（方便地，我们只能通过最后两位十进制数字来检查）。 如果我们相信这种蛮力枚举，这就完成了条件 (2) 的证明； 但让我们继续更好地理解为什么我们得到的确切数字。

让我们首先看看为什么我们首先得到 100 的倍数。 因子 9175 接近 17 位定点表示形式中 1/100 的倍数：
2<sup>17</sup> ⋅ 7/100 = 9175.04（精确）。
将 100 的倍数乘以 9175.04，将在位 17 及以上的位置得到一个整数（7 的倍数），并在位 17 以下的位置得到 17 个零位。 例如
9175.04 ⋅ 500 = 100011000000000000000002，其中 1000112 = 35 = 5 ⋅ 7。
将 100 的倍数乘以 9175 得到稍微少一点：
9175 ⋅ 500 = 100011000000000000000002 − 500 ⋅ 0.04 = 100010111111111111011002。
通常，从以很多零结尾的数字中减去一点点会创建一个以很多一结尾的数字，除了最后。 在这里，我们检查 _B_ 中的 5 位。 对于 _y_ 是 100 的倍数，这些保证都是 1，直到累积误差达到 _B_ 的下端，这仅在 _y_ = 2<sup>17</sup> / 0.04 = 102400 之后才会发生，因此这很合适。

现在 14、57 和 71 等其他数字来自哪里？ 让我们从不同的角度来看。 我们有 9175 = 2<sup>17</sup> ⋅ 0.06999969482421875（精确），并且 _B_ = 2<sup>17</sup> ⋅ 0.96875，因此
_p_ & _B_ = _B_ ⇔ {_y_ ⋅ 0.06999969482421875}  ≥ 0.96875，其中 {_x_} 是 _x_ 的小数部分 ⇔ 6.999969482421875 _y_ mod 100  ≥ 96.875

这是另一种看待为什么接受 100 的倍数的方式：对于它们，7 _y_ mod 100 为 0，因此 6.999969482421875 _y_ mod 100 略低于 100，并且仅在 _y_ = (100 − 96.875) / (7 − 6.999969482421875) = 102400 之后降至 96.875 以下。

为了理解出现在我们的序列中的其他数字，让我们首先考虑如果我们在此不等式中得到一个干净的 7，解决方案会是什么：
7y mod 100  ≥ 96.875 ⇔ 7y mod 100  ∈ {97, 98, 99}。
为了找到这个的解，我们首先需要 7 模 100 的模逆，即数字 _x_ 使得 7 _x_ mod 100 = 1。 我们可以使用扩展欧几里得算法来计算它，或者只是使用 [some online calculator](https://hueffner.de/falk/blog/<https:/www.wolframalpha.com/input?i=modular+inverse+of+7+modulo+100>)，它告诉我们结果是 43。 然后解是 43 ⋅ 97 (mod 100)、43 ⋅ 98 (mod 100) 和 43 ⋅ 99 (mod 100)，分别为 71、14 和 57 (mod 100)。 这解释了为什么我们最初看到 *14、*57 和 *71 形式的数字。 这也解释了为什么我们在 4071 之后停止看到例如 *71：虽然 7 ⋅ 4171 = 29197，但我们有 6.999969482421875 ⋅ 4171 = 29196.872711181640625，这（模 100）小于 96.875。 类似地，32843 弹出是因为累积误差 (7 − 6.999969482421875) ⋅ 32843 = 1.002288818359375 超过 1。 如果在这里付出更多的努力，我们可以手动重现上面 `Python` 程序的输出，并检查这些数字中没有一个是 4 的倍数。

## 扩展到其他位宽

现在我们了解了这个技巧是如何工作的，我们可以尝试找到其他位宽的参数。 变量部分是块 _B_ 的位置和 _f_ & _D_ 中 100 的分数。

uint64_t score(uint64_t f, uint64_t m, uint64_t t) { for (uint64_t y = 0; ; y++) if ((((y * f) & m) <= t) != is_leap_year(y)) return y; }

int main() { uint64_t best_score = 0; for (int k = 0; k < BITS; k++) { for (int k2 = 0; k2 < k; k2++) { uint64_t t = (1ULL << k) - (1ULL << k2); uint64_t m = (0b11ULL << (BITS - 2)) | t | 0b1111; for (int n = 0; n < 100; n++) { uint64_t f = (0b01ULL << (BITS - 2)) | (((1ULL << k) * n) / 100); uint64_t new_score = score(f, m, t); if (new_score > best_score) { printf("%llu %llu %llu: %llu (%d %d %d)\n", f, m, t, new_score, k, k - k2, n); best_score = new_score; } } } } return 0; }


对于 `BITS = 64`，在大约 7 分钟内，我们找到 _f_ = 4611686019114582671、_m_ = 13835058121854156815、_t_ = 66571993088，这对于 _y_ = 5965232499 是正确的。 这很好，因为 5965232499 > 2<sup>32</sup>，因此可以使用此变体正确测试任何 32 位年份。

这是 64 位可以达到的最高水平吗？ 也许还有一些其他常数效果更好？ 我无法立即找到一种方法来证明它，因此我采用了经过验证的“让其他人为我做”的方法，将其 [posting it to the Code Golf StackExchange](https://hueffner.de/falk/blog/<https:/codegolf.stackexchange.com/q/275505/116815>)。 事实上，仅在 1 小时后，用户 `ovs` 发布了一个非常好的结果，两天后，用户 `Exalted Toast` 发布了一个 [proof](https://hueffner.de/falk/blog/<https:/codegolf.stackexchange.com/a/275541>)，即 5965232499 确实是 64 位的最佳范围，也使用了 `z3` 求解器。

## 基准测试

获得有意义的基准测试很棘手，因为该函数仅占用很短的时间，而且对于分支版本，执行时间取决于输入的模式。 让我们尝试两个极端情况：始终是 2025 年，以及完全随机的年份。 这些是在 `i7-8700K` (`Coffee Lake`, 4.7 GHz) 上 [benchmark](https://hueffner.de/falk/blog/<benchmark.cc>) 的结果，使用 `g++ -O3 -fno-tree-vectorize` 编译：

2025 (ns) | 随机 (ns)  
---|---  
`is_leap_year` | 0.65 | 2.61  
`is_leap_year2` | 0.65 | 2.75  
`is_leap_year3` | 0.67 | 0.88  
`is_leap_year_fast` | 0.69 | 0.69  

这里突出显示了一些奇怪的事情：
  * 在随机情况下，`is_leap_year2` 实际上比 `is_leap_year` 稍慢。 这令人惊讶，因为 `y % 100` 需要比 `is_leap_year2` 中的技巧多一个指令。
  * 对于随机数据，`is_leap_year3` 比对于固定值稍慢。 这令人惊讶，因为它没有做任何分支。

除了基准测试很难之外，我对此没有任何解释。

对于随机数据，新函数 `is_leap_year_fast` 比标准实现快 3.8 倍，对于完全可预测的输入，速度慢约 6%。 总的来说，这似乎非常可靠。

## 结论

那么，这真的值得吗？ 我们是否应该用这个技巧替换例如 [CPython datetime implementation](https://hueffner.de/falk/blog/<https:/github.com/python/cpython/blob/5cdd49b3f4cbdcf0472a65fd0c723912c3d48211/Modules/_datetimemodule.c#L416>)？ 嗯，这取决于。 实践中查询的大多数年份将是当前年份，或者至少是相当可预测的，在这种情况下，我们没有很大的优势。 为了完全证明更改的合理性，理想情况下，我们将有一个使用实际数据的基准，该基准使用闰年检查作为子例程，而不仅仅是微基准测试。 我很乐意听到任何此类结果！