0. 緣由

去年某位前輩在網路上發表了一段文章,其中討論了圓周率的計算,當看到
方法3時小弟真是佩服的五體投地,深覺自己才疏學淺:

句末這位前輩留一下句話:「您覺得您要被訓練多久?才會有這種想法啊?!」

當時小弟心理想,哇賽!少說也要一百幾十年吧(話說絕世武功都要練個一百幾十年),
但是事情真的是這樣嗎?

 

1. 考證

後來,作者看了 Dr. Jack W. Crenshaw 在雜誌上發表的一篇文章 - More Strange Integer,
才知道6千年前埃及人就已經知道類似的作法,假如你跟我當初一樣覺得很神奇的話那表示
你跟小弟一樣該回家多唸書,埃及人的作法是:

哇賽,誤差居然只有 0.000008%,神奇吧!而且不只可以這麼做,連也可以如法炮製:

問題是,這些神奇的技巧是怎麼來的?有沒有法則可尋?

 

2. 實踐

假如你用 355/113 這個法子,前輩還是會笑你分子 355 超過 8 位元能表示的範圍,
113 也不是 2 的冪次,用在 8051 這種 8 位元系統上還是不適合,那 201/64 這個完美
的結果是如何找出來的?

Dr. Crenshaw 給了一個很直接的答案,用暴力法!就以為例,取到小數點後第八位:

我們試圖找出近似 x 藉由比值:

用國小就會的交叉相乘,我們可以得到:

接著我們用 C++ 寫一個很簡單的程式,把每個 q 都試試看(只取 2 的冪次),並且 p 超過
255 的都去掉不要,並且把結果排序:

#include <vector>
#include <stdio.h>
#include <math.h>
#include <algorithm>
#define LIMIT 256
struct Ratio{
   int p, q;
   double err;
};
   
bool less_err(const Ratio& lhs, const Ratio& rhs)
{
   return lhs.err < rhs.err;
}
typedef std::vector<Ratio> RatioVector;
int main(void)
{
   double x; 
   int q;
   RatioVector rv;
   
   x = 3.14159265f;
   
   for(q = 1; q < LIMIT; q<<=1)
   {
       Ratio r;
       r.q = q;
       r.p = int(x * r.q + 0.5);
       double tmp = (double)r.p/r.q;
       r.err = ::abs(x - tmp) * 100.0f;
       if(r.p < LIMIT)
           rv.push_back(r);
   }
   
   std::sort(rv.begin(), rv.end(), less_err);
   
   for(size_t i = 0; i < rv.size(); ++i)
   {
       printf("p = %d, q = %d, error = %f%%\n", rv[i].p, rv[i].q, rv[i].err); 
   }
   
   return 0;
}

 

最後我們看看輸出結果:

果然最佳結果是 201/64!現在各位應該知道201/64並不是亂試試出來的,而是用一個小學就知道
的公式叫電腦去幫我們算(也許前輩提到小學就學過的演算法就是指這個)。各位也可以拿來試
試看,應該也會得到類似的結果。

 

3. 結語

前輩在文中提到很多台灣工程師不事先好好分析就急著寫程式,除了這個毛病外,其實台灣工程師
也很少會找國外的技術刊物來看,各位知道我提的那篇文章是什麼時候刊出的嗎?1995年4月!也就是
12年前國外的雜誌就有寫了,而且這些出版社都很好心幫你把文章收集成一張 CD 讓你網路訂購,
現在您覺得看到問題就馬上硬幹還是應該先看看前人怎麼做的呢? ...(^_^)