Project Euler 解題筆記 (3) – #167

來把題號擺在標題上好了。原本不放的原因是打算一篇可以擠好幾個，不過後來想想會能寫筆記的都應該不是三言兩語能解決的東西的。

#167 Investigating Ulam sequences (難度 75%)

所謂 Ulam sequence 是指這樣的數列：給定開頭兩個正整數 a < b，每次都找大於最後一項的整數當中，能夠唯一地表示為之前數列當中的兩項和的數中最小的那一個。題目給的例子是最簡單的 a=1, b=2 所產生的數列 (OEIS A002858)：1, 2, 3, 4, 6, 8, 11, ……。因為 5 = 1+4 = 2+3 有兩種兩項和表示法，所以不在數列中；因此 6 就在數列中了 (只有 2+4 一種)。令這樣的數列第 \(k\) 項記為 \(U(a,b)_k\)。

試求 \(\sum_{n=2}^{10} U(2,2n+1)_k\)，其中 \(k=10^{11}\)。

一般化的 Ulam sequence 其實有很多很詭異的性質，例如單就上面提的 a=1, b=2 產生的數列來說：

• 它是一個完全序列，表示任何正整數都可以表示成數個這個數列中的項的和；
• 除了開頭三項 1, 2, 3，其他任何連續三項都滿足三角不等式；
• 在數值方面它有一個很大略的規律：如果把所有正整數排成 216 個一排，點出在數列中的數，會發現它們大致地排成十條帶狀圖，帶與帶之間有個還不錯明顯的空白，也就是數值上大約有個差 21.6 左右的「週期關係」在；
• 如果只考慮數列本身的話，它也被人發現其傅立葉轉換當中有兩個異常突出的頻率：對這數列的前十億項，存在一個數 \(\theta\approx2.5714474995\)，使得這十億項當中除了四個都滿足 \(\cos(\theta a_n)<0\)。

等等。不過這些都無助於我們求出第 10¹¹ 項－－這個項數實在是太後面了，漸近性質沒辦法精準點到特定的某一項來。

這題的線索其實在上面連結到的 Ulam number 條目裡就有提到了：

A sequence of (u, v)-Ulam numbers is regular if the sequence of differences between consecutive numbers in the sequence is eventually periodic. When v is an odd number greater than three, the (2, v)-Ulam numbers are regular.

— 英文維基百科 Ulam number 條目

正好，題目要求的這個序列有個很特別的性質：如果開頭兩項是 2 和比 3 大的一個奇數時，相鄰兩項的差在一定項數之後會循環。因此我們只要能求出這是什麼樣的循環，就能直接推算出第 10¹¹ 項是什麼數了。

這是什麼樣的循環？就先簡單一點，看 \(n=2\) 的序列長怎樣：

2, 5, 7, 9, 11, 12, 13, 15, 19, 23, 27, 29, 35, 37, 41, 43, 45, 49, 51, 55, 61, 67, 69, 71, 79, 83, 85, 87, 89, 95, 99, ……

……好像除了 2 和 12 之外都是奇數?!

2 之後的下一個偶數是 12 很容易看得出來，但這之後為什麼會一個偶數都沒有就很令人費解了。不過如果假定這之後真的沒有其他偶數，那奇數的出現樣式就可以用這個敘述描述：當奇數 \(p\) 滿足 \(p-2\) 和 \(p-12\) 恰有一個在數列中，則 \(p\) 就在數列中。

因此由 5 到 15 連續六個奇數都在其中；5 和 15 都在所以 17 不在數列中；7 在但 17 不在所以 19 在；9 和 19 都在所以 21 不在；……。

或許是因為上一篇文章在談偽隨機亂數吧，這樣的結構其實立刻讓我想到只有在文章中提過名字的 LFSR 產生器。

LFSR

LFSR，全名 Linear Feedback Shift Register (線性反饋移位暫存器)，維基百科上寫說它是「給定前一狀態的輸出，將該輸出的線性函數再用作輸入的移位暫存器」。這個「暫存器」，它當中的位元進行「線性」函數組合，其結果「反饋」回到原暫存器中，將暫存器中的原值「移位」以放入新結果，所以叫做這個名字。最常見的線性函數是 XOR，也就是反饋的公式就是暫存器中的某些位元進行 XOR 後作為新結果移回暫存器中。維基百科上舉了一個 16 位元的例子，公式取第 11、13、14、16 位元進行 XOR 再移回去。

數學上，這樣的 LFSR 可以證明是 \(GF(2^n)\) (大小為 \(2^n\) 的有限體) 的操作。事情是這樣的：LFSR 其實有兩種定義型式，上面提的稱做 Fibonacci 式 (因為它和費氏數列的遞迴式很像，下一項產生使用到前面幾項的運算)；另一種稱做 Galois 式，它的定義是在移位時把移掉的輸出位元和對應位置的位元進行 XOR。如果對利用 \(\mathbb{Z}_2\) 係數多項式定義 \(GF(2^n)\) 的運算很熟的話，就能看得出來 Galois 式 LFSR 其實就是在這種運算之下「乘」以多項式 \(x\) 的運算 (這時所取的位元變成了有限體的除多項式)；它和有限體 (Galois Field) 之間的這種直接關連正是這種型式被稱做 Galois 式的原因。

這兩種定義型式之間有個很有趣的關連：它們所求出的輸出位元除了在時間上平移之外是一模一樣的序列。這是可以由內部狀態的對應方式來證明的，不過有點繁瑣就略過不談；但這個對應給了我們使用有限體理論求得周期的機會：有限體理論告訴我們，當除多項式 \(p(x)\) 是在 \(\mathbb{Z}_2[x]\) 中的不可約多項式時，乘以 \(x\) 這個操作 (連帶地，LFSR 的操作) 會有最大週期 \(2^n-1\)。這就是我們想找的週期性質。

以上面 2, 5 開始的序列為例，若把連續奇數有沒有在數列中以 1 和 0 描述的話，則這個數列正好就是以 \(x^6+x+1\) 描述的 LFSR 從 000001 (表示 5 在數列中, 但它前面的奇數都不在) 開始所產生的序列。由於 \(x^6+x+1\) 在 \(\mathbb{Z}_2[x]\) 中不可約，因此我們立刻就知道 LFSR 週期是 \(2^6-1=63\)。這代表，奇數 \(p\) 和 \(p+63\times2\) 要嘛都在數列中，要嘛都不在數列中。

實際列出來看看吧：它是 OEIS A007300。可以看到 5+126=131 到 141 也是連續六個奇數出現，而如果把它跟 5 開始的奇數併排的話：

  5,   7,   9,  11,  13,  15,  19,  23,  27,  29,  35,  37,  41,  43,  45,  49,
      51,  55,  61,  67,  69,  71,  79,  83,  85,  87,  89,  95,  99, 107, 109, 119,
131, 133, 135, 137, 139, 141, 145, 149, 153, 155, 161, 163, 167, 169, 171, 175,
     177, 181, 187, 193, 195, 197, 205, 209, 211, 213, 215, 221, 225, 233, 235, 245,
...

上面這個併排可以看到前 32 個奇數加上 126 正好就是接下來 32 個奇數。這就是上面引文提到的「regular」性質。

這個性質其實還可以回頭說明為什麼偶數在 12 之後就都不會出現了：任意選大於 15 的六個連續奇數 \(p, p+2, \dots, p+10\)。如果這六個連續奇數裡有至少兩個在數列中，那它們和 5~15 之間就會出現至少兩個 \(p+15\) 的組合，因此不在數列中。如果只有一個，它只會在跨週期附近出現 (例如 119 附近的 111~121 或 121~131)；對應要找的偶數是 126、128、……、136，這些偶數可以被 107 或 109 抓住另一個組合：

• 126 = 119+7 = 107+19
• 128 = 119+9 = 109+19
• 130 = 119+11 = 107+23
• 132 = 119+13 = 109+23
• 134 = 119+15 = 107+27
• 136 = 131+5 = 109+27

因此這些偶數也都不會在數列當中。由於奇數有這樣的週期性，同樣的論證可以適用在其他週期上，因此就能得到結論：大於 12 的偶數通通不會出現在數列當中。¹

那麼，我們就能從這裡求出 \(U(2,5)_{10^{11}}\) 了：這 \(10^{11}\) 項中，扣掉兩個偶數，剩下的每 32 個一組。做除法得到 \((10^{11}-2)\div32=3,124,999,999\cdots30\)，即是所求為第 3,124,999,999+1=3,125,000,000 組，當中的第 30 個數。第一組的第 30 個數是 107，所以第 3,125,000,000 組的第 30 個數就是
$$U(2,5)_{10^{11}}=107+3,124,999,999\times126=393,749,999,981$$

可約多項式

嗯，鑑於不捏他原則，上面我既然寫出部份答案出來了就表示這題其實還有其他問題。

來看下一組 \(U(2,7)\) 吧。類似的方式可以發現它是以 \(x^8+x+1\) 描述的 LFSR (因為唯二的偶數是 2 和 16)，同樣由 00000001 開始 (表示 7 在數列中)。但是問題來了，在 \(\mathbb{Z}_2[x]\) 當中有：
$$x^8+x+1=(x^2+x+1)(x^6+x^5+x^3+x^2+1)$$
因此它的週期並不是 \(2^8-1=127\)。

那這要怎麼求週期？這裡給我們線索的是中國剩餘定理。還記得上面我們說過 Fibonacci LFSR 的操作可以對應到 Galois LFSR (及其對應的有限體中) 乘以 \(x\) 的操作嗎？如果除法群的除式 \(p(x)\) 可分解但沒有因式重覆的話，那除法群中除以 \(p(x)\) 取餘的操作其實可以分別對其因式取餘，然後再由中國剩餘定理合併起來。這樣一來，乘以 \(x\) 再除以 \(p(x)\) 取餘操作的循環就能夠「分解」成對因式取餘的循環，再合併起來。

以 \(x^8+x+1=(x^2+x+1)(x^6+x^5+x^3+x^2+1)\) 為例，這代表它的週期可以拆成除以 \(x^2+x+1\) 和除以 \(x^6+x^5+x^3+x^2+1\) 的兩個週期；這兩個都是不可約多項式了，所以可以套已知結論，前者週期 \(2^2-1=3\)，後者週期 \(2^6-1=63\)，因此全部的週期就是這兩者的最小公倍數，也就是 63。

那麼我們就只要對所有 \(x^{2n+2}+x+1\) 分解，抓出不可約因式的最高次求個別週期 (好在因式都沒有重覆的)，就能求得 LFSR 週期了。知道 LFSR 週期代表每一段差是已知，那就只要寫支程式算過去，求出這一段數字內有多少數有選中就是數列週期；有了數列週期，再用和上面一樣的方式就能求得所要求的第 \(10^{11}\) 項了。

提件稍微值得一提的事當結尾：只要上面拆分因式後的各週期沒有公因數，所求出來的第 \(10^{11}\) 項的值都會約略在 \(4\times10^{11}\) 左右。這其實是可以解釋的：在單一因式的全週期當中，0 會比 1 少一個 (因為這個全週期只有多項式 0 無法獲得)，也就是 1 會大約有一半，那要取到 \(10^{11}\) 個奇數，就要掃過大約 \(2\times10^{11}\) 個奇數，自然約略在 \(4\times10^{11}\) 附近了。多個週期如果沒有公因數，那所有週期的所有組合都會跑過一次，因此 1 大約有一半的這個估計仍然成立；但如果有公因數則就會有沒出現過的組合，總體來說 1 的個數就會少於週期的一半，因此結果就會比 \(4\times10^{11}\) 大很多了。

註腳

1. 對證明細節有些敏銳的人可能注意到我這裡的寫法其實有點怪怪的。理論上，正確的證明會需要用數學歸納法，假設這兩個性質 (偶數只有 2 和 12、奇數遵循這 32 個數一「循環」的方式) 對某個數以下成立，然後用其推出大一些些的範圍成立，最後用數學歸納法得到這兩個性質全部成立；一些細節像是如 119 這種「孤立」的數，也能證明它前兩項一定會是連續兩個奇數，進而完成證明。這個推論邏輯認真看的話會發現它其實是交錯式的：奇數性質推出大一些範圍的偶數性質、偶數性質推出大一點的奇數性質。要用這種行文方式的話會讓整串文章被這種交錯式的邏輯給卡住，還不如先講完比較麻煩的奇數性質，再回頭把比較簡單的偶數性質補完。 ↩︎