Matlab做project euler 14怎樣提速？

01-08

@吳運生@李海濤
@yellow 從數學上分析減小計算量，這是從最根本上解決問題。
不過即使不做數學分析，就題主本身的程序也有可改進的地方。下面的分析完全不用到對問題本身的數學方面的分析，即不對整個過程做任何篩選。
先說下李海濤說的python方法，其測試代碼如下（只修改了一個地方，把abc/2改為abc&>&>1了，這樣大概能快1/3）：

N = 1000000 - 1 A = {1: 1} def collatz(abc): if abc in A: return A[abc] elif abc % 2: A[abc] = collatz(3 * abc + 1) + 1 return A[abc] else: A[abc] = collatz(abc &>&> 1) + 1 return A[abc] imax = 1 nmax = 1 for x in range(2, N + 1): t = collatz(x) if t &> imax: imax = t nmax = x print(imax, nmax)
這段代碼我機器上跑 1.5 秒，如果換成 python 2.7 能跑 1.38 秒
然後是我說的只緩存比較小的數的寫法：
def collatz(n): l, lmax = {1: 1}, 0 for k in range(2, n + 1): a, lk = k, 0 while a &>= k: lk += 1 a = 3*a + 1 if a % 2 else (a &>&> 1) lk += l[a] l[k] = lk if lk &> lmax: f, lmax = k, lk return f, lmax r = collatz(1000000)
這段代碼跑 1.3 秒，在2.7上是 0.85 秒，如果 l 用 list 的話則分別是 1.2 秒和 0.76 秒
當然如果用 Pypy 之類的可以加速。而用 cython 的話可以快到和 C 差不多，不過那代碼也快和 C 差不多了
然後再說 MATLAB，題主的 MATLAB 代碼和李海濤的基本一樣的，所以就以題主的代碼為例，先運行一下（版本是windows 10 下的 R2016a），運行時間如下：
&>&> tic; pe14; toc 837799 時間已過 6.401452 秒。
然後用 profile 分析性能：

（這裡由於在 profile 模式下運行所以整個程序的運行時間也變長了，達到了 21 秒）
可以看到最慢的是
if mod(n1,2) == 0
執行該行佔據用了整個運行時間中相當多的比例。根據經驗我們知道 MATLAB 對大量重複對標量浮點數做 mod 操作應該是很慢的，所以考慮改成:
if floor(n1/2)*2 == n1
再試一下，發現提速還挺明顯的:
&>&> tic; pe14; toc 837799 時間已過 1.849706 秒。
如果緩存的話，同樣按照上邊的寫法：
function [f,lmax] = pe14(n) [l(n), l(1), lmax] = deal(0, 1, 0); for k = 2:n a = k; lk = 0; while a &>= k lk = lk + 1; if a == floor(a/2)*2, a = a*.5; else a = a*3 + 1; end end lk = lk + l(a); l(k) = lk; if lk &> lmax, lmax = lk; f = k; end end
結果：
&>&> tic; [f,lmax] = pe14(1e6); toc 時間已過 0.084118 秒。

不會Matlab，本人用的是C++實現：
#include& const int n=1000000; int f[n+1]; int main(){ int a=1,l=0; for(int i=2;i&<=n;i++){ long long x=i; int t=0; while(x&>=i){ if(x1)x+=x+1&>&>1,t+=2; else x&>&>=1,t++; } f[i]=t+=f[x]; if(t&>l)a=i,l=t; } printf("%d ",a); }
幾個優化：
1、用位運算代替部分算術運算。這個優化僅起到常數優化的效果（有時編譯優化會幫你做這個事，所以甚至不一定有用）。

2、當x是奇數時，3x+1一定是偶數，也就是下一次一定變成(3x+1)/2。所以當x為奇數時，可以直接模擬兩次操作得到x+(x+1)/2。同樣這也只是常數優化而已。
3、記憶化。注意到我們依次計算1,2,3,...,n的答案，在計算n的最少步數時，如果某一步x&

yellow最後用了Python，哈哈，戳到MATLAB的痛處了。
這題我之前做過，試了很多方法，最後發現還是直接用暴力法最快，為什麼呢？因為MATLAB沒有效率高的字典型數據結構。
MATLAB：
function [result,len_max] = euler014 len_max = 0; result = 0; for ii=1:1e6-1 n = ii; len = 1; while n~=1 if mod(n,2)==0 n = n/2; else n = n*3+1; end len = len+1; end if len&>len_max len_max = len; result = ii; end end disp(result); end
tic;euler014;toc; 837799 時間已過 8.519639 秒。
Python:
from time import clock
tbegin = clock() N = 1000000 - 1 A = {1: 1} def collatz(abc): if abc in A: return A[abc] elif abc % 2: A[abc] = collatz(3 * abc + 1) + 1 return A[abc] else: A[abc] = collatz(abc / 2) + 1 return A[abc]
imax = 1 nmax = 1 for x in range(2, N + 1): t = collatz(x) if t &> imax: imax = t nmax = x print(imax, nmax) print(clock() - tbegin)
結果：
(525, 837799) 2.0769729404
電腦配置：

總結：
MATLAB用了暴力法，因為沒有高效的字典型數據結構。Python用字典儲存了中間結果，速度快了很多。

瀉藥，粗略看了一下題主的代碼，應該是不帶任何分析的暴力循環。
其實這個問題稍微轉個彎就能將篩選次數壓縮到一半，方案就是直接從 $[500000, 1000000]$ 內循環查找就可以了。
就這點簡要分析一下，權當拋磚引玉了。
根據題目規則，無非就是當給定初始值 $a_0$ ，分兩種情形處理：

若 $a_0$ 為偶數，則 $a_0 ightarrow frac{a_0}{2}$ （該值可能為奇數，也可能為偶數）

若 $a_0$ 為奇數，則 $a_0 ightarrow 3a_0+1$ （該值一定是偶數）

所以最後產生的整個序列：
$a_0 ightarrow a_1 ightarrowcdots ightarrow 2 ightarrow 1$
無非就是在奇數與偶數間來換轉換，其中奇數經過一次變換為偶數，偶數經過若干次變換為奇數。

因此該序列中任意一個奇數（非序首）前面一定是偶數，而偶數前面則不確定。

所以，當你給定的初值 $a_0in[1, 500000]$ 時，顯然可以再對該序列擴展一個，
$2a_0 ightarrow a_0 ightarrow a_1 ightarrowcdots ightarrow 2 ightarrow 1$
--------------------------------------------&>&>&> 進一步篩選 &<&<&<--------------------------------------------
基於以上分析，我們已經粗略地鎖定了序首 $a_0$ 一定在 $[500000, 1000000]$ 內出現。那麼，我們不妨假設最長的序列 ${a_i}$ 其對應的 $a_0$ 為偶數，同時假定其具有如下形式：
$a_0 = 6n+4, quad n in N$

則 $a_0$ 可由一個奇數 $b_0=2n+1$ 變換得到： $a_0=3b_0+1$
不難發現該情形下， $b_0in[1, 333333]subseteq [1, 500000]$ ，於是可以構造更長序列：
$b_0 ightarrow a_0 ightarrow a_1 ightarrowcdots ightarrow 2 ightarrow 1$
矛盾於之前關於 ${a_i}$ 為最長序列的假設，因而但凡具有 $a_0 = 6n+4$ 形式的初值一定不是最優解.
於是在當你在 $[500000, 1000000]$ 範圍內搜索時，可以事先判斷 $mod(a_0, 6)$ 是否等於4，若是則直接進入下一次搜索。這種篩選方案還能進一步排除掉其中的 $83334$ 種情況。
那麼我們嘗試在該思路下測試代碼效率：
tic [seq_len, max_start] = deal(0, 1); for i = 5e5:1e6 a0 = i; if mod(a0, 6) == 4 continue end
n = 1; while a0 &> 1 if mod(a0, 2) == 0 a1 = a0 / 2; else a1 = 3*a0 + 1; end n = n + 1; a0 = a1; end
if n &> seq_len [seq_len, max_start] = deal(n, i); end end toc
發現運行時間基本穩定在8 sec左右，也並不是很快。
之後，只能換個思路了，嘗試利用緩存機制來優化一下：

在進行第n+1次迭代搜索時，藉助前面n次搜索的緩存結果

遞歸思路：

根據@觀察者的遞歸思路，再結合 @Falccm 提出的n/2改為右移操作n&>&>1，在其原有代碼上進行相應改動可以跑到2.43 sec左右：

# coding: utf-8
import time collatz_cache = {1: 1} # 初始化計算緩存字典 def collatz(n): """遞歸求解指定數字n對應的序列長度""" if n in collatz_cache: return collatz_cache[n] t = collatz(3*n+1) + 1 if n 1 else collatz(n&>&>1) + 1 collatz_cache[n] = t return t def max_collatz_len(n=1000000): """遍歷指定區間內最長序列對應的數字""" max_len, number = 1, 1 for i in range(n//2, n+1): #if i % 6 == 4: # continue i_len = collatz(i) if i_len &> max_len: max_len, number = i_len, i return number, max_len
if __name__ == "__main__": start_time = time.clock() result, max_len = max_collatz_len() print("number:".rjust(20, " "), result) print("max length:".rjust(20, " "), max_len) print("time elapsed:".rjust(20, " "), "%.2f sec" % (time.clock()-start_time))
原始代碼見：ProjectEuler/pr014.py at master · KNCheung/ProjectEuler · GitHub

直接緩存機制：

@Falccm提出的「小數緩存」機制（在計算n+1對應的序列長度時，當序列首次迭代到比n+1小的數時，則利用之前的緩存結果直接計算），在我電腦上跑基本穩定在2.84 sec附近（電腦配置確實老古董了-_-||），不比富帥們直接衝進&<1 sec區。
改用list緩存基本維持在2.83 sec附近，
# coding: utf-8
import time collatz_cache = [0] * 1000001 collatz_cache[1] = 1 def collatz(n): start, depth = n, 0 while start &>= n: depth += 1 start = 3*start + 1 if start 1 else start &>&> 1 depth += collatz_cache[start] collatz_cache[n] = depth return depth def max_collatz_len(n=1000000): result, max_len = 1, 1 for i in range(2, n+1): depth = collatz(i) if depth &> max_len: result, max_len = i, depth return result, max_len
if __name__ == "__main__": start_time = time.clock() result, max_len = max_collatz_len() print("number:".rjust(20, " "), result) print("max length:".rjust(20, " "), max_len) print("time elapsed:".rjust(20, " "), "%.2f sec" % (time.clock()-start_time))
這樣一種緩存思路需要連續緩存從 $1 ightarrow N$ 的全部計算結果，於是我之前提出的篩選規則 $[frac{N}{2}, N]$ 就可能很難直接起作用了，簡單的解決方案是在搜索迭代的同時緩存該搜索序列，例如在查找序列
$a_0 ightarrow a_1 ightarrowcdots ightarrow a_k ightarrow cdots ightarrow 2 ightarrow 1$ 時，當首次出現 $a_k$ 在我們的緩存結果中時立即終止，同時錄入 $a_0, a_1, dots, a_{k-1}$ 的計算結果。
當然，這意味著可能會把一些超出 $[0, N]$ 範圍的計算結算緩存進去，緩存的最大數量也不可知，不過這一點可以通過判斷後過濾掉。以下是一個簡單實現，速度基本在3.60 sec：
# coding: utf-8
import time collatz_cache = {1: 1} # 初始化緩存字典 def collatz(n): """通過緩存機制計算制定數字n對應的序列長度""" start, depth = n, 0 search_path = [start] while start not in collatz_cache: depth += 1 start = 3*start + 1 if start 1 else start &>&> 1 search_path.append(start) depth += collatz_cache[start] for i, num in enumerate(search_path): collatz_cache[num] = depth - i return depth def max_collatz_len(n=1000000): """遍歷指定區間內最長序列對應的數字""" max_len, number = 1, 1 for i in range(n//2, n+1): #if i % 6 == 4: # continue i_len = collatz(i) if i_len &> max_len: max_len, number = i_len, i return number, max_len
if __name__ == "__main__": start_time = time.clock() result, max_len = max_collatz_len() print("number:".rjust(20, " "), result) print("max length:".rjust(20, " "), max_len) print("time elapsed:".rjust(20, " "), "%.2f sec" % (time.clock()-start_time))
總結來看，執行效率大概是"遞歸" &> "連續緩存" &> "離散緩存"，需要特別指出的是，在"連續緩存"中我們並未加入篩選規則（沒法加入）。

既然你做pj,難道不知道用一種方法完成之後就可以去看討論帖
裡面有人們提供的各種解法

有很多人回答了問題，不過似乎都沒提出一些很激進的優化吧。。我來提一個:角谷猜想驗證是可以記憶化最後幾個二進位位的操作的，這樣可以大大加速。（不過會有點難寫

Mathematica湊個熱鬧
cf1 = Compile[{{n0, _Integer}}, Block[{n = n0}, Catch@Do[If[n == 1, Throw@i]; n = If[BitAnd[n, 1] == 0, Floor[n/2], 3 n + 1] , {i, n0}] ], RuntimeAttributes -&> Listable, CompilationTarget -&> "C", RuntimeOptions -&> "Speed" ];
Ordering[cf1@Range[10^6], -1] // AbsoluteTiming cf2 = Compile[{{n, _Integer}}, Module[{A, j, t}, A = ConstantArray[0, n]; A[[1]] = 1; Do[ j = i; t = 0; While[j &>= i, t++; If[BitAnd[j, 1] == 0, j = Floor[j/2], j = 3 j + 1]; ]; A[[i]] = A[[j]] + t; , {i, 2, n}]; Ordering[A, -1] ], CompilationTarget -&> "C", RuntimeOptions -&> "Speed" ];
cf2[10^6] // AbsoluteTiming

8s已經很快了。。有些題不會優化，智商不夠，時間來湊，跑個8h都是常事

沒寫過MATLAB，這個題可以把算好的結果存下來。Python實現如下：
def collatz(): memo = dict()
def iter(n, s): if n in memo: step = s + memo[n] return step else: if n == 1: return s else: result = 3*n+1 if n % 2 else n // 2 step = s + iter(result, s) memo[n] = step return step max_n, max_step = 0, 0 for i in xrange(1, 1000001): step = iter(i, 1) if step &> max_step: max_n, max_step = i, step return max_n, max_step
print(collatz())
i5 2.6GHz跑的

Python2 1.5～1.6秒之間
PyPy 1.0～1.1秒之間

推薦閱讀：