重访天数计算

# 数日子，重访 —— Tony Finch 来源：https://dotat.at/@/2026-08-09-rata-die.html 许多年前我写过如何将公历日期转换为儒略日数（https://dotat.at/@/2008-09-10-counting-the-days.html）或类似计数，例如*Calendrical Calculations*中使用的*rata die*（https://en.wikipedia.org/wiki/Julian_day）。这个算法是 C 语言`mktime()`（https://pubs.opengroup.org/onlinepubs/9799919799/functions/mktime.html）函数的核心，该函数将分解后的日期时间转换为线性的`time_t`。我最近从 Ben Joffe（https://www.benjoffe.com/fast-date-64）那里得知我错过了一些技巧，我的旧代码本可以写得更好。这是一个改进版本（使用常规编号而非 C 风格编号）： `` if m > 2 { m -= 2; } else { m += 10; y -= 1; } y*365 + y/4 - y/100 + y/400 + m*979/32 + d - 336 `` - 主要思路（https://dotat.at/@/2026-08-09-rata-die.html#the-main-idea） - 儒略年（https://dotat.at/@/2026-08-09-rata-die.html#julian-years） - 格里高利修正（https://dotat.at/@/2026-08-09-rata-die.html#gregorian-correction） - 月份模式（https://dotat.at/@/2026-08-09-rata-die.html#the-month-pattern） - 纪元（https://dotat.at/@/2026-08-09-rata-die.html#the-epoch） - 定义域与值域（https://dotat.at/@/2026-08-09-rata-die.html#domains-and-ranges） - 闰年测试（https://dotat.at/@/2026-08-09-rata-die.html#leap-year-test） - 月份长度（https://dotat.at/@/2026-08-09-rata-die.html#length-of-month） ## 主要思路（https://dotat.at/@/2026-08-09-rata-die.html#the-main-idea） 公历中有一个有用的巧合。虽然每个月的长度看起来并不规则，但存在一个重复的 5 个月模式，当从三月开始看时会更明显，如下表所示。这个模式在二月结束时重置，而此时恰好是其第三次重复的中途——*巧合的是*，这也正是闰日出现的位置。因此上面代码的第一行调整了月份和年份编号，使得一月和二月被计入前一年的末尾，而这种巧合对齐正好发生在调整后的年份编号的分界处。我会在讨论第二行中相关部分时详细解释这种调整。 三月 31 天 四月 30 天 五月 31 天 六月 30 天 七月 31 天 八月 31 天 九月 30 天 十月 31 天 十一月 30 天 十二月 31 天 一月 31 天 二月 28 或 29 天 ## 儒略年（https://dotat.at/@/2026-08-09-rata-die.html#julian-years) 主公式的第一部分计算年份`y`开始之前的天数，包括普通年和闰日。 - `y * 365 + y / 4` 调整部分在一月和二月时将年份减一。其效果是，第 4 年的闰日被计为调整后的第 4 年（即三月）开始之前的一天，也就是恰好在正确的位置。我以前将这部分表达式合并为一项： - `y * 1461 / 4` Ben Joffe 指出，这样写时，该函数只能利用其输出数据类型 25% 的范围，因为对于非常大的年份数字，乘法会溢出。而且在现代 CPU 上，消除加法实际上并不会更快。 ## 格里高利修正（https://dotat.at/@/2026-08-09-rata-die.html#gregorian-correction） 下一部分修正当前年份之前的闰年数量。 - `- y/100 + y/400` 它的工作原理与儒略闰年计算基本相同，只是`y/4`可以简单地编译为移位操作，而这里需要更多技巧。正如 *Hacker's Delight*（https://en.wikipedia.org/wiki/Hacker%27s_Delight）所解释的，现代编译器会将`y/100`转换为乘法和移位： - `y * (1 << (N+2)) >> (N+2)` 之类的形式 也就是说，编译器使用除数倒数的定点表示。然后它使用公共子表达式消除来抑制第二次乘以 1/25。因此这两个除法被转换为一次宽乘法加两次移位。很巧妙。 ## 月份模式（https://dotat.at/@/2026-08-09-rata-die.html#the-month-pattern) 下一部分计算在当前（调整后的）年份中，月份`m`开始之前的天数。 - `m * 979 / 32` 我以前是用 5 个月重复模式的天数来写的： - `m * 153 / 5` 但编译器不太容易将其优化得很好（clang 会使用两次乘法而不是一次），而且编译器不知道`m`的范围有限，所以我们可以手工将其转换为乘法和移位。979/32 == 30.59375，足够接近精确值 153/5 == 30.6。这两个表达式都能产生正确的 5 个月长短交替模式，但该模式不一定与正常的月份编号对齐。（上面两个表达式需要不同的调整。）我们将三月移到编号 1，正好在模式开始之前。当计算四月之前的天数时，我们得到比三月多 31 天；当计算五月之前的天数时，我们得到比四月多 30 天，以此类推。一月被调整为跟随十二月，以匹配调整后的年份编号。 `` m *979/32 diff ------------------- 1 30 March 2 61 31 April 3 91 30 May 4 122 31 5 152 30 6 183 31 7 214 31 8 244 30 9 275 31 10 305 30 December 11 336 31 January 12 367 31 13 397 30 `` ## 纪元（https://dotat.at/@/2026-08-09-rata-die.html#the-epoch) 由于日历从 1 开始计数，公历日期 0001-01-01 被编号为*rata die* 1。调整将一月变为月份 11，因此（如上表第二列所示）我们计算出调整后的第 0 年在一月之前有 336 天。我们需要减去这些额外的天数来抵消调整的影响。我们可以改变偏移量来选择不同的纪元，例如 MJD 纪元 1858-11-17 是*r.d.* 678576，Unix 纪元 1970-01-01 是*r.d.* 719163。 ## 定义域与值域（https://dotat.at/@/2026-08-09-rata-die.html#domains-and-ranges) 在我以前的 C 代码中，我随意使用了`int`，这误导性地暗示它适用于公元 1 年之前的格里高利历日期。然而，在常见 CPU 和低级编程语言上，有符号除法和取模（https://en.wikipedia.org/wiki/Modulo）是向零截断的，但这个算法需要欧几里得除法或向下取整除法（对于正除数，两者是等价的）。因此最好使用`u32`进行这些计算。（当代码使用无符号整数时，编译器也能做得更好。）要支持负年份，可以加上 400 年的倍数，将第 0 年移到`u32`范围的中间，然后从返回值中减去该倍数，以产生有符号的天数计数。 ## 闰年测试（https://dotat.at/@/2026-08-09-rata-die.html#leap-year-test) Ben Joffe 还研究了快速闰年测试（https://www.benjoffe.com/fast-leap-year）。我最喜欢的新版本是： `` if y % 25 == 0 { y % 16 == 0 } else { y % 4 == 0 } `` 我认为它即使不是最快的，也是最简洁的。注意 25 * 16 == 400，所以如果它同时能被 25 和 16 整除，那么就是闰年；否则，如果它能被 4 整除但不能被 25 整除，也是闰年。这种检查的经典版本首先测试能否被 100 整除。依赖分支预测的 CPU 在 99% 的情况下都能正确预测：比先测试能否被 4 整除的 75% 要好得多！但现在这个`if`会被编译成 CMOV 或 CSEL（因此分支预测无所谓），而且除以 25 比除以 100 更容易编译优化。 ## 月份长度（https://dotat.at/@/2026-08-09-rata-die.html#length-of-month) 促使我重新审视这段代码的原因是：当值具有有限范围（如上文`m*979/32`）和/或我们不依赖乘法的精确结果时，有可能得出更简单的乘法和移位优化。我以前写过这样的代码来计算月份长度： `` if m == 2 { 28 + is_leap_year(y) as u32 } else { 30 + (m * 275 % 9 > 3) as u32 } `` 乘法和移位的思路让我找到了这个针对二月以外月份的替代方案： `` 30 + (m * 7 % 16 < 9) as u32 `` 我是通过编写一个暴力程序找到它的，该程序尝试连续的位掩码宽度和乘数，直到找到一个情况，使得所有长月产生的结果都大于所有短月，或者反过来（就像这个获胜者一样）。但还有一个更简洁的表达式，据称出自 Matthias Kretz 博士： `` 30 | (m ^ (m >> 3)) `` 这用到了两个技巧： - 月份编号的 1 位与八月（月份 8）之前奇数/偶数长短模式相匹配，此时相位翻转。翻转是通过使用 8 位来切换 1 位实现的。 - 与 30 按位或会设置第 2、4、8、16 位，因此月份编号的更高位无关紧要。 它只编译成两条 ARM 指令： `` eor w0, w0, w0, lsr #3 orr w0, w0, #0x1e `` 它太棒了，我真的很喜欢它比我的尝试好多少！