0. 本文的初衷及蔡勒公式的用處
前一段時間，我在準備北郵計算機考研復試的時候，做瞭幾道與日期計算相關的題目，在這個過程中我接觸到瞭蔡勒公式。先簡單的介紹一下蔡勒公式是幹什麼用的。
我們有時候會遇到這樣的問題：看到一個日期想知道這一天是星期幾，甚至看到一個歷史日期或紀念日，我們想快速的知道這一天是星期幾。對於這個問題，如果用編程的方式，應該怎麼實現呢？你可能已經有思路瞭，比如你知道某個日期是星期幾，把這個日期作為原點，然後計算目標日期和這個原點之間相差多少天，再除以 7 求餘數，最後通過餘數判斷目標日期的星期數。通過這樣的過程，你確實可以得到正確的結果，但這不夠快速也不夠優雅。對於這個問題，如果你懂得蔡勒公式，那就變得異常簡單瞭，你隻需要將年月日等數據代入公式，然後計算出結果，這一結果就是目標日期對應的星期數。
當我知道蔡勒公式之後，我覺得它非常有趣也很酷，所以我不僅希望掌握公式的使用，也希望可以掌握公式背後的推導過程。然而，當我在網上搜索相關的文章時，我發現幾乎所有向我展示的博客（從零幾年到最近的 19 年）大多是轉載、復制於這篇文章：
這篇文章質量很不錯，講解過程自然流暢，但是在一些細節上存在錯誤，有些推導步驟讓人感到困惑。因此，當我掌握蔡勒公式後，很希望可以將我的理解輸出出來，讓想要學習蔡勒公式推導過程的人看到一些新的材料。好瞭，廢話少說，我們開始吧。
1. 蔡勒公式的形式
如果你對公式的推導過程不感興趣，隻是希望使用蔡勒公式，那麼隻看此小節即可。蔡勒公式的形式如下：
begin{align} D &= left[ frac{c}{4} right] – 2c + y + left[ frac{y}{4} right] + left[ frac{13(m+1)}{5} right] + d – 1 \[2ex] W &= D bmod 7 end{align} \
其中：
- 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，帶入公式：
begin{align} D &= left[ frac{19}{4} right] – 2 times 19 + 94 + left[ frac{94}{4} right] + left[ frac{13 times (12 + 1)}{5} right] + 13 – 1 \[2ex] &= 4 – 38 + 94 + 23 + 33 + 13 – 1 \[2ex] &= 128 \[2ex] W &= 128 bmod 7 = 2 end{align}\
由此可得 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 瞭：
D_2 = (y-1) times 365 + left[ frac{y-1}{4} right] – left[ frac{y-1}{100} right] + left[ frac{y-1}{400} right]\
簡單的解釋一下。因為一年最少有 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 ，並且我們可以得到這樣一張表格：
仔細觀察，可以發現 Excess 從 3 月份開始是 3、2、3、2、3 的循環，因此，當 mgeq 3 時，Accum(m) 的值的增幅也是 3、2、3、2、3 的循環。因為每 5 個月增加 13，所以把 frac{13}{5} 作為系數；因為 Accum(m) 的值是離散的（都是整數），所以我們用取整運算符，得到：
f(x) = left[ frac{13}{5}x right] \
我們將 x 的值取 1，2，3……，然後觀察 f(x) 的值，可得下面這張表格：
我們可以發現，當 x geq 4 時，f(x) 的值的增幅也是 3，2，3，2，3 的循環。也就是說 f(x) 的值的增幅（x geq 4）與 Accum(m) 的值的增幅（m geq 3）相同，這意味著 f(x) 與 Accum(m) 之間相差一個常數，我們隨便帶入一個具體的值計算：
f(4) – Accum(3) = 10 – 3 = 7\
可知相差的常數為 7。由此可得，當 m geq 3 時，Accum(m) 的值的序列，等於當 x geq 4 時， f(x) – 7 的值的序列。這樣我們就得到瞭 Accum(m), m geq 3 的函數形式：
begin{align} Accum(m) &= f(m+1) – 7 \[2ex] &= left [ frac{13(m+1)}{5} right ] – 7 end{align}\
進一步，我們可以得到 D_3 的函數形式：
D_3 = begin{cases} d, &m = 1 \[2ex] 31 + d, &m = 2 \[2ex] (m-1) times 28 + left[ frac{13(m+1)}{5} right] – 7 + d + i, &m geq 3 end{cases} \
其中，平年時，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_3 = (m-1) times 28 + left[ frac{13(m+1)}{5} right] – 7 + d, quad 3 leq m leq 14\
現在，我們已經得到瞭 D_2 和 D_3 的計算函數，由 D = D_2 + D_3，可知：
D = (y-1) times 365 + left[ frac{y-1}{4} right] – left[ frac{y-1}{100} right] + left[ frac{y-1}{400} right] + (m-1) times 28 + left[ frac{13(m+1)}{5} right] – 7 + d, quad 3 leq m leq 14\
註意！這個公式離正確形式還差一小步。因為在當前的公式中，每年的 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 = (y-1) times 365 + left[ frac{y}{4} right] – left[ frac{y}{100} right] + left[ frac{y}{400} right] + (m-1) times 28 + left[ frac{13(m+1)}{5} right] – 7 + d, quad 3 leq m leq 14 tag{1}
計算出 D 的值後，對 7 取模即可得到星期數。
根據同餘定理，D 除以 7 所得的餘數等於 D 式的各項分別除以 7 所得餘數之和（餘數之和大於等於 7 時，再對 7 取餘），因此可以消去 D 式中能被 7 整除的項，進行化簡：
begin{align*} D &= (y-1) times 365 + left[ frac{y}{4} right] – left[ frac{y}{100} right] + left[ frac{y}{400} right] + (m-1) times 28 + left[ frac{13(m+1)}{5} right] – 7 + d \[2ex] &= (y-1) times (364+1) + left[ frac{y}{4} right] – left[ frac{y}{100} right] + left[ frac{y}{400} right] + left[ frac{13(m+1)}{5} right] + d \[2ex] &= (y-1) + left[ frac{y}{4} right] – left[ frac{y}{100} right] + left[ frac{y}{400} right] + left[ frac{13(m+1)}{5} right] + d tag{2} end{align*} \
簡單說明一下：
begin{align*} (y-1) times 365 &= (y-1) times (364+1) \[2ex] &= (y-1) times 364 + (y-1) \[2ex] &= (y-1) times 52 times 7 + (y-1) end{align*} \
顯然，結果中的第一項是 7 的倍數，除以 7 餘數為 0，因此 (y-1) times 365 除以 7 的餘數其實就等於 (y-1) 除以 7 的餘數，我們隻保留 (y-1) 就夠瞭。化簡過程中，其他被銷去的項同理。
公式（2）還不是最簡練的形式，我們還可以對年份進行處理。我們現在用公式（2）計算出每個世紀第一年 3 月 1 日的星期數，得到如下結果：
可以發現，每隔 4 個世紀，星期數就會重復一次。因為在數學上，-2 和 5 除以 7 的餘數相同，所以我們不妨把這個重復序列中的 5 改為 -2。這樣，4、2、0、-2 恰好構成瞭一個等差數列。利用等差公式，我們可以得到計算每個世紀第一年的 3 月 1 日星期數的公式：
W = 4 – 2(c bmod 4) tag{3}
其中，c 是世紀數減一。我們把公式（2）和公式（3）聯立，代入特定的日期——3 月 1 日，可以得到：
((y-1) + left[ frac{y-1}{4} right] – left[ frac{y-1}{100} right] + left[ frac{y-1}{400} right] + 11) bmod 7 = 4 – 2(c bmod 4)\
利用同餘定理，經過變換得到：
(y-1) + left[ frac{y-1}{4} right] – left[ frac{y-1}{100} right] + left[ frac{y-1}{400} right] equiv -2(c bmod 4) quad (bmod 7) tag{4}
其中，equiv 是表示同餘的符號，括號中 bmod 7 的意思是指 equiv 兩邊的數除以 7 得到的餘數相同。根據公式（4），我們可以知道在每個世紀的第一年，equiv 可以被 -2(c bmod 4) 同餘替換。進而計算 D 的公式得到如下形式：
D = -2(c bmod 4) + left[ frac{13(m+1)}{5} right] + d tag{5}
註意！現在的計算公式隻能適用於每個世紀的第一年。但是，有個這個公式，再加上計算一個世紀中閏日的部分，我們就可以很容易地得到計算這個世紀其他年份的日期的星期數的公式瞭。令 c 等於世紀數減一，y 等於世紀中的年份數（如 1994 年，則 c = 19，y = 94）。因為一個世紀中隻有一百年，所以不用考慮“四百年又閏”的情況；因為每百年，即每個世紀最後一年的 y = 00，而 left[ frac{y}{4} right]_{y=0} = 0，所以 left[ frac{y}{4} right] 既可以計算四年一閏的情況，又滿足百年不閏的要求 。綜合這些情況，與得到公式（2）的過程類似，我們可以得到：
D = -2(c bmod 4) + (y-1) + left[ frac{y}{4} right] + left[ frac{13(m+1)}{5} right] + d tag{6}
在公式（6）中， y 是年份的後兩位。
最後，我們來把公式中的取模運算改成四則運算。設商為 q，餘數為 r，則：
4q + r = c \
又因為，
begin{align*} q &= left[ frac{c}{4} right] \[2ex] r &= c bmod 4 end{align*}\
可得：
c bmod 4 = c – 4 times left[ frac{c}{4} right]\
代入公式（6）可得：
D = left[ frac{c}{4} right] – 2c + y – 1 + left[ frac{y-1}{4} right] + left[ frac{13(m+1)}{5} right] + d tag{7}
至此，我們就得到瞭蔡勒公式的最終形式。