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?
[BCC16-B]
單維彰 (2012/10/10) ---
[Prev] [Next] [Up]