Maxima 教材：高位數的計算

Maxima 是一個「電腦代數系統」(Computer Algebra System)， 重點就是它做數學上真確的計算，而不是數值的估計。 因此，1 / 3; 的計算結果就是「三分之一」而不是 0.33333333。 如果我們想要將數學代數符號轉換成（十進位的）數值，已經知道要用 float 指令或者在計算過程中故意使用含有小數點的數字。 請看 平方根與十進制數值 教材。

但是，上述方法採用硬體的數值計算，也就是「雙精度浮點數」的計算， 所以只有大約 16 位的有效數字。而 Maxima 其實有軟體計算的功能， 使得它能處理很大整數的真確計算，例如

5!;
會得到「5階乘」的真確整數，也就是 5*4*3*2*1 的結果，120。同理
40!;
就是「40階乘」了，也就是 1*2*3*...*40 的大數：
815915283247897734345611269596115894272000000000
這是一般專門做硬體計算的程式（例如 Matlab）辦不到的事。

Maxima 會在整數運算時自動啟用軟體計算的功能， 但是對實數的估計值，就只會自動使用硬體計算。例如

float(1/3);
得到
0.33333333333333
而
sqrt(2.0);
得到
1.414213562373095

我們用 bfloat (big float: 大浮點數) 啟動軟體數值計算， 而有效位數就可以提高了。可是，bfloat 的內建位數還是 16， 因此看起來並沒有成效。例如

bfloat(sqrt(2));
結果是
1.414213562373095b0
結果和普通的硬體運算一樣，只是在最後添了 b0， 表示「乘以 10 的 0 次方」的意思。

Maxima 用一個參數 fpprec 來控制有效數字的位數 (Floating Point Precision)。它的初始值是 16。如果我們將它設定成 32， 則 bfloat 的計算結果就不同了：

fpprec : 32;
bfloat(sqrt(2));
得到的結果是
1.4142135623730950488016887242097b0

如果有必要使用高位數的數值計算，須注意：

  1. 過程中全用真確的代數計算，到了最後才將真確的解轉換成高位數的估計值。
  2. 從一開始就要使用高位數的估計值。
重點就是，一旦在過程中使用了低位數的數值計算，它帶進來的誤差， 是無法用 bfloat 來彌補的。例如說
bfloat(sqrt(2.0));
則因為 sqrt(2.0) 已經使用了硬體計算，產生了 16 位有效數字的數值， 就算再用高位數的 bfloat 去計算，結果還是錯的。 以上指令產生的輸出是
1.4142135623730951454746218587388b0
但是，其實紅色部分的數值都是錯的；正確部分的數值， 就跟 sqrt(2.0) 的結果一樣，大約只有 16 位。

習題

  1. 請算出 264。
  2. Maxima 用 %pi 表示圓周率。 試問圓周率的小數點下三十位以內，出現過幾次 0？幾次 1？幾次 9？
單維彰 (2012/10/10) --- [Prev] [Next] [Up]