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[]