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.

Mathematicaの糖衣構文

Mathematicaの関数へのアクセスは糖衣構文によって色々な人に受け入れやすい事になっています．
多分，同じ数学概念でも棲んでいる世界によって解釈が異なったりするからだと思います．ただ，私には甘すぎては苦手です．

可読性って概念はMathematicaにはないのかね…
他人のプログラム読む時，絶望を覚えるんですが・・・

例えば以下の例は全てSqrt[3]を返します．

In[1]:= f[x_, y_: 2] := x^(1/y)

In[2]:= f[3]
f[3, 2]
3~f~2
3 // f
3 // f[#, 2] &
f @ 3
f @@ {3}
f @@ {3, 2}
f[#, 2] &  @ 3
Apply[f, {3}]
Apply[f, {3, 2}]
Apply[f, {##}] &  @@ {3, 2}


Out[2]= Sqrt[3]

Out[3]= Sqrt[3]

Out[4]= Sqrt[3]

Out[5]= Sqrt[3]

Out[6]= Sqrt[3]

Out[7]= Sqrt[3]

Out[8]= Sqrt[3]

Out[9]= Sqrt[3]

Out[10]= Sqrt[3]

Out[11]= Sqrt[3]

Out[12]= Sqrt[3]

Out[13]= Sqrt[3]

最近Mathematicaに関する事ばかりだな・・・

Mathematica Fourier変換について

バグを疑い始めたので，プログラムを純化させてみる．

In[1]:= IdentityTr[x_] := InverseFourier[Fourier[x]]

In[2]:= a = {0, 0, 0.33333333333333333, 0}
b = {0, 0, 0.333333333333333333, 0}
Precision[a[[3]]]
Precision[b[[3]]]

Out[2]= {0, 0, 0.333333, 0}

Out[3]= {0, 0, 0.33333333333333333, 0}

Out[4]= MachinePrecision

Out[5]= 17.5229

In[8]:= NestList[IdentityTr, a, 20]
NestList[IdentityTr, b, 20]

Out[8]= {{0, 0, 0.333333, 0}, {0., 0., 0.333333, 0.}, {0., 0., 
  0.333333, 0.}, {0., 0., 0.333333, 0.}, {0., 0., 0.333333, 0.}, {0., 
  0., 0.333333, 0.}, {0., 0., 0.333333, 0.}, {0., 0., 0.333333, 
  0.}, {0., 0., 0.333333, 0.}, {0., 0., 0.333333, 0.}, {0., 0., 
  0.333333, 0.}, {0., 0., 0.333333, 0.}, {0., 0., 0.333333, 0.}, {0., 
  0., 0.333333, 0.}, {0., 0., 0.333333, 0.}, {0., 0., 0.333333, 
  0.}, {0., 0., 0.333333, 0.}, {0., 0., 0.333333, 0.}, {0., 0., 
  0.333333, 0.}, {0., 0., 0.333333, 0.}, {0., 0., 0.333333, 0.}}

Out[9]= {{0, 0, 0.33333333333333333, 0}, {0.*10^-18, 
  0.*10^-18 + 0.*10^-18 I, 0.33333333333333333, 
  0.*10^-18 + 0.*10^-18 I}, {0.*10^-17, 0.*10^-18 + 0.*10^-18 I, 
  0.3333333333333333, 0.*10^-18 + 0.*10^-18 I}, {0.*10^-17, 
  0.*10^-17 + 0.*10^-17 I, 0.333333333333333, 
  0.*10^-17 + 0.*10^-17 I}, {0.*10^-16, 0.*10^-16 + 0.*10^-16 I, 
  0.333333333333333, 0.*10^-16 + 0.*10^-16 I}, {0.*10^-15, 
  0.*10^-15 + 0.*10^-15 I, 0.33333333333333, 
  0.*10^-15 + 0.*10^-15 I}, {0.*10^-14, 0.*10^-15 + 0.*10^-15 I, 
  0.3333333333333, 0.*10^-15 + 0.*10^-15 I}, {0.*10^-14, 
  0.*10^-14 + 0.*10^-14 I, 0.333333333333, 
  0.*10^-14 + 0.*10^-14 I}, {0.*10^-13, 0.*10^-13 + 0.*10^-13 I, 
  0.33333333333, 0.*10^-13 + 0.*10^-13 I}, {0.*10^-12, 
  0.*10^-12 + 0.*10^-12 I, 0.33333333333, 
  0.*10^-12 + 0.*10^-12 I}, {0.*10^-11, 0.*10^-11 + 0.*10^-11 I, 
  0.3333333333, 0.*10^-11 + 0.*10^-11 I}, {0.*10^-10, 
  0.*10^-11 + 0.*10^-11 I, 0.333333333, 
  0.*10^-11 + 0.*10^-11 I}, {0.*10^-10, 0.*10^-10 + 0.*10^-10 I, 
  0.33333333, 0.*10^-10 + 0.*10^-10 I}, {0.*10^-9, 
  0.*10^-9 + 0.*10^-9 I, 0.33333333, 
  0.*10^-9 + 0.*10^-9 I}, {0.*10^-8, 0.*10^-8 + 0.*10^-8 I, 0.3333333,
   0.*10^-8 + 0.*10^-8 I}, {0.*10^-7, 0.*10^-8 + 0.*10^-8 I, 0.333333,
   0.*10^-8 + 0.*10^-8 I}, {0.*10^-7, 0.*10^-7 + 0.*10^-7 I, 0.33333, 
  0.*10^-7 + 0.*10^-7 I}, {0.*10^-6, 0.*10^-6 + 0.*10^-6 I, 0.3333, 
  0.*10^-6 + 0.*10^-6 I}, {0.*10^-5, 0.*10^-5 + 0.*10^-5 I, 0.3333, 
  0.*10^-5 + 0.*10^-5 I}, {0.*10^-4, 0.*10^-4 + 0.*10^-4 I, 0.333, 
  0.*10^-4 + 0.*10^-4 I}, {0.*10^-3, 0.*10^-4 + 0.*10^-4 I, 0.33, 
  0.*10^-4 + 0.*10^-4 I}}

Mathematicaでの任意制度計算

Mathematicaが手に入ったので，いろいろテストしてみました．
その結果， 任意精度計算を行うにあたって精度の保持はユーザーが監視しないと行けない事が分りました．
また，監視を怠ると正しい結果が得られない事も分りました．

簡単な例として離散Fourier変換(FFTか？)と逆変換を繰り返す演算をさせます．

Wolfram Mathematica Document Foureirの説明にある通り，
Inputのlistが厳密な数値である場合，
N関数を適用し精度はMachinePrecisionに落とされます．
そのためFoureir関数は厳密精度計算を行ってくれません．

しかしながらユーザーがInputのlistにN関数を適用し精度をユーザー側でコントロールしたい状況が発生しますが，安易な方法では正しい結果が得られません．

その結果が In[3]です．

In[1]:= func1[y_, iter_, dps_: 10] := Module[{i, x},
  x = y;
  For[i = 0, i < iter, i++;
   If[ Mod[i, 2] != 0, x = SetPrecision[Fourier[x], dps], 
    x = SetPrecision[InverseFourier[x], dps]];
   ];
  Return[x]]
func0[y_, iter_] := Module[{i, x},
  x = y;
  For[i = 0, i < iter, i++;
   If[ Mod[i, 2] != 0, x = Fourier[x], x = InverseFourier[x]];
   ];
  Return[x]]


In[3]:= a = {0, 0, 1/3, 0}
res = func0[a, 100]
Precision[res[[1]]]
$MachinePrecision

b = N[{0, 0, 1/3, 0}, 10]
res = func0[b, 100]
Precision[res[[1]]]


Out[3]= {0, 0, 1/3, 0}

Out[4]= {0., 0., 0.333333, 0.}

Out[5]= MachinePrecision

Out[6]= 15.9546

Out[7]= {0, 0, 0.3333333333, 0}

Out[8]= {0.*10^5, 0.*10^4 + 0.*10^4 I, 0.*10^5, 0.*10^4 + 0.*10^4 I}

Out[9]= 0.

listに厳密な数値を代入した場合，正しい結果がMachinPrecisionで帰ってきます．
一方listにN関数を作用させて作った初期条件の場合，全ての結果のPrecisionが0になっています．

これでは使い物になりません．１回のFourier関数の演算でどうやら10^2-10^10程度の桁落ちが発生し最終的にゼロが帰ってくる様です．

上記問題を回避する為にはSetPrecision関数を用いて，強制的に精度を維持する事です
(func1関数を参照)．
ただしこの方法が正当な回避方法かどうか私には分りませんが，少なくとも正しい答えが望んだ精度で帰ってくる様です(In[10]以降)．

In[10]:= a = {0, 0, 1/3, 0}
res = func1[a, 100, 10]
Precision[res[[3]]]
b = N[{0, 0, 1/3, 0}, 10]
res = func1[b, 100, 10]
Precision[res[[3]]]

Out[10]= {0, 0, 1/3, 0}

Out[11]= {0, 0.*10^-21 + 0.*10^-21 I, 0.3333333333, 
 0.*10^-21 + 0.*10^-21 I}

Out[12]= 10.

Out[13]= {0, 0, 0.3333333333, 0}

Out[14]= {0, 0.*10^-21 + 0.*10^-21 I, 0.3333333333, 
 0.*10^-21 + 0.*10^-21 I}

Out[15]= 10.

In[34]:= Exit[]