[ruby-dev:5776] Precision in Mathematica (Re: Bigfloat)

ごとけんです

遅くなってすみません．とりあえず，Mathematica に関する精度の
取り扱いに関する報告です．

さきに僕の結論をいっておくと，我々の目指す BigFloat は以下で
紹介する Real に相当するのがよく，これはオブジェクトごとに
精度(有効桁数)をもち，演算の結果は自動的に調節されるべきです．
んが，この精度の実現を厳密に行うならば実装は複雑です．とは
いうものの妥協案の提示も難しそうです．

このような精度概念を考慮した ruby の Precision と Bigfloat の
全体像はまだ考え終ってないのでのちほど．


0. 型と評価

  はじめに，Mathematica のオブジェクトについて説明します．
  Mathematica には普通に言われる意味の型はありません．
  全てのオブジェクトはS式に似た頭部(プレフィックス)つきの
  形式 f[...] で表され，これを表現と言います．また，
  表現 expr の頭部は Head[expr] で得ることができます．
  この Head[expr] の値がいわゆる型に相当します．
  例えば Head[1] の値は Integer です．頭部は自由に作る
  ことができます．頭部で識別されるものを型というようです．

  Mathematica にはパターンという概念がありパターンを使って
  関数の引数の型を実現します．そのため演算における型変換の
  一般則は明示できません．例えば頭部は _ で指定できます．
  例として任意係数整数ベキ項の微分の定義と動作を掲げます．

    d[c_ x^n_Integer] := c n x^(n-1)
    d[a x^3] ==> 3 a x^2      # マッチしたので評価される
    d[a x^3.0] ==> d[a x^3.0] # マッチしないのでそのまま

  Mathematica は与えられた表現を可能な限り知っている規則
  で変形し続けるシステムで，この変形のことを評価と言います．
  変形できなくなった時点でそれを値として返します．

1. 組み込みの数

  Mathematica には組み込みの4つの数の型があります．
  Real は任意精度であり無限精度ではありません．

   |  Integer 任意の大きさの整数
   | Rational integer/integer で表せる有理数
   |     Real 任意の精度の近似値で表された実数
   |  Complex number + number I で表現される複素数

2. 定数記号

  この他に Simbol として定義された定数が存在します．
  いくつか例を挙げます．ただし虚数単位 I は Complex[0,1] 
  で定義されますが，Degree は Pi/180 で定義されているわけ
  ではありません．

   |          Pi 円周率
   |           E 自然対数の底
   | GoldenRatio 黄金比
   |      Degree 度からラジアンのへ変換係数 Pi/180
   |  EularGamma オイラーのγ
   |     Catalan カタラン数

  方向の有無で2種類の無限大があります．実数の正負の無限大は
  それぞれ，DirectedInfinity[1], DirectedInfinity[-1] で，
  定義され，Infinity, -Infinity と表すこともできます．

   | DirectedInfinity 実数の無限大
   |  ComplexInfinity 複素数の無限大

3. 精度

  Mathematica において，数は有限の精度を持ち，記号は無限の
  精度を持ちます．数の精度は変更することができますが，記号
  の精度は変更することができません．数の精度には2種類あり
  全体の有効桁数を Precision で，小数以下の有効桁数を
  Accuracy で得ることができます．桁数は10進でのものです．
  Precision の値は一般に0以上ですが，Accuracy は負を返す
  こともあります：

    {Precision[10.^67], Accuracy[10.^67]} ==> {16, -51}

  任意精度の取り扱いは以下のように記されています．

   | Mathematica では，任意の精度の数をある桁数の数値は
   | 確定しているもののその他はまったく不確かな量を表わす
   | 数として取り扱います．一般に，任意の精度による数 x は
   | 確定した Precision[x] 桁の数にそれに続くまったく
   | 不確かな数とみなされます．
   |
   | いったん計算を行うと，Mathematica は不確かな要素が
   | 結果に与える影響がある桁を追跡していきます．そして
   | 影響を受けた部分が含まれないように精度を再設定します．
   | この手続きを採用することで Mathematica が与える結果は，
   | 不確かな桁がどのような数であっても，結果のすべての
   | 桁が常に正確であることが保証されます．

  以下の結果は奇妙にみえますが，有効桁の意味を考えれば妥当です．

    delta = N[3^-30, 20]
      ==> 4.856935749618861138 10^-15
    Precision[delta] ==> 20
    Precision[1 + delta] ==> 34

  不確かな桁の追跡はオブジェクトごとに独立に行われるので，次の
  ような結果を得ます．

    Precision[(1 + delta) - delta] ==> 34

  表現からある精度 pr の数を得るには N[expr,pr] を用います．
  第2引数を省略すると機械精度(倍精度浮動小数点数)になります．
  一旦精度が与えられると，その精度で計算されます．与えた
  精度が不十分な場合は，結果で精度を失うこともあります．

    Sin[N[10^30, 20]] ==> 0.
    {Precision[%], Accuracy[%]} ==> {0, 0}

  表現の精度はその中に現れる数の精度の最小値を意味します．
  例えば，1 == 1. という表現の精度は機械精度で，この比較は
  機械精度で行われ，この場合 True が返されます．ただし
  機械精度以下の精度の数が含まれる表現は機械精度とみなされます．

  次の例では計算過程で精度が失われています．

    In[1]:= x = N[1 - 10^-30,40]
    Out[1]= 0.999999999999999999999999999999000000000
    In[2]:= 1+x
    Out[2]= 1.999999999999999999999999999999000000000
    In[3]:= Precision[%]
    Out[3]= 40
    In[4]:= (x^2-1)/(x-1)
    Out[4]= 2.00000000
    In[5]:= Precision[%]
    Out[5]= 10

  一般には可能な限り高い精度で計算します．

    N[1-10^-30,40] + N[1-10^-30,50]
      ==> 1.999999999999999999999999999998000000000
    {Precision[%], Accuracy[%]}
      ==> {40, 40}

    N[1-10^-30,40] + 1
      ==> 1.999999999999999999999999999999000000000
    {Precision[%], Accuracy[%]}
      ==> {40, 40}

    1+Pi ==> 1 + Pi
    {Precision[%], Accuracy[%]}
      ==> {Infinity, Infinity}

    Sqrt[1/7] ==> 1 / Sqrt[1/7]
    {Precision[%], Accuracy[%]}
      ==> {Infinity, Infinity}

5. まとめ
  
  Mathematica の精度の取り扱いは誤差を含まない範囲で，可能な
  限り長い有効桁数を確保するのであって，単に大きな方に合わせる
  といったようなルースなものではありません．機械精度の数だけは
  例外で表現が機械精度に感染するとそれ以降を全て機械精度で評価
  します．これはパフォーマンスのための選択であり納得できます．

-- gotoken