COMPUTER ARITHMETIC パッケージ チュートリアル
コンピュータ演算パッケージ
Mathematica で使われる演算は,可変精度のソフトウェア演算と浮動小数点ハードウェアの製造元が提供したもの(浮動小数点ハードウェアがない場合はコンパイラの設計者が提供したもの)とが混ざったものである.コンピュータの浮動小数点計算一般についての基本事項について知りたい場合,または使用中のマシンでの機械計算を調べたい場合は,
を使うとよい.これを使うとさまざまな基底,精度,丸め規則の演算について調べることができる.
| ComputerNumber[x] | 通常の数 x を,現在使用中の演算におけるコンピュータ数に変換する |
| ComputerNumber[sign,mantissa,exp] | 符号 sign, 仮数 mantissa, 指数 exp を持つコンピュータ数を形成する |
|
| コンピュータ数を構成する完全なデータオブジェクト |
| NaN | 現行の演算で表示不能な数 |
におけるコンピュータ数と非数
コンピュータ数を構成するデータオブジェクトに含まれている情報のほとんどは不要なものである.特に,最初の3つの引数は4番目の引数と全く同じ情報を有している.この冗長さは効率のためでもあり,ユーザが多くのフィールドにアクセスできるようにするためでもある.第5引数はコンピュータ数自体とは何の関係もなく,コンピュータ数に含まれるようになったすべての丸め誤差の累積的な影響を除いた場合にその数の値がどのようなものになるかということを表す. これは Mathematica の高精度数学をもって計算されたもので,一般にその数の正確な値であると考えられる.コンピュータ数をこの数と比較するとコンピュータ数における誤差が分かる.
| Out[2]= |  |
次は

をコンピュータ数で表したものである.
| Out[3]= |  |
| Out[4]= |  |
次は

の完全な構造である.
Out[5]//InputForm= |
| |  |
| Out[6]= |  |
しかし,現在有効な演算で入力が意味を成さなければ
NaNが返される.ここでの問題は仮数が3桁しかない点である.
| Out[7]= |  |
使用する演算の型の変更
デフォルトの演算は基底10で桁数4,丸め規則は
RoundToEvenである.

から

までの数のみが許容される.混合モードの演算は許容されず割算は正しい丸め操作では行われない(これは逆数による掛け算という2つの操作で行われ,それぞれの操作に丸め誤差が含まれる).
| Out[8]= |  |
| Out[9]= |  |
| Out[10]= |  |
演算の設定にはいくつかのオプションがある.精度と基底を変更する他に,使用する丸めのタイプ,許容する数の大きさ,混合モードの演算を許容するか否か,割算を単独の丸め操作で行うか否か等が変更できる.
SetArithmeticのオプション
正確に丸められた割算は関数IdealDivideで実装されることに注目しなければならない.この関数は演算が正確な割算を自動的に使用するかどうかに関わらず使うことができる.
はパッケージに認識される前に
に変換されるので,これを使って正確な割算を得るのは難しい.パーサが
を
に変更してしまう前に,$PreReadを使って
を x~IdealDivide~y に変換させて,これを避ける.自分自身の理由で$PreReadを使いたい場合は,これによって自身の定義が干渉されることがある.
使用可能な丸めのタイプ
混合モードの演算(整数とコンピュータ数の両方を含む演算的操作すべて)は許容されない.
| Out[11]= |  |
これは混合モードの演算となる(そして演算は16進法の6桁に設定される).
| Out[12]= |  |
| Out[13]= |  |
計算機数学には普通の数学と異なる点が多数ある.その大部分は簡単に示すことができる.
| Out[14]= |  |
式はコンピュータ数に変換される前に数値的に計算される.
| Out[15]= |  |
別々の部分が先にコンピュータ数に変換されると,数が破滅的に簡約されることがある.
| Out[16]= |  |
項の最小から最大までの和を求めると答は1つになる.
Out[17]//FullForm= |
| |  |
原則としては最小から最大へと加算していく方がよい.仮数(第2引数)と正しい値(最後の引数)を比較すると誤差が分かる.
Out[18]//FullForm= |
| |  |
次は最大から最小へと加算した方結果がよい場合の例である.
Out[19]//FullForm= |
| |  |
Out[20]//FullForm= |
| |  |
基本的な演算はすべてこのパッケージに実装されている.初歩的な関数を含むように簡単に拡張することもできる.
| Out[21]= |  |
| Out[22]= |  |
演算が正しい場合,正しく丸められた小さな整数の平方根は2乗すると常にもとの整数になるという法則がある.
| Out[23]= |  |
| Out[24]= |  |
次は演算を基底が10で桁数が7,丸め規則は
RoundToEvenで-50から50の指数区間に変更したものである.
| Out[25]= |  |
| Out[26]= |  |
次の数は奇数(9999999)ではなく偶数 (1000000)の仮数に丸められるので端数を切り上げる.
| Out[27]= |  |
これも偶数の仮数(1000001ではなく1000000)に丸められる.
| Out[28]= |  |
逆数の逆数はもとの数ではない.実際,もとの数とは甚だしく異なることがある.
| Out[29]= |  |
次は逆数の乗算である.丸め操作が2回行われている.
| Out[30]= |  |
| Out[31]= |  |
Normalはコンピュータ数を全く同じ値の普通の有理数に変換し直す.
| Out[32]= |  |
| Out[33]= |  |
次は,例題中で非数値の値をプロットしようとして出されたエラーメッセージを表示しないようにする.
各数のコンピュータ数表示におけるエラーをプロットするのは簡単である.
| Out[35]= |  |
| Out[36]= |  |
Mathematica はマシンが提供する機械計算と独自の任意精度の計算の両方を使う.任意精度の計算はどんな機械を使っても同じである(表示可能な最大数の違いとマシンのメモリ量によっての違いはある).このパッケージは,マシンが提供する浮動小数点演算を扱う.
コンピュータ上の数値は離散集合である.数値の間にはギャップがあり,計算をしても結果が表示できないことがしばしばある.結果が2つの表示可能な数値の間に当たる場合は,正確な結果ではなく表示可能な数値のうちの近い方を使うことが最善である.浮動小数点数の形式でコンピュータ上で表現することのできる数の集合は,一般に機械数集合と言われる.このパッケージを使うと,計算をコントロールする特定の基本的なパラメータを変えることができる.
連続する2つの機械数の差は1 ulp(unit in the last place,最後の桁の単位)と呼ばれる.1 ulpの大きさは機械数集合のどの部分で演算を行っているかによって異なる.1と2の間なら1 ulpは$MachineEpsilonに等しい.2と4の間なら2$MachineEpsilonに等しい.どんな関数も2分の1 ulpを超える誤差を含む結果を返さないことが望ましい.なぜなら,真の結果から最も近くの機械数までの距離は常に2分の1 ulpより小さく,最悪の場合でも2つの機械数のちょうど中心にあるからである.四則演算と平方根関数ではこの理想に到達するのは比較的簡単である.初歩的な超越関数は1つの計算で評価することができないので,すべての丸め誤差を合わせて2分の1 ulpより小さく収めなければならないことになり,この理想を達成するのは少々難しくなる.しかし,これらの関数で誤差が1 ulpより大きくならないようなアルゴリズムを設計するのは可能である.
計算機演算パッケージは,使用しているマシンでの機械浮動小数点計算を検証するのに有益な4つの関数を提供する.
機械浮動小数点計算を検証する関数
MicroscopePlotとMicroscopicErrorPlotにおける点
のデフォルトの近傍は,a の左に30ulpから a の右に30ulpの間である.ここで使用されているulpは
で定義されたもの,すなわちUlp[a]である.2のベキ乗付近と0付近ではulpの大きさが変化するので,どうしたらよいか定かでない.2のベキ乗付近では,ulpは2つの値の小さい方になるように選ばれていて,丸めの影響によって結果の機械数が複数回含められるので,横座標にはそれ程問題はない.近傍のデフォルトサイズは近傍指定に3番目の値を含むことで変更できる.
MicroscopePlotとMicroscopicErrorPlotでは,関数の値が無限大になる付近で問題が起こる.関数は大きさの順に急速に変化しているので,知りたい情報を表示するのによいスケールを選ぶのは困難である.そのような場合は,自分の所望する点から何百ulpか離れたところで関数を検証する必要があるかもしれない.
近傍の大きさのコントロール
オプションJoinedはプロットの描画方法をコントロールする.現在検証中の関数は機械数を機械数にマップする,ということはいくら強調しても強調しすぎることはない.我々の目的のためには,間にある実数は存在しないと言ってもよい.ゆえに,MicroscopePlotとMicroscopicErrorPlotが生成するプロットはプロットされた一群の点の集合にすぎない.しかし,その点を直線で結んでみると何が起こっているかがよく分かるようになる.最後に,関数が実数から実数へのマップとして,すなわち実際の引数を最も近くの機械数に丸めてこの結果に関数を適用することで何を行うのかを知りたければ,そのようにしても構わないが,誤差が大きくなることが分かっても,それが関数の問題だと結論付けてはならない.関数は,それが引数として受け取る機械数を扱う以上のことはできないからである.
MicroscopePlotとMicroscopicErrorPlotのオプション
| Out[37]= |  |
| Out[38]= |  |
| Out[39]= |  |
これは,

を最も近くの機械数に丸め,次にこれの平方根を取った際の誤差を合わせたものである.
| Out[41]= |  |
| Out[42]= |  |

に近い別の数では,平方根関数の誤差が1 ulpの半分より小さい.
| Out[43]= |  |
関数

はコンピュータの製作者によって与えられたものではないので,与えられた操作のネストしたシーケンスとして評価する.これはシーケンスにどれほどの誤差が含まれるかは示すが,与えられる情報は余り役に立たない.テスト目的では,誤って条件付けされた関数が前の誤差を拡大することを避けるために,一般に機械数で評価された個々の関数に注目することが望ましい.
| Out[44]= |  |
次は
Log関数を7近くの機械数で評価したものである.
| Out[45]= |  |
次は
Log関数を7近くの10+1+10個の機械数で評価し,各点を直線で繋いだものである.
| Out[46]= |  |
次は実数を最も近くの機械数に丸めることによる効果と,結果の
Logを取ることによる効果を一緒にしたものである.これは
Log関数の有効なテストではない.
| Out[47]= |  |
これは16付近の31個の機械数で正弦関数を評価した際の誤差を示している.16のところでulpがどれほど変化するかに注目のこと.また,丸めの結果が同じ点になることがあるので,点の数が31よりも少なく見えることにも注目のこと.垂直軸上のスケールはulpを表している.
| Out[48]= |  |
次は正弦を31付近の21個の機械数で評価した際の誤差である.各点を直線で繋いである.
| Out[49]= |  |
次は実数を最も近くの機械数に丸め,結果の正弦を取ったものの効果をまとめたものである.これは正弦関数の有効なテストではない.正弦関数を31付近の機械数で評価した際の誤差は,最初に実数を最も近い機械数に丸めてからその機械数の正弦を取った結果生じる誤差を合わせたものよりもはるかに小さい.
| Out[50]= |  |
次は平方根関数が実際にはどれほどひどいかを示している.

付近で誤差がおよそ

ulpほどになっている.
| Out[51]= |  |
この関数は厳密に1の点が特異点になっているので,この点では評価できない.しかし,1近くの点でなら評価することができる.

のラベルは

と出力される.ラベルが6文字に制限されているからである.
| Out[52]= |  |
丸め誤差がしばしば誤解され,誤差とは何の関係もないアルゴリズムが不当に責められるということがよくあるので,1点ここではっきりさせておく必要がある.アルゴリズムは機械数の集合を機械数の集合にマップするように設計されているため,これ以外の数を考慮に入れても意味がないのである.コンピュータに関する限り,他の数というものは存在しないのである.例えば
の平方根を求めたければ次のように行う.
まず

を数値的に近似することから始める.
| Out[53]= |  |
次の式は

に最も近い機械数を求めるものである.
| Out[54]= |  |
次は,平方根のための機械精度のアルゴリズムを使って機械数

の平方根を求める.
| Out[55]= |  |
次は

の誤差をulpで計る.この結果は例題を実行した特定のマシンにおける誤差を与える.ユーザが使用しているマシンでは結果が甚だしく異なることも考えられる.
| Out[56]= |  |
以下では誤差が1 ulpの半分より小さくなっている.しかし,ここでの誤差には平方根関数と
という数を
という機械数に丸めたときの丸め誤差の両方が含まれている.実際のところ,平方根関数の誤差は0.1ulpより小さい.
次は

の平方根の正確な値である.
| Out[57]= |  |
| Out[58]= |  |