第一篇認真的文章就來記錄一下這個吧。
題目
原題目在這裡。這篇要記的是中間和另一個人有些討論的第三題：
觀察
簡單觀察一下題目所指定的這些操作：
- 平方根後除非是恰好為整數，否則無法馬上接階乘。
- 因此在連續的平方根操作後幾乎都是接向下取整，而在這樣的平方根當中插入向下取整不影響後一個向下取整的結果 (理由是向下取整函數值不會比輸入值大，以及整數的平方仍是整數。) ，因此這樣做不會是最佳。
- 除了 \(1!=1\) 之外，其他整數的階乘都不是完全平方數，因此開平方根不是整數，從而階乘後連續的平方根必須接上一個向下取整才可再接階乘。
- 如果中間結果出現 1 或 2 則後續就無法脫離，這是因為唯一可能將數字變大的操作是階乘，但 \(1!=1\) 和 \(2!=2\) 都沒有變大。
因此除了一開始的 2021 之外，後續的動作可以由階乘操作切開，每一段都會是〈階乘→多次平方根→向下取整〉這樣子的段落。2021 因為可以不先階乘，所以是獨立的〈2021→零或多次平方根→有平方根的話向下取整〉這樣的開頭段落。
階乘運算
這樣切開的好處是：開頭段落不看的話 (反正那邊數字夠小)，每一段只有一個階乘運算，而我們只需要在最後向下取整轉成整數時算出實際的整數值，階乘和平方根可以利用對數簡化計算。
不過問題來了：這樣一來我們必須計算 \(\log n!\) 的值，而且需要稍微精確一些，不能只是由積分推得的粗略估計值：$$\log n!=\sum_{i=1}^n \log i\leq\int_1^{n+1} \log x dx=x\log x-x\bigg\vert_1^{n+1}=(n+1)\log(n+1)-n$$ 對 \(n=100\) 時，\(\log 100!\approx363.739\)，但代入積分式得到 \(101\log101-100\approx366.127\)，自然對數值差了約 2.4，也就是實際數值差了約十倍……
就算用下和跟上和的平均估計，積分式代 \(n=99\) 得到約 361.517，兩者平均是 363.822，依然大了約 0.082，實際數值約大了 8.6%；這顯然不足以使中間計算得到正確結果。
搜尋了一下網路，在 mathstackexchange 上找到了這一篇回答，給出了一個拉馬努金提出過的近似式：$$\log n!\approx n\log n-n+\frac16\log\left(n(1+4n(1+2n))+\frac1{30}\right)+\frac12\log\pi$$ (加 \(\frac1{30}\) 是在其他回答裡有提到過的修正項，但在這題裡幾乎沒有影響)
這個近似式有多好？簡單試了一下在 \(n=100\) 時，給出的答案和正確值到小數下第 10 位都相同！只能說不愧是拉馬努金……
演算法
實際去搜尋的演算法則沒什麼特別的，使用的是經典的 Dijkstra 演算法，但並不直接把整張圖建立起來。原因很簡單：階乘的數字太大了。上面提到過 \(\log100!\approx363.739\)，這代表 \(100!\) 實際的數值約是 \(e^{363.739}\approx 9.332\times10^{157}\)，簡單的程式不用科學記號 (aka. 浮點數) 根本無法計算。
這個太大造成的影響之一是很多大數字其實在搜尋過程中幾乎踩不到。繼續以 \(100!\) 為例，由它開始連續開平方根會依序得到：
- \(100!\approx9.332\times10^{157}\)
- \((100!)^{1/2}\approx9.661\times10^{78}\)
- \((100!)^{1/4}\approx3.108\times10^{39}\)
- \((100!)^{1/8}\approx5.575\times10^{19}\)
- \((100!)^{1/16}\approx7466640272.6695\)
- \((100!)^{1/32}\approx86409.7232\)
- \((100!)^{1/64}\approx293.9553\)
- \((100!)^{1/128}\approx17.1451\)
- \((100!)^{1/256}\approx4.14066\)
所以如果這一步是 100，下一步可能會踩到 7466640272、86409、293、17、4；對小數字來說其他起點有可能會得到 16 或 18 或 292 或 294，但 86408、86410、甚至是 7466640271 就很難說會不會馬上就算得出來。這就是之所以為什麼不直接把整張圖建立起來的原因；在 relax 的時候再去求出從現在這一點出發走幾步會到哪些地方就行了。
大數上限
數字太大的影響還有一個：這種大數字的階乘只會是個更大一團東西。繼續上面的例子，如果把 \(100!\) 代進 \(\log n!\) 的式子裡，會求得一個結果約為 \(3.38\times10^{160}\)；注意到這是對數值，所以取指數會得到實際值是 \((100!)!\approx10^{1.47\times10^{160}}\)，已經出現指數塔了。不過好在這個問題影響不大：由於操作中讓數字減少的方式主要是平方根 (向下取整是冪等運算，只有一次有用；而由前述討論也不會取整之後再平方根，所以作用更有限)，而平方根一次只會讓指數對半，所以這種大數要減少到小數範圍的平方根次數也會很多－－例如 \((100!)!\) 需要平方根約 530 次才會變小到相對能掌握的範圍：\(1.47\times10^{160}\big/2^{530}\approx4.18\)，也就是 \((100!)!\) 開 530 次平方根的值約 \(10^{4.18}\approx15000\)。在已知有 \(K\) 步解的時候，如果出現這種其階乘需要更多次平方根才會減少到小範圍的數的話，對其用階乘是沒有意義的；也就是說，我們可以設定一個上限，如果中間結果超過這個數字就不對它使用階乘操作。
令這個上限是 \(M\)。於是這個很粗略的條件就是 \(M\) 的階乘要平方根 \(K\) 次以上才能回到比 \(M\) 小的數：$$\frac{\log(M!)}{\log M}=2^K$$ 若以 \(M\log M\)做為分子的近似，左邊就化簡為 \(M\)，也就是說所要的上限就約是 \(M\approx2^K\)。
(如果代幾個實際的 \(M=2^k\) 值進去會發現求出來的 \(K\) 會比 \(k\) 稍小，這由上面的近似公式即可看到原因：因為有一個 \(-n\) 的項。實際運用上為了補償這個差值，是可以在選定 \(K\) 之後將上限定為比 \(2^K\) 稍大一些些的數；因為這個差約 \(O(\frac{n}{\log n})\)，我們對指數 \(K\) 加個小值差不多就能涵蓋得到)
前一版程式
說了這麼多，還是把已經跑過的這前一版程式給貼出來好了：連結在此。它的執行結果是：
Found 110 at 88 steps: 110 <- 784453!, sqrt * 21 784453 <- 24427!, sqrt * 14 24427 <- 496!, sqrt * 8 496 <- 37!, sqrt * 4 37 <- 4062!, sqrt * 13 4062 <- 46!, sqrt * 4 46 <- 26!, sqrt * 4 26 <- 6!, sqrt * 1 6 <- 2021, sqrt * 2
這就是在原文推文裡我推的 88 步解。列出來的順序反過來的原因就只是 Dijkstra 演算法做出來的路徑本來就是倒過來由目標走回來的。
這前一版程式裡很多地方其實是因為只是簡單寫寫而已所以任意決定的，例如計算上限 MAX_LOG 原本是因為我的對數是用常用對數來算所以是正好切在 \(10^{15}\) 的，但後來想想這裡的對數果然還是用自然對數比較不會有意外，所以就改成了 \(e^{35}\approx1.58\times10^{15}\approx2^{50.5}\)，但還遠不到這裡給出來需要的 \(2^{88}\) 上限。
然而說到精確度，扣掉比較沒有資料的近似公式精確度之外，這裡其實還有一個我有點擔心的地方：中間的 \(\frac16\log\left(n(1+2n(1+4n)+\frac1{30}\right)\) 的計算，參數約 \(8n^3\) 的數值我直接用 double 計算了，因此數字可能會大到約 \(3\times10^{46}\)－－直覺上來說因為是用在對數，這裡的精確度損失應該影響不大，但究竟有沒有影響還真的不敢說…… (也就是因為這裡會這麼大的關係，上面列公式時才會說這裡有沒有 \(\frac1{30}\) 幾乎沒有影響。我是加進去了就是。)
近似公式的精確度我有用過 Mathematica 任意選了幾個 \(10^{14}\) 等級左右的數去試了，意外的精準：結果大概至少都有十幾位是正確的。也就是說，只要這公式有這種等級的準確，那上面的精確度損失的影響就真的不大。
其他部份就是很傳統的 Dijkstra 演算法，只不過因為用的是沒有 Decrease Key 的 std::priority_queue 的關係，用了一個也很常見的做法：原本 Decrease Key 的操作也照樣塞進 queue 裡，在彈出來時檢查這一項是不是已經做過了，是就忽略它。
LOG_ERR 是用來印出很大量的中間結果，這樣我能檢查到底中間的計算過程中有沒有什麼問題進而修正。輸出的 log 量大約在 2~3 MB 左右，所以不印 log 的話結果幾乎是馬上就出來了。
溢位 bug
於是這是某次我回頭檢查 log 時看到的片段：
可以看到有一個紀錄雖然值出來了但沒有後續運算。在這版程式的再前一版曾經有過類似的紀錄，後來查出來是 \(\log n!\) 的當中 \(n\) 的三次式計算溢位造成 log 吃到負數給出 NaN；這次一樣的現象所以回頭再查一次那段程式，結果還真的找到又一個溢位 bug：
在程式其他地方是 64-bit 整數的數變成 32-bit 傳進來了……這造成所有大數字的計算都不是正確的，等同於是把階乘操作的上限降到 \(2^{31}\) 了。
把它改掉，以及把小數字的階乘提出去不經過 log (避免一來一回連小數字都會歪掉)，改版後的程式在這裡。執行結果是：
Found 110 at 75 steps: 110 <- 22280837523899!, sqrt * 47 22280837523899 <- 1278!, sqrt * 8 1278 <- 120!, sqrt * 6 120 <- 5!, sqrt * 0 5 <- 6!, sqrt * 2 6 <- 2021, sqrt * 2
這就是我在後來的回文裡給出的 75 步解了。這裡出現了一個很可怕的 14 位數的階乘開 47 次平方根，還很精準的直接踩中目標 110 (大約值是 110.77)，所以有點怕怕的才會想要多方驗算，不過看起來這個計算確實足夠精確……
會有這種大數字直接踩中目標的狀況也就是我在開頭文章的推文說的，不確定允許中間過程是大數字有沒有可能找到更少步數解，哪天哪個大數字天外飛來一筆直接撞上目標或接近目標路線的數都有可能；這只能在我們找到至少一個解之後利用前一段提到的 \(M\) 值估計才能肯定地排除。
最後……
在我的回文那篇底下和我討論的那位版友用了另一個計算架構 Magma 據說是驗證了這個 75 步是最佳解了。理論上來說，用 Mathematica 應該是也能寫得出來啦，不過還沒怎麼研究在 Mathematica 裡寫 priority queue 比較有效率所以就暫且放著了；說不定之後想到了會回來寫？(笑)