BCC-16 (in Chinese) 計算機概論十六講 C -- more involved iterations

C 教材：再多點兒技巧的迭代範例

前一節，我們反覆地舉例，讓讀者觀察基本的迭代範例。透過模擬與練習，可能讀者已經掌握了最基本的語法。前一節的基本迭代有個特色，就是每一次迭代只需要計算一個單項。例如計算

$\sum_{n=1}^{1000} \frac1{n^2}$

則每次迭代過程只要計算 1.0/(n*n)，而 n 的數值在 for 迴圈中決定。

這一節，我們稍微增加一點曲折，再舉出幾個範例，考慮那些迭代過程中不只有單項計算的例子。

【例一】計算

$\begin{displaymath} \sum_{n=0}^{52} (\frac12)^n \end{displaymath}$

並以 %f 格式輸出。讀者一看就知道這是以 1/2 為公比的等比級數，因為如果有無窮多項的極限是 2，所以這個計算的結果也應該差不多是 2。這個問題和前一節的「單純」問題的主要差異，是每一次迭代過程中要加進來的項是 (0.5)ⁿ，好像要計算一個次方。其實不然，只要另外宣告一個變數來儲存被加的項，利用迭代過程將這個「被加項」逐次更新就行了。以下是計算的程式。

#include <stdio.h>
main() {
    int n;
    double sum, x;
    sum = 0.0;
    x = 1.0;
    for (n=0; n<=52; ++n) {
	sum = sum + x;
	x = x/2;
    }
    printf("%f\n", sum);
}

輸出應該是，毫無意外地：2.000000。

以上程式中，我們用 sum 來儲存和，所以它的初始值設定為 0.0；用 x 來儲存「被加項」，所以它的初始值是 1.0 (當 n=0 的時候)。而後，每次迭代的時候，利用 x = x/2 這句話將「被加項」更新，使得 n=1 的時候 x 是 1/2，而 n=2 的時候 x 是 (1/2)²，依此類推。

以上程式大概是最中規中矩的寫法了。稍微取巧一點，我們知道題目的第一項是 1，所以可以不做 n=0 那一項，而從 n=1 那一項做起。程式就可以修改成以下模樣。我把有變動的地方用藍色字型表現。

#include <stdio.h>
main() {
    int n;
    double sum, x;
    sum = 1.0;
    x = 1.0;
    for (n=1; n<=52; ++n) {
	x = x/2;
	sum = sum + x;
    }
    printf("%f\n", sum);
}

留意 x=x/2; 那句話搬到 sum=sum+x; 的前面了。這樣才能使得第 n 次迭代時，被加項是 (1/2)ⁿ。

再用上 C 語言的一點技巧，前述程式可以再改寫得比較短一點。

#include <stdio.h>
main() {
    int n;
    double sum, x;
    x = sum = 1.0;
    for (n=1; n<=52; ++n) {
	x = x/2;
	sum = sum + x;
    }
    printf("%f\n", sum);
}

這個程式還可以利用更複雜的 C 語言技巧讓它更短。但是對於初學者而言，沒必要學取太多技巧，先掌握基本思考方法比較重要。

【例二】計算以下連加，並以 %f 格式輸出。

$\begin{displaymath} \sum_{n=0}^{30} \frac1{n!} \end{displaymath}$

這個題目的「被加項」是 1/n!，看起來需要特別的函式來計算，其實也是可以利用另一個變數來逐次計算它。程式如下。

#include <stdio.h>
main() {
    int n;
    double sum, x;
    x = sum = 1.0;
    for (n=1; n<=30; ++n) {
        x = x/n;
        sum = sum + x;
    }
    printf("%f\n", sum);
}

學過微積分的讀者，知道上述題目的極限是常數 e，也就是 2.71828...。而以上程式的輸出是 2.718282。

因為我們知道題目中的第一項是 1 (0! 定義為 1)，所以從 n=1 開始計算。而 n=1 時，「被加項」x 是 1/1 也就是 1；到了 n=2 時，x 是 x/2 也就是 1/2；到了 n=3 時，x 是 x/3 也就是 (1/2)/3 亦即 1/6 或 1/3!。依此類推，x 就是 1/n!。

【例三】令 x 為 0.69314718 (差不多是 ln(2): 2 的自然對數)，計算以下連加並以 %f 格式輸出。

$\begin{displaymath} \sum_{n=0}^{30} \frac1{n!}x^n \end{displaymath}$

這一題也不過是把前面兩個例子合併在一起。我們先中規中矩地寫一個程式，讓「被加項」等於 1/n! 和 (0.69314718)ⁿ 相乘。我們用 coef 來儲存 1/n!，用 x 來儲存 (0.69314718)ⁿ。

#include <stdio.h>
main() {
    int n;
    double sum, x, coef;
    sum = x = coef = 1.0;
    for (n=1; n<=30; ++n) {
        coef = coef/n;
        x = x*0.69314718;
        sum = sum + coef*x;
    }
    printf("%f\n", sum);
}

如果讀者知道微積分，那麼知道這個計算的結果應該差不多是 2。而以上程式的輸出是 2.000000。

以上程式還有兩點改進的空間。第一，我們在 [原始碼的可讀性、符號常數] 這一節說過，程式中盡量不要有神秘的常數。這裡，0.69314718 是個神秘的常數。我們可以用一個符號常數來代表它。以後如果要修改這個程式，使它計算另一個 x 值的同樣連加，就只要修改那個符號常數。第二，仔細想想，「被加數」其實沒有必要分成兩項的乘積，可以刪除 coef 變數而只用 x 變數就行了。

#include <stdio.h>
#define VAR 0.69314718  /* VAR is the input variable */
main() {
    int n;
    double sum, x;
    sum = x = 1.0;
    for (n=1; n<=30; ++n) {
        x = x*VAR/n;
        sum = sum + x;
    }
    printf("%f\n", sum);
}

【例四】令 x 為 1.5707963 (大約是 Pi/2: 圓周率的一半)，計算以下連加並以 %f 格式輸出。

$\begin{displaymath} \sum_{n=0}^{30} (-1)^n \frac1{(2n+1)!}x^{2n+1} = x - \frac1{3!}x^3 + \frac1{5!}x^5 - \frac1{7!}x^7 + \cdots \end{displaymath}$

知道微積分的讀者，應該看得出來這是 sin(x) 的 Taylor 級數展開，因此答案應該差不多是 1。

這個問題的難處在於需要正負項交替出現。我們目前還沒有這種「若 n 為奇數則負、否則正」的技巧 (需要 if 語法)，但是其實也不需要：我們可以運用「正負兩項一起做」的技巧。也就是說，每次迭代就計算兩項，例如 n=0 和 n=1 一起做， n=2 和 n=3 一起做。但是，要在迭代過程中計算 n=0 那一項會有困難，我們不在這裡做錯誤的示範，請讀者自己比照以下程式想想看，如果有 n=0 那一項，有何麻煩？因此我們決定從 n=1 開始做，而把 n=0 的數值當作 sum 和 x 的初始直。因為「一次做兩項」，所以 n=1 和 n=2 一起做、... n=29 和 n=30 一起做。

#include <stdio.h>
#define VAR 1.5707963  /* VAR is the input variable */
main() {
    int n;
    double sum, x;
    x = VAR;
    sum = x;
    for (n=1; n<=30; n=n+2) {
        x = x*VAR/(2*n);
        x = x*VAR/(2*n+1);
        sum = sum - x;
        x = x*VAR/(2*n+2);
        x = x*VAR/(2*n+3);
        sum = sum + x;
    }
    printf("%f\n", sum);
}

輸出的答案是 1.000000。

以上程式雖然正確，但是未免太過於平鋪直敘，以至於明顯地缺乏效率。譬如說 2*n 居然計算了四遍。仔細看看「被加項」的變化，是

-x³/3!, +x⁵/5!, -x⁷/7!, +x⁹/9!, ... -x⁵⁹/59!, +x⁶¹/61!

何不直接讓 for 迴圈產生 n = 3, 7, ... 59 這些參數？而「被加項」第一項 (負項) 就成了

x * VAR² / ((n-1)*n)

第二項 (正項) 就接著變成

x * VAR² / ((n+1)*(n+2))

以下是修改的程式。

#include <stdio.h>
#define VAR 1.5707963  /* VAR is the input variable */
main() {
    int n;
    double sum, x;
    sum = x = VAR;
    for (n=3; n<60; n=n+4) {
        x = x*VAR*VAR/(n*n-n);
        sum = sum - x;
        x = x*VAR*VAR/(n*n+3*n+2);
        sum = sum + x;
    }
    printf("%f\n", sum);
}

輸出的答案當然還是 1.000000。

上述兩種程式寫作的範例，一種比較平鋪直敘，稍微沒有效率，但是容易寫也容易看懂。另一種比較多技巧和事前的代數推演，因此比較有效率，但是比較不容易看懂 (其實也比較容易寫錯)。因此，在「效率」和「方便」之間，程式設計者要有所取捨。這種取捨並沒有公式或規則可循，這表現了個人的經驗和風格。

習題

令 x 為 3.1415926535 (大約是圓周率)，計算以下連加並以 %f 格式輸出。
$\begin{displaymath} \sum_{n=0}^{30} (-1)^n \frac1{(2n)!}x^{2n} = 1 - \frac1{2!}x^2 + \frac1{4!}x^4 - \frac1{6!}x^6 + \cdots \end{displaymath}$
令
$\begin{displaymath} a_n(x) = \frac{1\cdot 3\cdot 5\cdots (2n-1)}{2\cdot 4\cdot 6\cdots(2n)}\, \frac{x^{2n+1}}{2n+1} \end{displaymath}$
定義
$\begin{displaymath} f(x) = x + \sum_{n=1}^N a_n(x) \end{displaymath}$
其中 N 是使得 $\vert a_N(x)\vert\leq 10^{-12}$ 的最小正整數。用 C 語言寫程式計算
$\begin{displaymath} f(\frac12),\quad f(\frac1{\sqrt3}),\quad f(\frac1{\sqrt2}) \end{displaymath}$
用 double 資料型態，以 %.8f 格式輸出答案。您能否看出來，若 N 趨向無限大，f(x) 趨近於哪個函數？


	單維彰 (2001/01/06) ---