蔡勒(Zeller)公式及其推导：快速将任意日期转换为星期数
0. 本文的初衷及蔡勒公式的用处
前一段时间，我在准备北邮计算机考研复试的时候，做了几道与日期计算相关的题目，在这个过程中我接触到了蔡勒公式。先简单的介绍一下蔡勒公式是干什么用的。
我们有时候会遇到这样的问题：看到一个日期想知道这一天是星期几，甚至看到一个历史日期或纪念日，我们想快速的知道这一天是星期几。对于这个问题，如果用编程的方式，应该怎么实现呢？你可能已经有思路了，比如你知道某个日期是星期几，把这个日期作为原点，然后计算目标日期和这个原点之间相差多少天，再除以 7 求余数，最后通过余数判断目标日期的星期数。通过这样的过程，你确实可以得到正确的结果，但这不够快速也不够优雅。对于这个问题，如果你懂得蔡勒公式，那就变得异常简单了，你只需要将年月日等数据代入公式，然后计算出结果，这一结果就是目标日期对应的星期数。
当我知道蔡勒公式之后，我觉得它非常有趣也很酷，所以我不仅希望掌握公式的使用，也希望可以掌握公式背后的推导过程。然而，当我在网上搜索相关的文章时，我发现几乎所有向我展示的博客（从零几年到最近的 19 年）大多是转载、复制于这篇文章（链接）：
星期制度是一种有古老传统的制度。据说因为《圣经·创世纪》中规定上帝用了六天时间创世纪，第七天休息，所以人们也就以七天为一个周期来安排自己的工作和生活，而星期日是休息日……
这篇文章质量很不错，讲解过程自然流畅，但是在一些细节上存在错误，有些推导步骤让人感到困惑。因此，当我掌握蔡勒公式后，很希望可以将我的理解输出出来，让想要学习蔡勒公式推导过程的人看到一些新的材料。好了，废话少说，我们开始吧。
1. 蔡勒公式的形式
如果你对公式的推导过程不感兴趣，只是希望使用蔡勒公式，那么只看此小节即可。蔡勒公式的形式如下：
其中：
- W 是星期数。
- c 是世纪数减一，也就是年份的前两位。
- y 是年份的后两位。
- m 是月份。m 的取值范围是 3 至 14，因为某年的 1、2 月要看作上一年的 13、14月，比如 2019 年的 1 月 1 日要看作 2018 年的 13 月 1 日来计算。
- d 是日数。
- [] 是取整运算。
- mod 是求余运算。
注意：这些符号在后面的推导中还会使用。举一个实际的计算例子：计算 1994 年 12 月 13 日是星期几。显然 c = 19，y = 94，m = 12，d = 13，带入公式：
由此可得 1994 年 12 月 13 日是星期二，通过查询日历可以验证正确性。
最后关于蔡勒公式，还需要做两点补充说明：
- 在计算机编程中，W 的计算结果有可能是负数。我们需要注意，数学中的求余运算和编程中的求余运算不是完全相同的，数学上余数不能是负数，而编程中余数可以是负数。因此，在计算机中 W 是负数时，我们需要进行修正。修正方法十分简单：让 W 加上一个足够大的 7 的整数倍，使其成为正数，得到的结果再对 7 取余即可。比如 -15，我可以让其加上 70，得到 55，再除以 7 余 6，通过余数可知这一天是星期六。
- 此公式只适用于格里高利历（也就是现在的公历）。关于历法的问题，大家有兴趣可以自行查阅。
下面是公式的推导。
2. 推导过程
推导蔡勒公式之前，我们先思考一下，如果我们不知道这一公式，我们如何将一个日期转化为星期数呢？
我们可能会很自然地想到：先找到一个知道是星期几的日子，把这个日期作为“原点”，然后计算目标日期和这个原点相差几天，将相差的天数对 7 取余，再根据余数判断星期数。举一个实际例子，比如我们知道 2019 年 5 月 1 日是星期三，把这一天当作原点，现在我们想知道 2019 年 5 月 17 日是星期几，显然这两个日期之间相差 16 天，用 16 除 7 余 2，因为原点是星期三，加上作为偏移量的余数 2，可知 2019 年 5 月 17 日是星期五。
从这个思路出发，经过优化扩展，我们就可以得到神奇的蔡勒公式了。首先，如果我们仔细观察一下可以发现，这个思路中比较麻烦的是计算相差天数（设为 \(D\)），我们可以把 \(D\) 的计算过程分解成三部分：
- \(D_1\)：从原点到原点所在年份末尾的天数。
- \(D_2\)：原点所在年份和目标日期所在年份之间所有年份的天数和。
- \(D_3\)：目标日期所在年份的第一天到目标日期的天数。
显然，\(D = D_1 + D_2 + D_3\)。如果我们把原点选择在某一年的 12 月 31 日，那么就可以省去 \(D_1\) 的计算了，因为原点恰好就是所在年份的最后一天。现在，\(D = D_2 + D_3\)。
我们再去观察上述思路中的实际例子，可以发现，因为原点是星期三，所以得到余数后，我们需要加上 3 才是正确的星期数。这启示我们可以把原点选定在星期日，这样算出来的余数是几就是星期几（余数 0 代表星期日）。
经过这样的分析。我们希望可以优化原点的日期，使其满足下面两个条件：
- 是某一年的 12 月 31 日。
- 是星期日。
我们按照现在使用的公历的规则逆推，可以发现公元元年的前一年的 12 月 31 日恰好是星期日，满足我们想要的两个条件，是一个完美的原点。
现在假设目标日期是 y 年 m 月 d 日，我们已经可以很容易的计算 \(D_2\) 了：
简单的解释一下。因为一年最少有 365 天，所以 \(D_2\) 至少是 \((y-1) \times 365\)。此外，因为闰年比平年多一天，我们还需要加上这些年份中闰年的数量。按照闰年的规则：每 4 年一闰，但每 100 年不闰，每 400 年又闰。可知闰年的数量为 \(\left[ \frac{y-1}{4} \right] - \left[ \frac{y-1}{100} \right] + \left[ \frac{y-1}{400} \right]\)。
现在，我们需要得到 \(D_3\) 的计算公式，这块要复杂一些。首先，不考虑闰年的话，每年中 2 月份天数最少，为 28 天。因此，我们不妨把每个月的天数看作 “28 + Excess”的模式，m 月之前所有月份的 Excess 之和为 Accum(m)，则 \(D_3 = 28 \times (m-1) + Accum(m) + d\)，并且我们可以得到这样一张表格：
| 月份 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 天数 | 31 | 28 | 31 | 30 | 31 | 30 | 31 | 31 | 30 | 31 | 30 | 31 |
| Excess | 3 | 0 | 3 | 2 | 3 | 2 | 3 | 3 | 2 | 3 | 2 | 3 |
| Accum | 0 | 3 | 3 | 6 | 8 | 11 | 13 | 16 | 19 | 21 | 24 | 26 |
仔细观察，可以发现 Excess 从 3 月份开始是 3、2、3、2、3 的循环，因此，当 \(m\geq 3\) 时，\(Accum(m)\) 的值的增幅也是 3、2、3、2、3 的循环。因为每 5 个月增加 13，所以把 \(\frac{13}{5}\) 作为系数；因为 \(Accum(m)\) 的值是离散的（都是整数），所以我们用取整运算符，得到：
我们将 \(x\) 的值取 1，2，3……，然后观察 \(f(x)\) 的值，可得下面这张表格：
| x | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 |
|---|---|---|---|---|---|---|---|---|---|---|---|
| f(x) | 10 | 13 | 15 | 18 | 20 | 23 | 26 | 28 | 31 | 33 | 36 |
我们可以发现，当 \(x \geq 4\) 时，\(f(x)\) 的值的增幅也是 3，2，3，2，3 的循环。也就是说 \(f(x)\) 的值的增幅（\(x \geq 4\)）与 \(Accum(m)\) 的值的增幅（\(m \geq 3\)）相同，这意味着 \(f(x)\) 与 \(Accum(m)\) 之间相差一个常数，我们随便带入一个具体的值计算：
可知相差的常数为 7。由此可得，当 \(m \geq 3\) 时，\(Accum(m)\) 的值的序列，等于当 \(x \geq 4\) 时，\(f(x) - 7\) 的值的序列。这样我们就得到了 \(Accum(m), m \geq 3\) 的函数形式：
这里多说两句，实际上，\(Accum(m)\) 的函数形式是不唯一的，使用其他的构造方法，可以得到形式不同的 \(Auccm(m)\)，只要符合要求即可。
进一步，我们可以得到 \(D_3\) 的函数形式：
其中，平年时，\(i = 0\)；闰年时，\(i = 1\)。这还不是 \(D_3\) 最完美的形式。我们继续分析，从 3 月份到 12 月份的 Excess 正好是两个 3、2、3、2、3 的循环，那么假如有第 13 月，想要继续保持这种循环规律，13 月的 Excess 值应该是 3。而 1 月份的 Excess 的值恰好是 3，所以我们不妨变通一下，把每年的 1 月、2月当作上一年的 13月、14 月。这样不仅仍然符合公式，而且 2 月份变成了上一年的最后一个月，也就是公式中 \(d\) 的部分，于是平闰年的影响也去掉了，\(D_3\) 的公式简化成了：
现在，我们已经得到了 \(D_2\) 和 \(D_3\) 的计算函数，由 \(D = D_2 + D_3\)，可知：
注意！这个公式离正确形式还差一小步。因为在当前的公式中，每年的 1 月和 2 月被当作上一年的 13 月和 14 月计算，因此当前公式中计算闰日的部分（\(\left[ \frac{y-1}{4} \right] - \left[ \frac{y-1}{100} \right] + \left[ \frac{y-1}{400} \right]\)）存在错误。举一个具体的例子，比如计算公元 4 年（闰年）3 月 1 日的星期数。在当前公式中，公元 4 年的 2 月被算作了公元 3 年的 14 月（换句话说公元 3 年变成了闰年），而公式中计算闰日的部分没有考虑这点，依然将公元 3 年当成平年计算，因此少算了一天。因此，计算闰日的部分应当改进，公式如下：
计算出 D 的值后，对 7 取模即可得到星期数。
根据同余定理，D 除以 7 所得的余数等于 D 式的各项分别除以 7 所得余数之和（余数之和大于等于 7 时，再对 7 取余），因此可以消去 D 式中能被 7 整除的项，进行化简：
简单说明一下：
显然，结果中的第一项是 7 的倍数，除以 7 余数为 0，因此 \((y-1) \times 365\) 除以 7 的余数其实就等于 \((y-1)\) 除以 7 的余数，我们只保留 \((y-1)\)就够了。化简过程中，其他被销去的项同理。
公式（2）还不是最简练的形式，我们还可以对年份进行处理。我们现在用公式（2）计算出每个世纪第一年 3 月 1 日的星期数，得到如下结果：
| 年份： | 1, 401, 801, … , 2001 | 101, 501, 901, … , 2101 | 201, 601, 1001, … , 2201 | 301, 701, 1101, … ,2301 |
|---|---|---|---|---|
| 星期： | 4 | 2 | 0 | 5 |
可以发现，每隔 4 个世纪，星期数就会重复一次。因为在数学上，-2 和 5 除以 7 的余数相同，所以我们不妨把这个重复序列中的 5 改为 -2。这样，4、2、0、-2 恰好构成了一个等差数列。利用等差公式，我们可以得到计算每个世纪第一年的 3 月 1 日星期数的公式：
其中，\(c\) 是世纪数减一。我们把公式（2）和公式（3）联立，代入特定的日期——3 月 1 日，可以得到：
利用同余定理，经过变换得到:
其中，\(\equiv\) 是表示同余的符号，括号中 \(\bmod 7\) 的意思是指 \(\equiv\) 两边的数除以 7 得到的余数相同。根据公式（4），我们可以知道在每个世纪的第一年，\((y-1) + \left[ \frac{y-1}{4} \right] - \left[ \frac{y-1}{100} \right] + \left[ \frac{y-1}{400} \right]\) 可以被 \(-2(c \bmod 4)\) 同余替换。进而计算 \(D\) 的公式得到如下形式：
注意！现在的计算公式只能适用于每个世纪的第一年。但是，有个这个公式，再加上计算一个世纪中闰日的部分，我们就可以很容易地得到计算这个世纪其他年份的日期的星期数的公式了。令 c 等于世纪数减一，y 等于世纪中的年份数（如 1994 年，则 c = 19，y = 94）。因为一个世纪中只有一百年，所以不用考虑“四百年又闰”的情况；因为每百年，即每个世纪最后一年的 y = 00，而 \(\left[ \frac{y}{4} \right]_{y=0} = 0\)，所以 \(\left[ \frac{y}{4} \right]\) 既可以计算四年一闰的情况，又满足百年不闰的要求 。综合这些情况，与得到公式（2）的过程类似，我们可以得到：
在公式（6）中，\(y\) 是年份的后两位。
最后，我们来把公式中的取模运算改成四则运算。设商为 \(q\)，余数为 \(r\)，则：
又因为，
可得：
代入公式（6）可得：
至此，我们就得到了蔡勒公式的最终形式。
