Project Euler 解題筆記 (2)

#152 Writing 1/2 as a sum of inverse squares (難度 65%)

這題躺在我的未解題裡躺很久了：有一陣子我在嘗試從頭開始答題，寫一寫就撞上這題這個軟牆，所以就開始往後亂跳題目 XD 也就是說，除了在約三四百題那附近時有比較集中在一出題就去解之外，其他的解題進度都是這樣散散的，看到有點難的題目就給它拼命跳，所以這題才留著這麼久。

題目很簡單：就如同標題說的，把 1/2 寫成不同整數平方的倒數和。題目敘述給了一個範例：

$$\begin{align}\dfrac{1}{2} &= \dfrac{1}{2^2} + \dfrac{1}{3^2} + \dfrac{1}{4^2} + \dfrac{1}{5^2} +\\
&\quad \dfrac{1}{7^2} + \dfrac{1}{12^2} + \dfrac{1}{15^2} + \dfrac{1}{20^2} +\\
&\quad \dfrac{1}{28^2} + \dfrac{1}{35^2}\end{align}$$

並且給出如果只考慮 2 到 45 的平方倒數的話，連同這個範例一共只有三組解。題目要問的是：如果放大範圍到 2 到 80 的話有幾組解？

之所以是軟牆是因為，2 到 80 這個範圍很微妙的稍微大了一點；這個題目本質上是一個子集和問題，要我們在給定的 79 個元素中選出數個來加到目標 1/2；由於子集和問題是著名的 NP-Complete 問題，所需要的工作量不會比把 79 個數個別選不選所得到的 \(2^{79}\) 種選項逐一嘗試來得好多少。

如果使用這些數字的數學性質呢？自然數平方的倒數和也是個很有名的問題，叫做巴塞爾問題，其總和約是 1.644934，也就是由 2 到無窮的平方倒數和是這和減去 1 約是 0.644934，和我們的目標 0.5 的差距又離的有一點遠，不太能夠提前過濾選項；能簡單推得的只有 \(\frac1{2^2}\) 不能不選，以及 \(\frac1{3^2}\) 跟 \(\frac1{4^2}\) 不能一起不選；要再下去就不容易推條件了。似乎還是需要深挖一點其他的數學條件來使用。

突破口：加起來吧！

說起來，題目給的這組答案很有趣：中間甚至出現了 28 和 35 這種 7 的倍數分母，但全部加起來之後卻能通通消光光，最後分母只剩下一個 2。到底這是怎麼消的？

容易看到這組解裡，分母是 7 的倍數的有三個：\(\frac1{7^2},\ \frac1{28^2},\ \frac1{35^2}\)，而它們加起來是：

$$\frac1{7^2}+\frac1{28^2}+\frac1{35^2}=\frac{20^2+5^2+4^2}{140^2}=\frac{441}{140^2}=\frac{9\times7^2}{140^2}=\frac9{20^2}$$

於是我們發現了：原來是因為通分相加之後，分子出現了 \(7^2\) 的倍數，因此可以把通分後分母的 \(7^2\) 給消掉。這個條件其實很大地限制了一些大因數的出場機會：例如 31，在 2 到 80 的範圍內的 31 的倍數只有 31 和 62，平方倒數和不會加出 \(31^2\) 的倍數出來，因此如果它們一出場，分母的 \(31^2\) 沒有機會消掉，結果也就不可能是 1/2 了。

照理來說，對奇質數來說所選的該奇質數的倍數之平方倒數和都要有這個性質，但這樣一來合數到底要給誰算？如果要歸給最大質因數的話，那小一點的質因數就有可能會加入一些和其他人長得不太一樣的值，這樣搜尋也很麻煩。以上面為例，如果接下來考慮 5 的倍數的話，這裡將會是：

$$\frac1{5^2}+\frac1{15^2}+\frac1{20^2}+\frac9{20^2}=\frac{12^2+4^2+3^2+9\times3^2}{60^2}=\frac{250}{60^2}=\frac{10}{12^2}=\frac5{72}$$

這個 \(\frac5{72}\) 就不像上面的結果一樣是個有理數完全平方。雖然 5 的因數消掉了，但因為形式變化太大，這並不容易混進 3 的一般討論裡求算。

大質因數搜尋

不過，至少對像上面舉的 31 這種大一點的質因數我們可以利用這個性質加以排除。很自然的問題就是還有沒有一些像這樣的組合？這樣一來我們只需要對只有小質因數的數 (這種數被稱為光滑數) 去搜尋即可。結果很有趣：如果考慮最大質因數大於等於 11 的數之組合的話，在題目問的 2 到 80 的範圍中只有一個組合有我們要的性質：

$$\frac1{13^2}+\frac1{39^2}+\frac1{52^2}=\frac{12^2+4^2+3^2}{156^2}=\frac{169}{156^2}=\frac1{12^2}$$

這個分母甚至連 5 和 7 都沒有，所以我們可以進一步搜尋最大因數是 7 的組合，不會受這一組到底選不選對搜尋的影響。以下就是範圍內能夠使用的 7 的組合，一共 19 組：

這裡的 19 組的分母很多都有帶數個 5，所以不容易代入 5 的倍數討論。不過到這裡至少我們篩到剩下所有的 5-光滑數了；小於等於 80 的 5-光滑數只有 29 個，比起原來的 79 個少太多了。

再深入思考就會發現，我們這樣做其實已經把一個 79 元素的子集和問題，轉變成了 40 個 29 元素的子集和問題。40 個是因為，7 的選項可以選 19 種組合或不選，然後上面 13 的那組也能選和不選，總計就有 40 個目標；由於這 40 個目標在問的母集合是同一個 29 元素集合，我們可以用同樣的動態規劃解法全部處理一遍之後，詢問這 40 個目標有沒有人有踩中就行了。分數的處理則可以簡單地把所有要處理的分數全部通分成這 29 個數的最小公倍數 \((2^6\cdot3^3\cdot5^2=64\cdot27\cdot25=43200)\) 的平方即可：

/* Final target:

1/2 = 933120000/43200^2

*/

const int64_t half = 933120000;

/* 7 and 13 terms:

{21, 42, 63}						 5760000/43200^2
{21, 35, 56}						 6350400/43200^2
{28, 35, 42, 56, 63, 70}			 6408000/43200^2
{14, 56, 70}						10497600/43200^2
{14, 35, 56, 63}					12110400/43200^2
{14, 28, 42}						12960000/43200^2
{14, 21, 35, 63, 70}				16128000/43200^2
{14, 21, 42, 56, 63, 70}			16257600/43200^2
{14, 21, 28, 35, 42, 56}			19310400/43200^2
{7, 63, 70}							38937600/43200^2
{7, 28, 35}							41990400/43200^2
{7, 28, 42, 56}						42120000/43200^2
{7, 21, 35, 56, 63, 70}				45288000/43200^2
{7, 21, 28, 42, 70}					46137600/43200^2
{7, 21, 28, 35, 42, 63}				47750400/43200^2
{7, 14, 21}							51840000/43200^2
{7, 14, 28, 42, 63, 70}				51897600/43200^2
{7, 14, 28, 35, 56, 70}				52488000/43200^2
{7, 14, 21, 28, 35, 42, 56, 63, 70}	58248000/43200^2
{13, 39, 52}						12960000/43200^2

*/

const int64_t sevens[] = {
       0,
 5760000,  6350400,  6408000, 10497600, 12110400, 12960000, 
16128000, 16257600, 19310400, 38937600, 41990400, 42120000, 
45288000, 46137600, 47750400, 51840000, 51897600, 52488000, 
58248000
};

constexpr int SEVENS_COUNT = _countof(sevens);

const int64_t thirteens[] = {0, 12960000};

constexpr int THIRTEENS_COUNT = _countof(thirteens);

/* 5-smooth numbers:

1/80^2 =   291600/43200^2
1/75^2 =   331776/43200^2
1/72^2 =   360000/43200^2
1/64^2 =   455625/43200^2
1/60^2 =   518400/43200^2
1/54^2 =   640000/43200^2
1/50^2 =   746496/43200^2
1/48^2 =   810000/43200^2
1/45^2 =   921600/43200^2
1/40^2 =  1166400/43200^2
1/36^2 =  1440000/43200^2
1/32^2 =  1822500/43200^2
1/30^2 =  2073600/43200^2
1/27^2 =  2560000/43200^2
1/25^2 =  2985984/43200^2
1/24^2 =  3240000/43200^2
1/20^2 =  4665600/43200^2
1/18^2 =  5760000/43200^2
1/16^2 =  7290000/43200^2
1/15^2 =  8294400/43200^2
1/12^2 = 12960000/43200^2
1/10^2 = 18662400/43200^2
1/9^2 =  23040000/43200^2
1/8^2 =  29160000/43200^2
1/6^2 =  51840000/43200^2
1/5^2 =  74649600/43200^2
1/4^2 = 116640000/43200^2
1/3^2 = 207360000/43200^2
1/2^2 = 466560000/43200^2

*/

const int64_t fivesmooth[] = {
   291600,    331776,    360000,    455625,    518400,    640000, 
   746496,    810000,    921600,   1166400,   1440000,   1822500, 
  2073600,   2560000,   2985984,   3240000,   4665600,   5760000, 
  7290000,   8294400,  12960000,  18662400,  23040000,  29160000, 
 51840000,  74649600, 116640000, 207360000, 466560000
};

當然實際處理並沒有真的開一個長度為 \(43200^2/2\) 的陣列在那裡一直算 0，而是就只記 DP 陣列裡哪個位置有值再去推進，但概念是一樣的。這裡還能因為我們要問的目標最大就是 1/2，所以所有超過它的總和都可以忽略。

這樣寫出來的程式在我的電腦上差不多踩在 Project Euler 的「一分鐘規則」線上，所以應該算是勉強過關的程式碼了。稍加改寫，花久一點的時間是可以跑出所有的解個別是什麼，這裡就列出範圍內兩組選擇數目最多的解 (選了 24 項) 做為結束吧：

$$\begin{align}\dfrac{1}{2} &= \dfrac{1}{2^2} + \dfrac{1}{3^2} + \dfrac{1}{5^2} + \dfrac{1}{6^2} + \dfrac{1}{7^2} + \dfrac{1}{8^2} + \dfrac{1}{10^2} + \dfrac{1}{13^2} + \dfrac{1}{14^2} + \dfrac{1}{18^2} + \dfrac{1}{21^2} + \dfrac{1}{24^2} + \\
&\quad \dfrac{1}{28^2} + \dfrac{1}{30^2} + \dfrac{1}{35^2} + \dfrac{1}{39^2} + \dfrac{1}{40^2} + \dfrac{1}{42^2} + \dfrac{1}{45^2} + \dfrac{1}{52^2} + \dfrac{1}{56^2} + \dfrac{1}{60^2} + \dfrac{1}{63^2} + \dfrac{1}{70^2}\\
&= \dfrac{1}{2^2} + \dfrac{1}{3^2} + \dfrac{1}{5^2} + \dfrac{1}{6^2} + \dfrac{1}{7^2} + \dfrac{1}{9^2} + \dfrac{1}{10^2} + \dfrac{1}{13^2} + \dfrac{1}{14^2} + \dfrac{1}{15^2} + \dfrac{1}{20^2} + \dfrac{1}{21^2} + \\
&\quad \dfrac{1}{24^2} + \dfrac{1}{28^2} + \dfrac{1}{30^2} + \dfrac{1}{35^2} + \dfrac{1}{39^2} + \dfrac{1}{40^2} + \dfrac{1}{42^2} + \dfrac{1}{52^2} + \dfrac{1}{56^2} + \dfrac{1}{63^2} + \dfrac{1}{70^2} + \dfrac{1}{72^2}
\end{align}$$