39405 我家$\pi$的故事—兼談電子計算機在數學上的威力

終極密碼

遊戲規則:本遊戲為猜密碼的遊戲。密碼為0到100之間的其中1個整數,電腦會提示密碼的所在範圍,玩家必須在6次之內猜到密碼才能過關。
★ 終極密碼為0到100之間 ★
您共有六次機會

一、前奏 --- 圓周率的發展歷史

圓周率, 一般用希臘字母 $\pi$ 來表示, 它是計算圓周長、 圓面積、 球體積等幾何量的關鍵常數。 $\pi$ 的定義是圓的周長與直徑的比值, 大約等於 3.14; 它也等於圓的面積與半徑平方的比值。 $\pi$ 的計算, 歷經四千年, 最新的結果是計算到六十兆位小數。 下表是一個簡略的介紹。

日期 計算者 $\pi$ 的值
前20世紀 巴比倫人 3.125
前20世紀 印度人 3.160493$\cdots$
前12世紀 中國人 3
前6世紀 聖經列王記上7章23節 3
前3世紀 阿基米德 3.1418
480年 祖沖之 $3.1415926\lt \pi\lt 3.1415927$
1400年 Madhava 3.14159265359
1596年 魯道夫$\cdot$范$\cdot$科伊倫 20位小數
1699年 Abraham Sharp 71位小數
1706年 John Machin 100位小數
1706年 William Jones 引入希臘字母 $\pi$
1761年 Johann Heinrich Lambert 證明 $\pi$ 是無理數
1853年 William Shanks 527位小數
1882年 Lindemann 證明 $\pi$ 是超越數
1949年 J. W. Wrench 和 L. R. Smith 2037位小數 (首次使用計算機)
1961年 IBM 7090電晶體計算機 20,000位小數
1981年 金田康正 2,000,000位小數
1994年 楚諾維斯基兄弟 4,044,000,000位小數
2010年 近藤茂 5,000,000,000,000位小數
2011年 IBM 藍色基因/P超級計算機 60,000,000,000,000位小數

二、幕起 --- 我的外孫叫小 $\pi$

2013年我們家熱烈討論的議題之一是, 住在美國德州的大女兒年底要生的外孫的名字。 其中討論最熱烈、最難達到共識的是他的中文名字。幾經轉折, 最後以「信恆」定案, 取其信心與恆心。 而他的英文名字 Andrew (簡稱 Andy), 並無太多討論。 我立刻附議, 認為這是個好名字, 請看, 華人出色的理論計算機專家姚期智, 就用這個英文名子, 當然很佳。 在決定外孫的中英文名字之前, 女兒已經替他取了一個小名, 叫做小 $\pi$, 後來我們習慣叫他的小名, 常常記不起他的真正名字。 小 $\pi$ 這個名字是有來由的, 這個小傢伙在媽媽肚子裡就很好動, 有事沒事就踢媽媽, 「踢踢」寫成英文字母就是「TT」, 合起來就像希臘字母的 $\pi$。

不曉得這和我從小就講一些圓周率 $\pi$ 的故事給家人聽有沒有關係。 但我立刻想起來, 當大女兒小學 6 年級的時候, 我們全家到美國新澤西州的 Rutgers 大學訪問一年。 到了一個新地方, 一切都要從頭適應, 包括辦公室的電腦系統也是新的。 為了適應新電腦, 我試寫了幾個程式當作練習, 其中有一個程式就是利用反正切 arctan 的 Taylor 展開式, 計算圓周率 $\pi$ 到一千位小數。 記得當時還展示給女兒們看, 並且背出五十位給她們聽。 我有愛收集東西的習慣, 所以就到我的舊日檔案, 找到當年的程式 (如附錄)* * 附錄請見本刊網站: http://w3.math.sinica.edu.tw/mathmedia/default.jsp 以及印出來到一千位小數的 $\pi$, 如下所示。

3.14159 26535 89793 23846 26433 83279 50288 41971 69399 37510
58209 74944 59230 78164 06286 20899 86280 34825 34211 70679
82148 08651 32823 06647 09384 46095 50582 23172 53594 08128
48111 74502 84102 70193 85211 05559 64462 29489 54930 38196
44288 10975 66593 34461 28475 64823 37867 83165 27120 19091
45648 56692 34603 48610 45432 66482 13393 60726 02491 41273
72458 70066 06315 58817 48815 20920 96282 92540 91715 36436
78925 90360 01133 05305 48820 46652 13841 46951 94151 16094
33057 27036 57595 91953 09218 61173 81932 61179 31051 18548
07446 23799 62749 56735 18857 52724 89122 79381 83011 94912
98336 73362 44065 66430 86021 39494 63952 24737 19070 21798
60943 70277 05392 17176 29317 67523 84674 81846 76694 05132
00056 81271 45263 56082 77857 71342 75778 96091 73637 17872
14684 40901 22495 34301 46549 58537 10507 92279 68925 89235
42019 95611 21290 21960 86403 44181 59813 62977 47713 09960
51870 72113 49999 99837 29780 49951 05973 17328 16096 31859
50244 59455 34690 83026 42522 30825 33446 85035 26193 11881
71010 00313 78387 52886 58753 32083 81420 61717 76691 47303
59825 34904 28755 46873 11595 62863 88235 37875 93751 95778
18577 80532 17122 68066 13001 92787 66111 95909 21642 01989

我把這個檔案寄到美國, 算是我給未出生的小 $\pi$ 的第一件禮物。 2013 (去)年底, 我們到德州看小外孫, 回台灣以後, 還是經常透過視訊, 看著他一天天長大, 內心充滿喜悅。 3月14日, 女兒很高興的告訴我們, 美國把這一天定為「$\pi$日」, 鼓勵科學發展。 因此, 3月14日在我們家中變得與6月28日一樣重要了。 有關6月28日的故事, 等一下再說。

三、高中時期 --- 與 $\pi$ 有約

我與圓周率 $\pi$ 的淵源始自高中時期。 我出生在離台中市二十公里的南投縣小鎮草屯, 初中就讀台中市立一中, 高中就讀台中省立一中。 從鄉下到城裡讀書, 有許多要適應的地方。 例如, 初中剛開始, 一口台灣國語, 被同學笑了很久才慢慢改進。 高中的時候, 我和坐在隔壁的同學倒是相處得很好, 在我這個鄉巴佬眼中, 他算是博學多聞, 常常告訴一些我聞所未聞的常識。

有一次他告訴我:「你知道嗎? 我哥哥說美國有很多數學天才。」 老實說, 除了地理和歷史課本的一些印象之外, 我對美國所知為零。

「是嗎? 那他們有多厲害?」 談到數學我就有興趣了。

「譬如說, 他們能背圓周率 $\pi$ 到五十位小數。」 從小時候老師就教導我們圓周率是 3.14, 我也聽說過比較精確的 3.1416, 但是到五十位小數倒是不一樣。

「背到五十位有何困難? 我們在國文課都能背那麼長的文章了, 只有五十個阿拉伯數字一定很容易。」

我回家後翻箱倒櫃, 找到一本大哥的數學書, 從中找出五十位小數的圓周率。 接下來的一星期裡, 我一有空就背這個數值, 直到隨心所欲為止。 我於是到學校將這五十位數值背誦給我隔壁的同學聽, 不過他還是認為我不是天才。 只是, 從那時候開始, 我就對 $\pi$ 產生了很大的感情。 在後來的日子裡, 我有時就會述說這一段故事給女兒們聽, 並且背 $\pi$ 娛樂她們。

四、大學時期 --- 杜詩統老師的教導

我大學聯考填志願是把各學校的物理系、化學系、數學系, 依照前一年的排行榜次序填志願, 結果考上中央大學物理系。 中央大學當時剛剛在台灣復校不久, 師資缺乏, 拜託了中原大學杜詩統老師來兼課, 教物理系的微積分和數學系的高等微積分。 物理系的微積分課本是 Johnson 第四版, 杜老師是一位很稱職的老師, 我們學得很紮實, 到現在我都還記得當時他把微積分基本定理很確實地證明的過程。

不過杜老師給我最大的印象是, 他第一天上課就告訴我們, 他要在高等微積分的課裡證明 $\pi$ 是無理數。 他並且說, 其實我們只要學會一些積分公式就可以做這件事情了。 我在高中時學過如何證明 ${\sqrt 2}$ 是無理數, 也聽說過 $\pi$ 是無理數, 但是從來不曾想像要怎樣證明這件事。 因為對 $\pi$ 的特殊感情, 我很急著想知道如何證明它是無理數。 下了課, 我就到辦公室找杜老師, 問他證明 $\pi$ 是無理數的方法。

「你會算 $\sin$ 和 $\cos$ 的積分嗎?」他試探的問我。

「報告老師, 我還不知道積分是什麼。」他看了看我, 大概覺得我是個無厘頭的小孩, 這種程度還敢來問問題。 不過他很有耐心, 寫了一個積分式子, 要我回家算清楚之後再去找他。

我回家之後, 把微積分課本直接翻到三角函數積分的地方, 在還沒有讀積分定義的情況下, 將三角函數積分的公式做了整體操作的瞭解, 把老師出的習題做出來。 我很有信心的又去找老師, 這次他稍微覺得孺子可教, 又給了我第二道習題。 如此來回了幾趟, 我越來越會做積分的操作, 最後把老師出的幾道習題合起來, 很神奇的居然證明了 $\pi$ 是無理數。

也許是這個機緣, 加上後來讀了王九逵教授翻譯的一篇「幾何與直觀」的文章, 深覺有趣, 二年級我就轉到數學系就讀。

五、電子計算機算 $\pi$ 的威力

我雖然轉到數學系, 但還是常常回物理系找同學聊天。 有一次, 聽說物理系要找附近電信研究所的一位研究員來開授計算機程式語言的課。 當時中央大學還沒有計算機, 所以我也好奇的去選這門課。 因為老師要上班的關係, 我們上課的時間是在每星期三晚上的三個小時。 我們學的是 IBM 1130 FORTRAN, 那是一個相當古老的編譯器, 例如, 正整數算到超過 2048 就會變成負整數。 我們的第一個程式作業是要解一元二次多項式的兩個根。 我們把程式寫在 coding sheet 上, 老師帶回去請人幫忙打卡, 因為每一行程式要打在一張卡片上, 一個程式就是一疊卡片。 隔一個禮拜, 老師帶來卡片以及印出程式和答案的報表紙, 大部份人的結果都是程式有誤, 於是在報表紙上修改程式, 老師再帶回去修改卡片, 重新跑程式。 有些人來來回回, 改了許多次才改對。 一學期下來, 總共才寫了兩個程式, 感覺很沒有成就感。

後來我們打聽到, 台大電機系開授暑期電腦班, 我們幾位同學就相約報名上課去。 我和兩位物理系同學租屋住在新生南路的巷子, 他們去學英文和打字, 我到台大上電腦課, 過了一個充實的暑假。

電腦班教好幾門課, 讓我們學到許多新東西。 最讓人高興的是, 我們可以自己打卡片, 然後在程式的卡片組最上頭放上一張「控制卡」, 就可以拿到計算機中心, 請他們上電腦跑程式, 兩天後就可以看到答案。我很快就做完所有程式習題, 剩下的時間就寫一些自己認為有趣的程式, 例如, 算 $e$ 和 $\pi$ 到一千位小數等。

這是我第一次感受到電腦的威力, 深深覺得, 它簡直是為數學家設計的最佳機器。 只是, 後來我在數學界發現, 數學家幾乎是最後會使用電腦的一批人。

計算圓周率的基本道理是 $\dfrac{\pi}4=\arctan 1$, 並且利用 $\arctan x$ 的泰勒展開式 $$\arctan x=x-\frac{x^3}3+\frac{x^7}7-\frac{x^9}9+\frac{x^{11}}{11}-\cdots\hbox{。}$$ 但是如果直接將 $x=1$ 帶入上式, 要求到比較多準確位數, 一定要算非常多項, 這是一個不智的作法。 反過來, 如果能用一個比較小的 $x$ 得到一樣的效果, 計算量將會大大降低。 這件事情可以靠下面的公式達到: $$\arctan x+\arctan y=\arctan \frac{x+y}{1-xy}\hbox{。}$$ 利用這個公式, 我們可以得到 $$2 \arctan \frac 15=\arctan \frac 5{12},$$ 再利用同樣公式, 進一步得到 $$4 \arctan \frac 15=2 \arctan \frac 5{12}=\arctan \frac{120}{119},$$ 其中 $\dfrac{120}{119}$ 比 1 略大, 事實上, $$\arctan\frac{120}{119}-\arctan 1=\arctan\frac{120}{119}+\arctan(-1)=\arctan \frac 1{239},$$ 因此 $$\frac {\pi}4=\arctan 1=4 \arctan \frac 15-\arctan \frac 1{239}\hbox{。}$$ 利用 $\arctan x$ 的泰勒展開式來算 $\arctan \frac 15$ 和 $\arctan \frac 1{239}$, 算比較少項就可以得到比較高的準確度。 這就是著名的 Machin 公式 (1706年), 我的程式就是利用這個公式計算 $\pi$。 後來有許多人推導出利用更小的 $x$ 計算 $\pi$ 的各種公式, 舉例如下。 \begin{eqnarray*} \frac {\pi}{4}&=&22 \arctan \frac{24478}{873121}+17 \arctan\frac{685601}{69049993}\hbox{。}\\ \frac {\pi}{4}&=&12 \arctan \frac 1{49}+32 \arctan \frac 1{57}-5 \arctan \frac 1{239} +12 \arctan \frac 1{110443}\hbox{。 (Takano, 1982)}\\ \frac {\pi}{4}&=&44 \arctan \frac 1{57}+7 \arctan \frac 1{239}-12 \arctan \frac 1{682} +24 \arctan \frac 1{12943}\hbox{。 (Stormer, 1896)}\\ \frac {\pi}{4}&=&183 \arctan \frac 1{239}+32 \arctan \frac 1{1023}-68 \arctan \frac 1{5832}+12 \arctan \frac 1{110443}\\ &&-12 \arctan \frac 1{4841182}-100 \arctan \frac 1{6826318}\hbox{。 (黃見利, 1994)} \end{eqnarray*}

新一代計算圓周率 $\pi$ 採取收斂更快速的公式。 這一類的公式, 始於印度天才數學家 Ramanujan 在 1914 年發現的公式 $$\frac 1\pi=\frac{2\sqrt{2}}{9801}\sum_{k=0}^\infty\frac{(4k)!(1103+26390k)}{(k!)^4 (396^{4k})}\hbox{。}$$ 到了 1987 年, Chudnovsky 兄弟得到下列公式, 那就是最近用來計算圓周率 $\pi$ 到幾兆位小數的公式。 電腦協助數學計算, 到此可以說已經達到了神奇的境界。 $$\frac 1\pi=\frac{12}{640320^{3/2}} \sum_{k=0}^\infty \frac{(6k)!(13591409+545140134k)}{(3k)!(k!)^3(-640320)^{3k}}\hbox{。}$$

六、與君有約 --- 證明無理數之爭

物換星移, 大學畢業忽忽過了四十年。 在台灣的同學利用五月份的中央大學校慶, 到雙連坡舉辦四十周年同學會。 在美國的同學則利用7月4日美國國慶, 到舊金山舉辦同學會, 我和內人代表台灣的同學們到美國參加聚會。 大女兒全家特意從德州飛來加入我們, 好讓我們看看小 $\pi$。 相聚的第一天晚上, 大家互相介紹多年來的發展, 我也介紹了小 $\pi$ 名字的由來。

我們班上的同學, 除了兩個人在數學系任教之外, 轉行到資訊界的人很多, 其他各種行業都有, 所以聚會時很少談數學。 倒是才華洋溢的 J, 他除了當銀行的 CEO 遊刃有餘之外, 下班時, 還研究宗教、哲學、相對論等等; 精力充沛的他, 經常夫婦倆人到世界各地旅遊。 這次, 他如往常一樣展現出他的博學多聞, 這邊跟學佛的同學大談心經, 一下又到另一處談康德哲學、儒家道理。

對於 $\pi$, 他也有興趣談談無理數。 我跟他說當年杜老師教我證明 $\pi$ 是無理數的故事, 但誠實的招認, 當時一定沒有瞭解透徹, 現在已經想不起當時是如何證明的了。 我順便向大家複習了 $\pi$ 是超越數, 以及 $\sqrt 2$ 雖然也是無理數, 但卻只是代數數。 J 才思敏捷, 雖然是幾十年前學的知識, 他還能立刻展現如何證明 $\sqrt 2$ 是無理數。 我們陷入下面的爭論, 我認為他會證明 $\sqrt 2$ 是無理數, 那是因為以前學過, 他則以為, 只要思想透徹、 抓到重點, 不管以前學過與否, 都能證明出來。 這話也有相當道裡, 事實上, J 證明 $\sqrt 2$ 是無理數的方法, 與一般的證明方法略有不同, 有可能是抓到重點, 自己證出來的。

一般證明 $\sqrt 2$ 是無理數的方法是利用反證法。假設 $\sqrt 2=a/b$, 其中 $a$ 和 $b$ 是兩個互質的正整數, 將式子平方並整理得到 $a^2=2b^2$ 是偶數, 所以 $a$ 也是偶數 (不然的話, 如果有整數 $c$ 使得 $a=2c+1$, 則 $a^2=2(2c^2+2c)+1$ 是奇數, 矛盾), 也就是有整數 $d$ 使得 $a=2d$, 這樣就有 $4d^2=2b^2$, 也就是 $b^2=2d^2$, 如前相同論證, $b$ 也是偶數, 所以 $a$ 和 $b$ 這兩個數有一個大於 1 的公因數, 矛盾。

J 證明 $\sqrt 2$ 是無理數的方法如下。 假設 $\sqrt 2=a/b$, 其中 $a$ 和 $b$ 是正整數, 將式子平方並整理得到 $a^2=2b^2$, 式子左邊因數分解有 $2^{2m}$ 這一項, 式子右邊因數分解有 $2^{2n+1}$ 這一項, 但是 $2m\not=2n+1$ (偶數不會是奇數), 矛盾。

要註明的是, J 的證明看起來更簡短, 這是因為它利用了因數分解的唯一性, 而其證明是要花一些功夫的。

為了測試他「只要思想透徹、抓到重點, 不管以前學過與否, 都能證明出來」的論點, 我給他一道習題, 要他證明 $$e=\frac 1{0!}+\frac 1{1!}+\frac 1{2!}+\frac 1{3!}+\cdots$$ 是無理數, 這個證明的難度介於證明「$\sqrt 2$ 是無理數」與證明「$\pi$ 是無理數」之間, 應該比較靠近前者。 不過我猜想這應該是J沒有見過的問題, 用來測試他算是恰當。 我給他一個月的時間證明, 但是他十分不以為然, 精明的要求我證明 $\pi$ 是無理數, 與他競賽, 看誰先證明出來, 而且規定, 回家以後才能開始想。

七、證明 $\pi$ 是無理數

J 很快就寄來 $e$ 是無理數的證明, 展現了他的思慮敏捷。 如果 $e$ 是有理數, 設 $e=a/b$, 其中 $a$ 和 $b\gt 1$ 為正整數, 則 $b!e$ 是整數。 但是, $b!e=$整數+餘項, 其中餘項為 $1/(b+1)+1/(b+1)(b+2)+1/(b+1)(b+2)(b+3)+\cdots$, 但是餘項是一個比首項及公比都是 $1/(b+1)$ 的等比級數小的正數, 此等比級數和為 $1/b\lt 1$, 這和 $b!e$ 是整數矛盾。 因此 $e$ 是無理數。

當時我還在調適時差中, 看起來要自己證明 $\pi$ 是無理數應該有困難, 於是上網搜尋尋找解救, 發現很多證明 $\pi$ 是無理數的方法。 下面是 Niven 證明的簡化。 假設 $\pi$ 是有理數, 也就是 $\pi=a/b$, 其中 $a$ 和 $b$ 為正整數。 對於任意正整數 $n$, 定義多項式函數 $$f(x)=x^n (a-bx)^n/n!, \qquad x \in {\bf R};$$ 以及其前面 $n$ 個偶次微分的交錯和 $$F(x)=f(x)-f^{(2)}(x)+\cdots+(-1)^i f^{(2i)}(x)+\cdots+(-1)^{(2n)} f^{(2n)}(x)\hbox{。}$$

我們先證明「$F(0)+F(\pi)$ 是整數」。首先, 因為 $f(x)=\sum_{k=n}^{2n}c_k x^k$, 其中 $c_k$ 是整數, 所以, 當 $k\lt n$ 時 $f^{(k)}(0)=0$, 當 $k\ge n$ 時 $f^{(k)}(0)=c_k k!/n!$ 是整數, 因此 $F(0)$ 是整數。 其次, 因為 $f(x)=f(\pi-x)$, 所以利用微分的連鎖率 $k$ 次得知 $f^{(k)}(x)=(-1)^k f^{(k)}(\pi-x)$, 代入 $x=\pi$ 得到 $f^{(k)}(\pi)=(-1)^k f^{(k)}(0)$ 也是整數, 因此 $F(\pi)$ 是整數。 這就證明了「$F(0)+F(\pi)$ 是整數」。

再來, 我們要證明「$\int_0^\pi f(x)\sin x dx=F(0)+F(\pi)$」。 因為 $f^{(2n+2)}$ 是零多項式, 所以 $F''+F=f$, 又由於 $\sin'=\cos$ 以及 $\cos'=-\sin$, 經由微分的乘法率可知 $(F'\cdot \sin-F\cdot \cos)'=f\cdot \sin$, 由微積分基本定理可以得到 $$\int_0^\pi f(x)\sin xdx=(F'(x)\sin x-F(x)\cos x)\Big|_0^\pi =F(0)+F(\pi)\hbox{。}$$

有了前面兩段的結論, 再加上, 當 $0\lt x\lt \pi$ 時 $f(x)\gt 0$ 且 $\sin x\gt 0$, 可以知道 $\int_0^\pi f(x)\sin x dx=F(0)+F(\pi)$ 是一個正整數。 另一方面, 當 $0\le x\le \pi$ 時 $0\le x(a-bx)\le \pi a$ 且 $0\le \sin x\le 1$, 所以 $\int_0^\pi f(x)\sin x dx \le \pi(\pi a)^n/n!$。 但是, 當 $n$ 很大時 $\pi(\pi a)^n/n!\lt 1$, 這和 $\int_0^\pi f(x)\sin x dx$ 是正整數矛盾。

八、中場 --- 神秘的數字

接下來, 我們來說說6月28日的故事。 這可以從數字的故事講起。

數字的神秘性產生於原始時代。 在人類的蒙昧時期, 由於對自然、 世界的無知和不可把握, 人們往往在認知外界事物的時候, 把「數」看得比「質」更重要, 也把「數」看作人神溝通的重要途徑。

1 是個既簡單又複雜的數字, 它既是所有數中最小的一個, 又是其中最有包容性的一個, 所有的事物綜合之後可以劃歸為一個整體。

2 的神秘反映在二元對立的思維上, 由此衍生出陰陽對立, 人們意識到自然界許多的對立, 例如: 晝和夜、晴和雨、高和低、冷和暖、雄和雌、生和死等等。

3 在符號象徵體系中, 幾乎不含任何反面意義, 例如, 中國民間故事中有「三度走運」之說, 基督教三位一體的教義使得基督教唯一的真神能通過聖父聖子與聖靈接受眾人的頂禮膜拜。

世界各地都有關於 4 的說法。 美洲人認為, 四是宇宙的核心原則。 瑪雅人有四種顏色和自然年的四種記法。 在阿茲台克人的宇宙意識中, 四株世界樹支撐天穹。 「四方位據說是風的發源地, 懸著四個巨大水罐, 倒出的水就是雨」。

在中國的神秘文化中, 「五行」分別處於東、西、南、北、中五個方位, 又有各自的金、木、水、火、土五種元素。 在西方, 當五角星形中兩點朝下、一點朝上的時候, 便可以看作直立的人形, 有頭、胳膊和腿。 在聖經《新約》中, 耶穌用五個餅就餵飽了五千人。

6 是一個有趣的數字。 中國人說六六大順。 西方人說, 上帝在六天裡創造世界後, 第七天就安息了。 聖奧古斯丁認為 6 含義非同尋常, 因為它是頭三個數字 1、 2、 3 的和。 福音書中描寫的六次寬恕構成了基督教傳統中與六有關的少數幾個象徵系列之一。 一個包括「六」的重要視覺標誌是由兩個等邊三角形反向疊成的六角星, 即大衛之盾。

古希臘時代, 畢氏學派的數學家用數學的語言來談論 6。 他們定義完全數 (perfect number, 事實上稱呼完美數也許較恰當, 不過約定俗成, 我們仍採用完全數) 為除了自身以外的所有因數 (稱為真因數) 的和恰好等於它本身的數。 例如: 第一個完全數是 6, 它有因數 1、 2、 3、 6, 除去 6 本身之外, 其餘的 3 個數相加, $1+2+3=6$, 恰好等於本身。 第二個完全數是28, 它有因數 1、 2、 4、 7、 14、 28, 除去 28本身之外, 其餘的 5 個數相加, $1+2+4+7+14=28$, 也恰好等於本身。 第二個完全數的神祕性在於, 月亮繞地球一周約為28天。 接下來的完全數是 496, 再接下來則是 8128。

九、計算完全數 --- 電腦再顯威力

我對完全數的喜好也是從高中時代開始。那時候, 我是一個不聽話的學生, 不喜歡數學老師的教課方式, 上課時自己讀課本、 看一些課外數學書, 有一次讀到一本叫《大眾數學》的通俗讀物, 介紹了完全數, 其中最引發我興趣的是, 關於偶完全數的刻劃, 但是書中並未給出證明。 經過一番苦思, 我終於獨立完成證明。 後來還寫了一篇文章投到剛創刊的科學月刊, 到我上中央大學時正式刊登出來。 從此我就和「完全數」結下了終生的緣分。

定理1: 若 $2^{n+1}-1$ 是質數, 則 $2^n (2^{n+1}-1)$ 是完全數。 而且, 偶完全數必定形如 $2^n (2^{n+1}-1)$, 其中 $2^{n+1}-1$ 是質數。

證明: 假設質數 $p=2^{n+1}-1$, 則 $2^n p$ 的真因數有 1, 2, $\ldots,2^n,p,2p,\ldots,2^{n-1}p$, 由於 $1+2+\cdots+2^n=2^{n+1}-1=p$ 且 $1+2+\cdots+2^{n-1}=2^n-1$, 所以這些因數的總和為 $p+p(2^n-1)=2^n p=2^n (2^{n+1}-1)$。 因此, $2^n (2^{n+1}-1)$ 是完全數。

(以上為歐基里得的證明。)

考慮偶完全數 $\alpha=2^n p$, 其中 $n$ 為正整數, $p$ 為奇數。 $\alpha$ 的所有因數和為 $$S=(2^n\ \hbox{的所有因數和}) (p \ \hbox{的所有因數和})=(2^{n+1}-1)(p+d),$$ 其中 $d$ 是 $p$ 的所有真因數的總和。 因為 $\alpha$ 是完全數, 所以 $\alpha=S-\alpha$, 由此可以知道 $(2^{n+1}-1)(p+d)=S=2\alpha=2^{n+1} p$, 化簡可以得到 $p=(2^{n+1}-1)d$, 其中 $2^{n+1}-1\gt 1$。 如果 $d\gt 1$, 則 $d$ 和 1 是 $p$ 的兩個不同的真因數, 這和 $d$ 是 $p$ 的所有真因數的總和矛盾, 所以 $d=1$, 由此可以知道 $p=2^{n+1}-1$ 為質數, 而偶完全數 $\alpha=2^n (2^{n+1}-1)$。

(以上為尤拉的證明。)

尋找偶完全數的工作, 到此變成決定 $M_m=2^m-1$ 是否為質數的工作, 尋求這種型式的質數, 最早以法國數學家莫仙尼 (Marin Mersenne, 1588$\sim$1648)算出最多, 故我們用他的名字稱呼此等質數, 叫莫仙尼質數。

要找尋形如 $2^m-1$ 的質數, 首先 $m$ 先要是質數才可以。 因為, 如果 $m$ 不是質數的話, 一定存在兩個大於 1 的整數 $a$ 和 $b$ 使得 $m=ab$, 因此 $$M_m=2^m-1=(2^a)^b-1=(2^a-1)[(2^a)^{b-1}+(2^a)^{b-2}+\cdots+2^a+1],$$ 是兩個大於1的整數的積, 必定為合成數。所以, 要尋找莫仙尼質數 $M_m$ 時, 只需要考慮 $m$ 是質數的情形。 但是這個命題反過來並不成立, 也就是說, $m$ 是質數並不保證 $M_m$ 也是質數, 例如, 11 是質數, 但是 $M_{11}=2^{11}-1=23\times 89$ 卻不是質數。 這是 $m$ 在 20 之前的唯一例外, 其他 $M_2=3$、 $M_3=7$、 $M_5=31$、 $M_7=127$、 $M_{13}=8191$、 $M_{17}=131071$、 $M_{19}=524287$ 都是質數。

當 $m$ 越來越大時, $M_m$ 的位數快速增大, 驗證其是否為質數的工作越來越不容易。 由「數論」的知識, 我們知道「合成數 $a$ 必定有小於或等於 $\sqrt a$ 的質因數」。 這提供我們檢查某一個數是否為質數的方法:用小於或等於 $\sqrt a$ 的所有質數去試驗, 如果這些數都不是 $a$ 的因數, 那麼 $a$ 就是質數。 但是, 使用這種方法, 要算出 $M_{19}$ 是質數, 已經十分不容易, 於是又有下面的定理。

定理2: 若 $p$ 為質數, 則 $M_p$ 的質因數必定形如 $ap+1$, 其中 $a$ 為正整數。

證明: 假設 $M_p=2^p-1$ 有一質因數 $x=ap+b$, 其中 $0\lt b\lt p$。 因為 $2^p\equiv 1 (\mod x)$, 所以 $$2^{x-1}\equiv 2^{ap+b-1}\equiv (2^p)^a\cdot 2^{b-1}\equiv 2^{b-1} (\mod x),$$ 根據費瑪小定理 $2^{x-1}\equiv 1 (\mod x)$, 所以 $2^{b-1}\equiv 1 (\mod x)$。 如果 $b\gt 1$, 則因為 $p\gt b\gt b-1\gt 0$, 所以 $p$ 和 $b-1$ 互質。 於是, 存在兩個正整數 $\alpha$ 和 $\beta$ 使得 $$\alpha p=\beta(b-1)+1\quad \hbox{或是}\quad \alpha(b-1)=\beta p+1\hbox{。}$$ 當 $\alpha p=\beta(b-1)+1$ 時, $$1\equiv(2^p)^\alpha\equiv 2^{\alpha p}\equiv 2^{\beta(b-1)+1}\equiv (2^{b-1})^\beta\cdot 2\equiv 2 (\mod x),$$ 得到矛盾。 同樣地, 當 $\alpha(b-1)=\beta p+1$ 時也矛盾。 因此 $b=1$, 證明完畢。 $\Box$

利用這個定理, 檢驗工作又減輕不少。 但是莫仙尼質數是呈幾何級數增大, 當 $p$ 越大, 檢驗 $M_p$ 是否為質數就越難。 尤拉在 1750 年算出 $M_{19}$ 的下一個莫仙尼質數 $M_{31}$, 早期的尋找工作就算告一段落。 1876 年法國數學家盧卡斯 (Lucas) 用筆算出一個最大紀錄 $M_{127}$, 這是一個長達 39 位的數, 是人類用紙筆算出的 12 個莫仙尼數中之最大者。 現在將這十二個質數列出如下 : 2、 3、 5、 7、 13、 17、 19、 31、 61、 89、 107、 127。

再下去的數越來越大, 不是人類有限的精力所能擔當得了的。 尤其莫仙尼質數的密度越來越小, 要找到確實不容易。 自從人類發明電腦以後, 計算能力逐日改良, 使得這項工作得以繼續下去。 最新的資料顯示, 到目前為止, 人們找到 48 個莫仙尼數, 最大的一個是 $M_{57,885,161}$, 高達 17,425,170 位。 最近的 14 個莫仙尼數都是透過從 1996 年起的一個計畫 GIMPS (Great Internet Mersenne Prime Search), 利用下列的盧卡斯 --- 來默判定法 (Lucas-Lehmer primality test) 算出來的。 數學加上電腦的威力, 再度展現無遺。

Lucas-Lehmer定理: 如果 $s_0=4$ 且當 $i\ge 1$ 時 $s_i=s_{i-1}^2-2$, 則對於任意質數 $p$, $M_p$ 是質數的充分必要條件是 $s_{p-2}\equiv 0 (\mod M_p)$。

十、幕落 --- 小女兒的失落

由於對完全數的喜好, 當女兒們小的時候, 我也常常讓她們「玩」完全數。 當她們還只會加減法時, 我教她們「玩」乘法和除法, 我讓她們玩的方法, 並不是背九九乘法表, 而是利用累加去實現乘法、 累減去實現除法。 為了增加趣味性, 我順便告訴她們因數的定義, 開始計算第一個完全數。 當時她們確實算出第一個完全數 6, 我也告訴她們一些數字的故事。 等她們大一點以後, 我有時也會告訴她們更大的完全數。

大女兒很受感染, 深深愛上完全數。 她參加交通大學女籃校隊時, 選的球衣號碼就是 6 號。 2008 年大女兒的結婚日期堅持選在 6 月 28 日 (星期六), 這一個日子就成為我們家中的一個重要節日了。 而 3 月 14 日是因為小 $\pi$ 的關係, 成為我們家中的另一個重要節日。 有一次, 我告訴同事 C 君有關我們家中的這些故事, 他一眼就看出這兩組數字的關係, 套用他郵件的話:

6 and 28 are perfect numbers.

6 月 28 日 = 2$\pi$ 圓滿

Amazing that she can find such perfect and 圓滿 date!!

至此, 這兩組重要的日子找到了它們的關係。 小女兒預定明年結婚, 她就想在 3 月 14 日結婚, 喜帖的詞也想好了 : 「兩個一半圓滿, 湊成一個圓滿。」可惜, 男朋友的媽媽經過算命, 認為這一天不是好日子, 兩個比較好的日子之一是 3 月 1 日。 小女兒傷心了一陣子, 只好選擇 3 月 1 日。 經過我努力分析, 3.1 雖然沒有 3.14 準確, 但是比起古時候用 3 估計 $\pi$, 已經正確多了, 而她也接受這個新日期了。

參考文獻

1.張鎮華, 完全數與莫仙尼質數, 科學月刊第二期第三卷, 民國60年。

2.張鎮華, 閒話尤拉的絕招, 數學傳播, 第2卷第3期, 民國67年2月, 第42-45頁。

3.王九逵, 怎麼算 $\pi$, 數學傳播, 第10卷第2期, 民國75年6月, 第2-9頁。

4.余文卿、葉國榮, 關於圓周率 $\pi$, 數學傳播, 第13卷第3期, 民國78年9月, 第55-59頁。

5.黃文璋, 完全數與梅仙尼質數, 數學傳播, 第21卷第3期, 民國86年9月, 第14-30頁。

6.文光耀、江邵祥, 試算表與完全數的探究, 數學傳播, 第22卷第3期, 民國87年9月, 第77-82頁。

7.石厚高, 祖沖之計算圓周率之迷, 數學傳播, 第25卷第1期, 民國90年3月, 第82-86頁。

8.黃見利, 如何衍生Machin型公式, 數學傳播, 第29卷第1期, 民國94年3月, 第18-29頁。

9.黃見利, 兩種近似圓周率的算法, 數學傳播, 第32卷第1期, 民國97年3月, 第77-83頁。

10.黃見利, 一個有趣的計算圓周率方法 --- 雙反正切項尤拉型公式, 數學傳播, 第33卷第1期, 民國98年3月, 第57-64頁。

---本文作者任教台灣大學數學系---

附錄

/** pi/4=4*arctan(1/5)-arctan(1/239)  **/
/** arctan(x)=x-x^3/3+x^5/5-x^7/7-...  **/
#include "stdio.h"
#define n 1010
#define nn 1000

main( )
{
  int a[n],pi[n],i,j,k;

  arctan(5,pi);
  arctan(239,a);
  mul(pi,4);
  sub(pi,a);
  mul(pi,4);
  output(pi);
}
arctan(xi,a)
int xi,a[n];
{
  int xp[n],xbeg,xpsign,xi2,i,t[n],j;
  xi2=xi*xi;
  xp[0]=1;
  xbeg=0;
  for (i=1; i\lt =n; i++)
      xp[i]=0;
  div(xp,xi,xbeg);
  for (i=0;i\lt =n;i++)
      a[i]=xp[i];
  j=3;
  xpsign=-1;
  while (xbeg\lt =n)
    {  div(xp,xi2,xbeg);
       xbeg=begin(xp,xbeg);
       divt(t,xp,j,xbeg);
       if (xpsign==1)
           add(a,t);
       else
           sub(a,t);
       xpsign=-xpsign;
       j=j+2;
     }
}
add(a,b)
int a[n],b[n];
{
  int i,r;
  r=0;
  for (i=n; i\gt =0; i--)
    { a[i]=a[i]+b[i]+r;
      if (a[i]\gt 9)
	{ a[i]=a[i]-10; r=1; }
      else
        r=0;
    }
}
sub(a,b)
int a[n],b[n];
{
  int i,r;
  r=0;
  for (i=n; i\gt =0; i--)
    { r=b[i]+r;
      if (a[i]\lt r)
        { a[i]=10+a[i]-r; r=1; }
      else
        { a[i]=a[i]-r; r=0; }
    }
}
mul(a,b)
int a[n],b;
{
  int i,r;
  r=0;
  for (i=n; i\gt =0; i--)
    { a[i]=a[i]*b+r;
      r=a[i]/10;
      a[i]=a[i]%10;
    }
}
div(a,b,beg)
int a[n],b,beg;
{
  int i,r;
  r=0;
  for (i=beg; i\lt =n; i++)
    { a[i]=a[i]+r;
      r=a[i]%b*10;
      a[i]=a[i]/b;
    }
}
begin(a,beg)
int a[n],beg;
{
  while ((beg\lt =n) && (a[beg]==0))
    beg++;
  return(beg);
}
divt(t,a,b,beg)
int t[n],a[n],b,beg;
{
  int i,r;
  r=0;
  for (i=0;i\lt beg;i++)
    t[i]=0;
  for (i=beg; i\lt =n; i++)
    { t[i]=a[i]+r;
      r=t[i]%b*10;
      t[i]=t[i]/b;
    }
}
output(a)
int a[n];
{ int i,j,k;
  printf("3.");
  j=1;  k=1;
  for (i=1;i\lt =nn;i++)
    {  printf("%d",a[i]);
       if (j==5)
	 { j=1; printf(" "); }
       else
         j++;
       if (k==50)
         { k=1; printf("\n  "); }
       else
         k++;
    }
  printf("\n");
}