Matlab 教材：for-if-break 技術介紹 3x+1 問題

令 x0 是一個正整數 (> 0 的整數)，定義一個數列

\begin{displaymath} x_{n+1} = \left\{\begin{array}{ll}x_n/2&\mbox{if $x_n$ is even}\\ 3x_n+1& \mbox{else}\end{array} \right. \end{displaymath}

此數列必然或進入 1-4-2 循環。 這就是俗稱的 3x+1 問題，學名是 Collatz 問題。 而按照上述規則產生的數列，稱為 Collatz 數列。

所謂 1-4-2 循環，就是說如果某個 xk = 1，則按照規則 xk+1 = 3*xk+1 = 4，然後按照規則 xk+2 = xk+1/2 = 2，然後按照規則 xk+3 = xk+2/2 = 1，然後按照規則 xk+4 就又是 4，然後又是 2，然後又是 1，如此循環不已。

因為只要有某個 xk = 1，就會進入 1-4-2 循環， 所以我們只要遇到第一個 1 就該停止計算，否則再算下去也沒什麼意義。 如果我們把規則訂成

在第一次遇到 xk = 1 時停止產生 Collatz 數列
那麼 3x+1 問題又可以改編成：
不論 x0 是什麼正整數，Collatz 數列一定會停下來

我相信絕大部分的讀者都感到疑惑：真的是這樣嗎？ 直覺上實在很難想像為什麼上述產生 Collatz 數列的程序一定會停下來？ 我們利用 Matlab 來實驗看看。

為了害怕萬一它永遠不停下來 (無窮迴圈) 就麻煩了， 所以我們固定讓它只迭代 1000 次，也就是最多只產生到 x1000 而已； 如果在 x1000 之前就遇到 xk=1 就提早停下來。 讓我們先試試看，若 x0 = 3 會怎樣？ 紅色數目字代表指令的編號，不是指令的一部分。 

 1 x = 3;
 2 for n=1:1000
 3     if (rem(x,2) == 1)
 4         x = 3*x+1
 5     else
 6         x = x/2
 7     end
 8     pause(1);
 9     if (x == 1)
10         break;
11     end
12 end
以上程式可以直接寫在 Matlab 介面內，但是最好寫成一個腳本程式。 在 1 我們設定 x0，然後在 2 -- 12 之間我們依序產生 x1, x2, ... 最多到 x1000。 而 3 -- 7 就是按照規則決定新的 xn，其中 3 中的 rem(x,2) 是 x 被 2 除的餘數 (remainder)， 也就是說「若 x 是奇數」的意思。我們故意在 4 和 6 的指令結束後不加分號， 讓 Matlab 把計算結果輸出在螢幕上，以便觀察。 在 8 我們要 Matlab 暫停 (pause) 1 秒鐘， 好讓我們有時間看看計算的結果是什麼。 接著 9 -- 11 就是標準的 if-break 技術：如果有某個 xn=1，那就不必再繼續下去了。 break 的意思就是離開包含這個 if 指令的 for 迴圈。

為了讓讀者方便剪貼上述程式到 Matlab 工作環境之中，以下我要取消那些紅色的編號。 順便，我要再簡化 3： 因為 rem(x,2) 的結果不是 0 就是 1，而 Matlab 本來就把 0 視為 FALSE、 非 0 視為 TRUE。因此如果 x 是奇數，則 rem(x,2) 必然為 TRUE； 如果 x 是偶數，則 rem(x,2) 必然為 FALSE。所以

if(rem(x,2) == 1)
這句話是「脫了褲子放屁」：多此一舉，只要寫
if(rem(x,2))
就足夠表達「若 x 是奇數」這個意思了。簡化後的程式如下： 
x = 3;
for n=1:1000
    if (rem(x,2))
        x = 3*x+1
    else
        x = x/2
    end
    pause(1);
    if (x == 1)
        break;
    end
end

現在，請您剪貼程式到 Matlab 環境中，試試看， 由 x0=3 所產生的 Collatz 數列是

3, 10, 5, 16, 8, 4, 2, 1
讀者可以自行修改 x0 的值 (也就是程式的第一列)， 試試看您能不能發現一個會使得 Collatz 數列在 1000 項以內還沒出現 1 的 x0？ (不要來問我答案)

習題

  1. 定義 Collatz 數列的「長度」為使得 xk=1 的最小 k。 令 x0 取值 3, 4, ..., 15， 哪一個 x0 導致的 Collatz 數列最長？有多長？ (您可能要適度修改這一節展示的程式，讓它可以顯示 Collatz 數列的「長度」。)
單維彰 (2003/06/15) --- 04/04/28 (單) [Prev] [Next] [Up]