計算 pi 的小程式

多年以前，曾經對如何用電腦計算 π 的小數展開值感到興趣。

例如，要如何算出 π，精確到小數點後 1000 位的值呢？課堂上，似乎只知道可以運用數值方法（例如牛頓法），或者 π 的無窮級數展開式（Wikipedia 關於 π 的詞條，有列出數十個這樣的展開式）來計算。

Wiki 雖然有列出數學式，但很可惜沒有相對應的程式碼，因此總覺得有些遺憾。前些時日，在這裡找到一份 Tcl 寫的程式碼，可以計算出 π 的前 2400 個小數位數。

將程式碼轉為 PHP script 如下：

<?php
   $digits = 2400;
   $size = (int) ($digits * 14 / 4);
   
   $e = 0;
   for ($b=0; $b<=$size; $b++) {
      $f[$b] = 2000;
   }
   
   for ($c=$size; $c>0; $c-=14) {
      $d = 0;
      for ($b=$c; $b>0; $b--) {
         $g = 2 * $b - 1;
         $d = ($d * $b) + ($f[$b] * 10000);
         $f[$b] = ($d % $g);
         $d = (int)($d / $g);
      }
      printf("%04d", $e + $d / 10000);
      $e = ($d % 10000);
   }
?>

頗難想像，計算 π 到 2400 位小數，竟然只需要如此簡短的程式碼。根據網頁上的說明，這份程式碼是依據 Stanley Rabinowitz 與 Stan Wagon 在 1995 年的論文寫成的，背後用到的數學式子是 π 的連分數展開式：

π = 2 + 1/3( 2 + 2/5( 2 + 3/7( ... ( 2 + k/(2k+1)(...)))))

我花了一些時間，嘗試了解程式是如何運用這個連分數展開式，來精確地計算出 π 的小數位數。可惜，直到現在，這份嘗試還是失敗的。或許，原始的論文有說明如何用程式實作；或許，我就是缺乏這樣的能力，能從數學展開式裡，直觀地看出如何漂亮地實作。

然而，有程式的原始碼與背後的數學式，可以跑出結果、程式碼也很簡短，但卻看不懂程式背後的演算法，說來是不是也有些諷刺呢？這是不是也暗示，程式如果失去了註解說明、失去了解讀的指引，那麼想要直接從原始碼來了解演算法，或許也還是一項艱困的事情？