一個有趣的計算圓周率方法

一個有趣的計算圓周率方法

一一雙反正切項尤拉型公式

黃見利

本文謹獻給作者的恩師 Dar´ıo Castellanos 教授 (生於1937年12月4日; 逝於1995年11 月 23日)。 作者永遠懷念他的教誨。

一 . 前言

在 1671 年, 蘇格蘭數學家 James Gregory 發現了下列的反正切型級數 tan⁻¹x=

∞

X

n=0

(−1)ⁿx²ⁿ⁺¹

2n + 1 , −1 ≤ x ≤ 1.

當 x = 1 時, 另一位偉大的數學家, 也是微積分的創造者之一, Gottfried Wilhelm Leibniz 在三年後的 1674 年發現了所謂的 Leibniz 級數 π

4 = 1 − 1 3+ 1

5 −1 7 + 1

9− 1

11+ · · · “不 過, 一位反正切型級數的發現者, 而且還著手於圓周率超越性方面的研究, 會忽略了在他的級數 中替代 x = 1 這麼如此明顯的舉手之勞, 還真是令人不可思議的一件事。 較可能的情況似乎是, 他不認為替代 x = 1 在他的級數中是重要的, 因為它的收斂速度實在太慢了, 因而無法使用在 圓周率的計算上。”^[1] 如果他真的如此考慮過, 那真是絕對正確且嚴謹周密的想法! 現在就讓我 們來瞧一瞧這條級數的收斂速度到底有多慢。

首先, 我們知道 Leibniz 級數是一種所謂的交替級數 (alternating series)。 因此, 我們可 以利用下列著名的定理來估計誤差:

交替級數估計定理 (Alternating Series Estimation Theorem) 假設 s =

∞

P

n=1

(−1)ⁿ⁻¹a_n 為一交替級數之和, 而 s_n =

n

P

k=1

(−1)^k−1a_k 為其部份項之和, 且滿足下列兩條件:

57

(i) 0 ≤ an+1 ≤ an

(ii) lim

n→∞a_n= 0

=⇒ |Rn| = |s − sn| ≤ an+1.

由這個定理, 我們得到了一個令人驚異的事實。 假設我們現在想得到德國數學家 Ludolph van Ceulen 在 1615 年死後才發表的 35 位圓周率小數

3.14159265358979323846264338327950288。 由於第 36 位小數是4, 所以 |Rn| ≤ 25·10⁻³⁷ 不會影響到第 35 位小數的實際值。 所以 |Rn| ≤ 25 · 10⁻³⁷= 1

4 · 10³⁵, 亦即利用 Leibniz 級 數計算 π ≈ 4 − 4

3 +4 5 − 4

7 + 4 9− 4

11 + · · · + 4

16 · 10³⁵− 1 時, 絕對有 35 位圓周率小數 的正確值。 不過, 話說回來, 僅僅因為想求得 35 位圓周率小數位數, 我們利用 Leibniz 級數竟 然需計算 8 · 10³⁵ 項才能得到!!! 我們所要下的苦工比 Ludolph van Ceulen 當年好不到那裏 去, 甚至可能更糟! 因此, 在 1983 年之前, 幾乎所有圓周率小數位數的“狩獵者”使用了各種不 同的方法將 Gregory 級數加以修飾從而增加它的收斂速度, 並且尋求有效的反正切型公式和 其配合以計算圓周率小數位數。

附帶一提的是, Archimedes 利用多邊形逼進圓形以求得圓周率小數位數的方法在 Lu- dolph van Ceulen 手中達到最高頂點。 他的一生大部份在現今的荷蘭渡過, 而且是獻給了圓 周率小數位數的計算。 他的 35 位圓周率小數位數是死後才發表在其著作 De Arithmetische en Geometrische fondamenten(Leyden, 1615) 上。 德國人為了敬佩其毅力和成就, 因此將 圓周率稱為 Ludolph 數 (Ludolphsche Zahl)。

在 1983 年之前, 搜尋“好的”或“有效率的”反正切型公式是一件重要的數學工作。 但是, 何謂“好的”或“有效率的”反正切型公式? 在 1938 年, D. H. Lehmer 首先提出“測度值”的 概念來回答此問題^[2,3]。 有興趣的讀者可參閱本文作者的另一篇文章^[4]。 我們在此列出一些“好 的”或“有效率的”反正切型公式讓讀者驚奇一下:

π

4 = 83 tan⁻¹ 1 107



+ 17 tan⁻¹ 1 1710



− 22 tan⁻¹ 1 103697



− 24 tan⁻¹ 1 2513489



− 44 tan⁻¹ 1 18280007883



+ 12 tan⁻¹ 1

7939642926390344818

 + 22 tan⁻¹ 1

3054211727257704725384731479018

 . measure: 1.34085 (Wetherfield, 2004)

這是目前已知測度值最小的 Machin 型公式, 由本文作者的恩師 Wetherfield 在 2004 年 9 月 所發現。

π

4 = 183 tan⁻¹ 1 239

+ 32 tan⁻¹ 1 1023

− 68 tan⁻¹ 1 5832

+ 12 tan⁻¹ 1 113021



− 100 tan⁻¹ 1 6826318

− 12 tan⁻¹ 1 33366019650

 + 12 tan⁻¹ 1

43599522992503626068

. measure: 1.50840 (黃見利, 2003)

π

4 = 183 tan⁻¹ 1 239

+ 32 tan⁻¹ 1 1023

− 68 tan⁻¹ 1 5832

+ 12 tan⁻¹ 1 110443



−12 tan⁻¹ 1 4841182

− 100 tan⁻¹ 1 6826318

. measure: 1.51244 (黃見利, 1994)

π

4 = 83 tan⁻¹ 1 107

+ 17 tan⁻¹ 1 4443

+ 68 tan⁻¹ 1 11343

− 5 tan⁻¹ 1 110443



− 34 tan⁻¹ 1 595667

+ 5 tan⁻¹ 1 4841182

 . measure: 1.53457 (Wetherfield, 1995)

π

4 = 38 tan⁻¹ 1 47

− 14 tan⁻¹ 1 610

− 12 tan⁻¹ 1 247057

+ 6 tan⁻¹ 1 740943

 + 14 tan⁻¹ 1

15628852

− 7 tan⁻¹ 1 54608393

 . measure: 1.58112 (Koepp, 1994)

π

4 = 44 tan⁻¹ 1 57

+ 7 tan⁻¹ 1 239

− 12 tan⁻¹ 1 682

+ 24 tan⁻¹ 1 12943

 . measure: 1.58604 (St¨ormer, 1896)

π

4 = 22 tan⁻¹ 1 28

+ 2 tan⁻¹ 1 443

− 5 tan⁻¹ 1 1393

− 10 tan⁻¹ 1 11018

 . measure: 1.63435 (Escott, 1896)

π

4 = 44 tan⁻¹ 1 109

+ 95 tan⁻¹ 1 239

− 12 tan⁻¹ 1 682

+ 24 tan⁻¹ 1 12943



− 44 tan⁻¹ 1 6826318

 . measure: 1.65366 (Arndt, 1993)

二 . 一個有趣的新方法

使用

tan⁻¹ 1 a− b

= tan⁻¹1 a

+ tan⁻¹ b a²− ab + 1



, (4)

和

tan⁻¹1 a

= 2 tan⁻¹ 1 2a

− tan⁻¹ 1 4a³+ 3a

,

兩條反正切型公式, 並且適當地選擇 a 和 b 的值, 讓我們可以轉換已經存在了的公式或發現新 的公式以計算圓周率的數值。 不過, 這兩條公式僅是下列公式的變型而已:

tan⁻¹(z1) + tan⁻¹(z2) = tan⁻¹ z₁+ z2

1 − z1z₂



. (5)

當我們用 a = 2, b = 1 代入 (4), 便可得到 Euler 在 1738 年所發現的公式 π

4 = tan⁻¹1 2



+ tan⁻¹1 3



. (6)

如果將此公式和下列 Euler 在 1755 年所發現的反正切級數^[5] 一起配合使用:

tan⁻¹(x) =

∞

X

n=0

2²ⁿ(n!)² (2n + 1)!

x²ⁿ⁺¹

(1 + x²)ⁿ⁺¹, (7) 我們可得到

π 4 = 4

10

 1 + 2

3

 2 10

+2 · 4 3 · 5

 2 10

2

+ · · ·

 + 3

10

 1 + 2

3

 1 10

+2 · 4 3 · 5

 1 10

2

+ · · ·

 . 若以現代的電子計算機的“迴圈” (loop) 觀點來看, 這一條級數每做一個迴圈計算大約可 以增加 2

3 位圓周率小數位數。 其中最特別的是, 這一條級數的每一個小括弧中的分母僅僅出現 10!!! 這是由於 (6) 的關係使得 2 · (2²+ 1) = 10 和 3²+ 1 = 10。

重點來了。 我們有可能在其他的公式轉化成上述級數型式時也保留這個性質, 也就是每一 個小括弧中的分母僅僅出現 10 的次方嗎?

這個答案可以是肯定的, 只要符合下述情形。 我們現在將 (5) 重新寫成下列型式 tan⁻¹1

a

= tan⁻¹1 b

+ tan⁻¹ b− a ab+ 1



, (8)

則我們會發現到, 假使 a²+ 1 可以被 5^k 整除且 b²+ 1 可以被 5^s 整除, k 和 s 皆為正整數, 則公式

ab+ 1 b− a

2

+ 1 = (a²+ 1)(b² + 1) (b − a)² 顯示出ab+ 1

b− a

2

+ 1 可以被 5^k+s 整除! 更甚的是, 由公式 (8), 配合著公式 (7), 將產生如 同上述級數所具有的每一個小括弧中的分母僅僅出現 10 的次方這種特殊性質!

以下我們就開始用例子來說明這個奇妙的性質。 例如, 我們可以增進 (6) 的收斂速度, 只 需將 tan⁻¹1

2

替換成 tan⁻¹1 3

 加上另一個反正切值即可。 將 a = 2 和 b = 3 代入 (8), 我們得到

tan⁻¹1 2

= tan⁻¹1 3

+ tan⁻¹1 7

 ,

結合 (6) 後, 產生了另一條新公式 π

4 = 2 tan⁻¹1 3

+ tan⁻¹1 7



, (9)

C. Hutton 首先在 1776 年發現它, 並且發表在 Philosophical Transactions of the Royal Society。 約 70 年後的 1847 年, Thomas Clausen 利用此公式將圓周率計算至 248 位小數 位數並且發表在 Astronomischen Nachrichten。

公式 (9), 配合著公式 (7), 給出下列級數 π

4= 6 10

 1 + 2

3

 1 10

+ 2 · 4 3 · 5

 1 10

2

+ · · ·

 + 14

10²

 1 + 2

3

 2 10²

+ 2 · 4 3 · 5

 2 10²

2

+ · · ·

 . 現在, 可以注意到, 正如前文所預測的, 我們已經增進了收斂速度; 這一條級數每做一個迴圈計 算大約可以增加 1 位圓周率小數位數, 而且維持每一個小括弧中的分母僅僅出現 10 的次方這 種特殊性質!

為了增進 (9) 的收斂速度, 我們將 tan⁻¹1 3

 替換成 tan⁻¹1 7

 加上另一個反正切值, 亦即將 a = 3 和 b = 7 代入 (8), 得到

tan⁻¹1 3



= tan⁻¹1 7



+ tan⁻¹ 2 11



; 結合 (9) 後, 又產生了另一條新公式

π

4 = 2 tan⁻¹ 2 11

+ 3 tan⁻¹1 7



. (10)

配合著 (7), 此公式給出下列級數 π

4=352 10³

 1 +2

3

32 10³

+2 · 4 3 · 5

32 10³

2

+ · · ·

 + 42

10²

 1 +2

3

 2 10²

+2 · 4 3 · 5

 2 10²

2

+ · · ·

 .

這一條級數每做一個迴圈計算大約可以增加 11

2 位圓周率小數位數。

再一次, 令 a = 11

2 和 b = 7, 我們有 tan⁻¹ 2

11



= tan⁻¹1 7



+ tan⁻¹ 3 79



; 結合 (10) 後, 又產生了另一條新公式

π

4 = 5 tan⁻¹1 7

+ 2 tan⁻¹ 3 79



. (11)

這是 Euler 在 1779 年所發現的, 並發表在 Nova Acta Petropolitanae。

公式 (11), 配合著公式 (7), 給出下列級數 π

4= 7 10

 1 +2

3

 2 10²

+2 · 4 3 · 5

 2 10²

2

+ · · ·



+7, 584 10⁵

 1 +2

3

144 10⁵

+2 · 4 3 · 5

144 10⁵

2

+ · · ·

 .

每做一個迴圈計算大約可以增加 12

3 位圓周率小數位數。 利用此級數, Euler 竟然在一個小時 內算出了 20 位圓周率小數位數!

現在, 讓我們做一些簡化符號的手續. 令

χ(n) = 2²ⁿ(n!)² (2n + 1)!

C₁為第一個反正切項的係數, N₁ 為第一個反正切項的分子, D₁ 為第一個反正切項的分母, C₂ 為第二個反正切項的係數, N2 為第二個反正切項的分子, D2 為第二個反正切項的分母, 則任 何一條此種 Euler 型的公式皆可寫成下列型式

π

4 = C1tan⁻¹N₁ D₁

+ C2tan⁻¹N₂ D₂

;

配合著 (7), 我們得到了下列一般式 π

4 = C₁N₁D₁ N₁² + D²₁

∞

X

n=0

χ(n) N₁² N₁² + D²₁

n

+ C₂N₂D₂ N₂²+ D₂²

∞

X

n=0

χ(n) N₂² N₂²+ D₂²

n

.

由公式 (8) 我們知道 N₁² + D²₁ = 5^k 或 N₁² + D²₁ = 2 · 5^k, 以及 N₂² + D₂² = 5^s 或 N₂²+ D₂² = 2 · 5^s, 而且 k = C2, s = C1。 因此, 給出任何一條 Euler 型的公式, 相關的級數 都可以不用再經過上述那些繁瑣的計算而馬上得到!!!

緊接著 (11), 我們有 π

4 = 5 tan⁻¹ 29 278

+ 7 tan⁻¹ 3 79

 , 相關的級數為

π

4 = 515, 968 10⁶

∞

X

n=0

χ(n)2⁷· 29² 10⁷

n

+ 26, 544 10⁵

∞

X

n=0

χ(n)2⁴· 3² 10⁵

n

,

每做一個迴圈計算大約可以增加 2 位圓周率小數位數。

再下來, 我們有 π

4 = 5 tan⁻¹ 1, 457 22, 049

+ 12 tan⁻¹3 79

 ,

相關的級數為 π

4 = 32, 896, 402, 432 10¹¹

∞

X

n=0

χ(n)2¹¹· 1, 457² 10¹²

n

+45, 504 10⁵

∞

X

n=0

χ(n)2⁴· 3² 10⁵

n

,

每做一個迴圈計算大約可以增加 21

3 位圓周率小數位數。

再來一次, 我們有 π

4 = 17 tan⁻¹ 3 79



+ 5 tan⁻¹ 24, 478 873, 121

 , 相關的級數為

π

4=64, 464 10⁵

∞

X

n=0

χ(n)2⁴· 3² 10⁵

n

+24, 478 · 873, 121 ·2¹⁶ 10¹⁶

∞

X

n=0

χ(n)2¹⁹· 12, 239² 10¹⁷

n

, (12) 每做一個迴圈計算大約可以增加 3 位圓周率小數位數。

這個過程可以無止境地繼續下去而得到這種只有兩個反正切項的 Euler 型公式, 而且, 在 理論上, 收斂速度可以隨意增加。 但是, 在此我們也必須提醒讀者的是, 每一個反正切項的分子 和分母, 其數值成長非常地快速。 例如, 為了要每做一個迴圈計算大約可以增加 4 位圓周率小 數位數, 我們必須有

π

4 = 61 tan⁻¹ 685, 601 69, 049, 993

+ 22 tan⁻¹ 17, 014, 395, 978, 641, 495, 719 2, 082, 431, 079, 356, 889, 478, 358

 , 相關的級數為

π

4 = 61 · 2²¹· N1· D1

10²²

∞

X

n=0

χ(n)2²¹· N₁² 10²²

n

+ 22 · 2⁶¹· N2 · D2

10⁶¹

∞

X

n=0

χ(n)2⁶¹· N₂² 10⁶¹

n

. 最後, 讀者可能會問到, 我們在此所提出的計算圓周率的方法有那些優點? 以級數 (12) 為例, 它的收斂速度很快, 僅需計算五個迴圈就可得到 15 位圓周率小數位數。 而且, 由於這是只有兩 個反正切項的 Euler 型公式, 程式的設計非常簡單。 最重要的是, 除了 χ(n) 之外, 我們沒有用 到除法這個無效率的運算, 這就是為何我們說這個方法有趣的由來。

本文作者已過世的恩師 Dar´ıo Castellanos 曾經在私人往來的信件中告訴本文作者:

“Combinations of arctangent constitute a method which is mostly of historical interest. Large-scale precision calculations of π are done today with recourse to the theory of elliptic modular functions. Nonetheless, the classical methods provide series such as (12) which comprise interesting approaches for the calculation of this most unique of all numbers, and will no doubt, continue to fascinate generations to come.”

這段話的後半段敘述是真的: 這個最著名的超越數至今仍吸引無數的數學家和數學愛好者 成為“digit hunter”。 但是, 在 2002 年的十二月六日, 日本東京大學金田康正教授和早稻田大 學後保範教授等人組成的團隊, 在考慮了速度, 效率, 記憶體和數學運算等因素後, 決定不再使 用 Gauss-Legendre-Brent-Salamin 法和 Borwein-Borwein 法, 而改用反正切型圓周率公 式, 利用一臺 Hitachi SR8000/MPP 型的超級電腦, 將這個最著名的數值計算至不可思議的 1,241,100,000,000位小數。 使用的公式為

π

4 = 44 tan⁻¹ 1 57

+ 7 tan⁻¹ 1 239

− 12 tan⁻¹ 1 682

+ 24 tan⁻¹ 1 12943

 和

π

4 = 12 tan⁻¹ 1 49

+ 32 tan⁻¹ 1 57

− 5 tan⁻¹ 1 239

+ 12 tan⁻¹ 1 110443

 . 第一條公式是由 St¨ormer 在 1896 年發現, 第二條公式是由高野喜久雄先生在 1982 年發現。

本文作者亦獨立地於 1994 年發現此公式 [6], 不但與有榮焉, 而且亦可告慰恩師在天之靈。

參考文獻

1. P. Beckmann, A History of π, 2nd edition, St. Martin’s, New York, p.132-133, 1971.

2. M. Wetherfield, The enhancement of Machin’s formula by Todd’s process, Math. Gaz., Vol. 80, p.333-344, 1996.

3. D. H. Lehmer, On arccotangent relations for π, Amer. Math. Monthly, Vol. 45, p.657- 664, 1938.

4. 黃見利, 如何衍生 Machin 型公式, 數學傳播季刊, 29 卷 1 期 (113), 民 94 年 3 月。

5. 黃見利, 勒貝格單調收斂定理一個有趣的應用: 証明尤拉的反正切公式, 數學傳播季刊, 31 卷 2 期 (122), 民 96 年 7 月。

6. Hwang Chien-lih, More Machin-type identities, Math. Gaz. Vol. 81, p.120-121, 1997.

———本文作者現就讀於國立臺灣大學數學研究所碩士班三年級———