作者LPH66 ( )
看板Math
標題Re: [其他] 2021 新年數學題目
時間Fri Jan 8 16:14:07 2021
※ 引述《TimcApple (肥鵝)》之銘言：
: https://reurl.cc/N6zekQ
: 一共 10 題和 2021 有關的題目
: 祝各位新年快樂ow o
: 推 MisatoMitumi: Q3先給一個157步的參考路徑（不過確定有更佳解） 01/03 19:57
: → MisatoMitumi: （只列出取階乘的數字）2021->26->2116->973->268-> 01/03 19:58
: → MisatoMitumi: 124->1726->53226->3059->36191->35910->32820->110 01/03 19:59
: → AnnaOuO : 看不懂樓上數字變化.. 01/05 11:37
: 推 LPH66 : 以第一個箭頭為例, 2021 取階乘後開十二次平方根得 01/05 12:11
: → LPH66 : 約 26.13, 在這裡向下取整得 26; 下一步是 26! 01/05 12:13
: → ejialan : M大列的數列經過取階乘->取n次平方根->高斯記號 01/05 12:13
: → LPH66 : 其開三次平方根後得約 2116.91 向下取整得 2116 01/05 12:13
: → LPH66 : 依此類推 01/05 12:13
: → ejialan : n依序為12 3 11 10 8 6 10 16 11 15 15 16 01/05 12:18
: 推 LPH66 : OK, 第三題寫程式跑個開頭就看到為何確定有更佳解了 01/05 18:39
: → LPH66 : 26 可以不用上面說的 14 步獲得: 2021 不先做階乘 01/05 18:40
: → LPH66 : 而是開兩次根號再取整得 6, 階乘 720 開方取整得 26 01/05 18:41
: → LPH66 : 這樣只需要 6 步就能到 26 了 01/05 18:41
: 推 LPH66 : 然後這支程式找到以下這個 88 步解: 接 26, 走 01/05 18:52
: → LPH66 : 26-[4]>46-[4]>4062-[13]>37-[4]>496-[8]>24427 01/05 18:53
: → LPH66 : -[14]>784453-[21]>110; -[n]> 表開 n 次平方再取整 01/05 18:54
: → LPH66 : 不過這依然無法保證是最佳解, 因為程式計算關係 01/05 18:56
: → LPH66 : 中間取階乘只在 < 10^15 的整數上用, 無法確定當 01/05 18:57
: → LPH66 : 允許更大數時有沒有可能有更少步數的解 01/05 18:58
: → LPH66 : 啊更正一下, -[n]> 表示階乘後開 n 次平方再取整 01/05 19:20
: → LPH66 : 漏提了階乘步驟 XD 01/05 19:20
: 推 MisatoMitumi: 88比預期的要低很多啊... 我是從110逆推的，所以202 01/05 20:05
: → MisatoMitumi: 1->26可以更短是結束前才發現的XD 01/05 20:05
: → MisatoMitumi: LPH大有辦法處理到(10^15)!附近嗎？我自己的方法是 01/05 20:08
: → MisatoMitumi: 類似試著建立(10^5)內的有向圖，然後跑Dijkstra最短 01/05 20:08
: → MisatoMitumi: 路徑，不過找到路徑就停了 01/05 20:08
: 推 MisatoMitumi: 原來如此，計算Log(n!)的複雜度不是O(n)，看來我太 01/05 20:16
: → MisatoMitumi: 沒效率了 01/05 20:16
: 推 LPH66 : 對, 我原本以為大概也要 O(n), 不過想說這總和看著 01/05 21:50
: → LPH66 : 像是個積分, 那似乎能用∫log(n) ~ n log(n) 估計 01/05 21:50
: → LPH66 : 然後一查才發現拉馬努金有給出過很好的近似公式 01/05 21:51
: → LPH66 : 大致形狀是 n log(n) - n + log(n 的三次式) 01/05 21:53
: → LPH66 : 所以能算的 n 就變成 double 能表達的整數 ~9e15 了 01/05 21:54
: → LPH66 : 不算到極限是因為不太能掌握精確度... 01/05 21:55
: → LPH66 : 然後我沒有明確建圖, 而是把 階乘→開方 n 次→取整 01/05 21:58
: → LPH66 : 這樣的 n+2 步對不同的 n 值散出去 01/05 21:59
: → LPH66 : 核心還是 Dijkstra 但就是有踩到哪些點再記那些點 01/05 22:00
: 推 MisatoMitumi: 原來如此，公式有趣！不過wiki給的誤差大約是1/1400 01/05 23:08
: → MisatoMitumi: n^3，算是有點怕... 不知道現在找到最佳解頂多88步 01/05 23:08
: → MisatoMitumi: 的話，爆出最佳解有沒有機會... 01/05 23:08
糟糕, 這兩天回頭看這支程式發現它有一個溢位 bug
改正之後它找到一個更少步的解了:
2021 開兩次根號取整得 6, 然後
6
－[2]→ 5
→ 120
─[6]→ 1278
─[8]→ 22280837523899
─[47]→ 110
總計 75 步
問題在於最後這一步: 22280837523899! 開四十七次平方根
我用的 log(n!) 的近似公式是在這裡寫的這一條：
https://math.stackexchange.com/a/138326
(有其他回答有提到那個 log 裡的 n 的三次式應該要加 1/30, 不過這裡差別不大)
因為有點掌握不住這個式子的誤差所以一直不敢完全信任 (雖然是拉馬努金提出的但..)
但十四位數的階乘又不可能真的直接去求...
是有試用過 Mathematica 的 NSum 驗證 (用 Euler-Maclaurin formula 由積分推算的)
給出的 log(22280837523899!) 值確實是在 2^47*110 和 2^47*111 之間:
2^47*110 ~ 6.615338*10^14, 2^47*111 ~ 6.6280745*10^14
而 NSum 給出的 log(22280837523899!) (和拉馬努金的公式) 都是 ~ 6.6251509*10^14
所以...這個四十七次根號的計算應該是對的吧?
====
至於那個溢位 bug
是因為我的 log(n!) 函數的參數寫錯了, 會把 64-bit 整數傳成 32-bit 整數進來
所以像上面這種數就會被當成負數傳進去了 orz
====
說起來,
: → MisatoMitumi: 不知道現在找到最佳解頂多88步 01/05 23:08
: → MisatoMitumi: 的話，爆出最佳解有沒有機會... 01/05 23:08
關於這一點我有想過了
因為太大的數字要使用平方根降下來也不容易
因此在已知有個總步數上限的狀況下確實能夠給出一個能使用階乘的上限值 M
當開這麼多次平方根還降不回到這個上限值範圍就表示不會是最佳解
就拿現在已知的 75 步來算
如果用 log(M!) ~ M*log(M) 來估的話, 降回到 log(M) 要在 75 次除以 2 內完成
(這還是個很粗略的條件)
但這差不多就是除以了 M, 也就是說 M 最大可以到 2^75 都是可能範圍
這個範圍好像只比常見的程式整數 / 浮點數的精確度多了一些些而已...
--
1985/01/12 三嶋鳴海 1989/02/22 優希堂悟 1990/02/22 冬川こころ 1993/07/05 小町
つぐみ 歡迎來到 1994/05/21 高江ミュウ 1997/03/24 守野いづみ 1997/03/24 伊野瀬
チサト 1998/06/18 守野くるみ 打越鋼太郎的 1999/10/19 楠田ゆに 2000/02/15 樋口遙
2002/12/17 八神ココ 2011/01/11 HAL18於朱倉岳墜機 ∞與∫的世界 2011/04/02 茜崎空
啟動 2012/05/21 第貮日蝕計畫預定 2017/05/01~07 LeMU崩壞 2019/04/01~07 某大學合宿
--
※ 發信站: 批踢踢實業坊(ptt.cc), 來自: 106.1.234.196 (臺灣)
※ 文章網址: https://www.ptt.cc/bbs/Math/M.1610093649.A.080.html
※ 編輯: LPH66 (106.1.234.196 臺灣), 01/08/2021 16:17:00
推 sunev : 啊，mathematica有內建LogGamma，應該是不用NSum 01/08 20:51
推 MisatoMitumi: 看來是最佳解了，我難過 01/09 05:44
→ MisatoMitumi: 用Magma驗證過了，有內建的LogGamma也有內建的自訂 01/09 05:44
→ MisatoMitumi: 有興趣的可以自行組裝XD，執行很快der 01/09 05:47
推 MisatoMitumi: 不過程式寫得很醜~ 01/09 05:49
→ TimcApple : 推推 01/09 10:55