Maximaによる任意精度計算

Mathematicaはいつ使えなくなるか(ライセンス的な関係で)分らないので，Maximaを任意精度数値計算に使えないか試してみる．


精度(precision)のコントロールは

fpprec:50

演算による多倍長浮動小数点の有理数への変換を制御するため

bftorat: true

にする必要がある．デフォルトではfalse．
falseの場合多倍長浮動小数点は小さい有理数で近似される．

また，出力する桁のコントロールは

fpprintprec:50;

で可能．
maxima 5.30.0 だとload("fft")が出来なかったので，古いバージョンを使う．
参考：http://maxima.sourceforge.jp/maxima_5.html

$ maxima
Maxima 5.21post http://maxima.sourceforge.net
using Lisp SBCL 1.0.39
Distributed under the GNU Public License. See the file COPYING.
Dedicated to the memory of William Schelter.
The function bug_report() provides bug reporting information.

(%i1) fpprec;
(%o1) 16
(%i2) fpprec:50;
(%o2) 50
(%i5) bftorat:true;
(%o5) true
(%i6) x:bfloat([0,1/3,0,0]);
(%o6) [0.0b0,3.3333333333333333333333333333333333333333333333333b−1,0.0b0,0.0b0]
(%i7) load("fft");
(%o7) "/Applications/Maxima.app/share/maxima/5.21post/share/numeric/fft.lisp"
(%i10) res:inverse_fft(fft(x));
(%o10) [0.0,.3333333333333333,0.0,0.0]
(%i11) bfloatp(res[1]);
(%o11) false

FFTは任意精度計算やってくれないっぽい．一方で行列mat=((i,pi),(pi,i))の対角化は任意精度でやってくれるっぽい．上のmatの固有値は(±pi, i)である．

(%i12) mat:bfloat(matrix([%i,%pi],[%pi,%i]));
(%o12) matrix([%i,3.1415926535897932384626433832795028841971693993751b0],[3.1415926535897932384626433832795028841971693993751b0,%i])
(%i13) load("eigen");
(%o13) "/Applications/Maxima.app/share/maxima/5.21post/share/matrix/eigen.mac"
(%i14) res:eigenvalues(mat);
`rat' replaced -9.8696044010893586188344909998761511353136994072408B0 by -41494057121218467653886768/4204226981644629723547549 = -9.8696044010893586188344909998761511353136994072408B0
(%o14) [[(4204226981644629723547549*%i−4*sqrt(10903152157933135742303986703335184117163956870727))/4204226981644629723547549,(4204226981644629723547549*%i+4*sqrt(10903152157933135742303986703335184117163956870727))/4204226981644629723547549],[1,1]]
(%i24) bfloat(realpart(res[1][1])); imagpart(res[1][1]);
(%o24) −3.1415926535897932384626433832795028841971693993751b0
(%o25) 1
(%i22) -bfloat(%pi);
(%o22) −3.1415926535897932384626433832795028841971693993751b0
(%i26) bfloat(realpart(res[1][2])); imagpart(res[1][2]);
(%o26) +3.1415926535897932384626433832795028841971693993751b0
(%o27) 1
---------- Mathmatica res ----------
In[2]:= N[Pi, 50]
Out[2]= 3.1415926535897932384626433832795028841971693993751

maxima対角化は時間が掛かりすぎるから現実的ではない・・・か

(%i1) f[i,j]:=random(1.0);
(%o1) f[i,j]:=random(1.0)
(%i6) showtime:true;
Evaluation took 0.0000 seconds (0.0000 elapsed) using 0 bytes.
(%o6) true
(%i12) mat:genmatrix(f,4,4);
Evaluation took 0.0000 seconds (0.0000 elapsed) using 0 bytes.
(%o12) matrix([.9138095996128959,.1951846177977887,.4823905248516196,.3554400571592862],[.4155627227214829,.6125738892538246,.8418126553359213,.5055894446438118],[0.0363933158999219,.6252975183418072,.9937314578694818,.04158924615444981],[.5931521815169765,.6546577661145074,.4707106740336644,.05192355366451751])
(%i13) eigenvalues(mat);
Evaluation took 9.2310 seconds (14.0640 elapsed) using 2764.567 MB.