math − 数学ライブラリ関数の導入

こ の 関 数 は、C の 数学ライブラリ libm を構成します。リンクエディタ は、“−lm” オプションが指定されると、このライブラリを検索します。この 関 数の宣言は、インクルードファイルで入手できます。

以 下のそれぞれの倍精度浮動小数点数関数は、単精度関数と対になっており、 名前が f で終わる関数は単制度浮動小数点 数 を 扱 い ま す。 た と え ば acos(double x) は、浮動小数を扱う acosf(float x) と対になっています。

acos

sin.3m

逆三角関数

3

acosh

asinh.3m

逆双曲線関数

3

asin

sin.3m

逆三角関数

3

asinh

asinh.3m

逆双曲線関数

3

atan

sin.3m

逆三角関数

1

atanh

asinh.3m

逆双曲線関数

3

atan2

sin.3m

逆三角関数

2

cabs

hypot.3m

複素絶対値

1

cbrt

sqrt.3m

立方根

1

ceil

floor.3m

以上の整数

0

copysign

ieee.3m

符号ビットのコピー

0

cos

sin.3m

三角関数

1

cosh

sinh.3m

双曲線関数

3

erf

erf.3m

誤差関数

???

erfc

erf.3m

補足誤差関数

???

exp

exp.3m

指数

1

expm1

exp.3m

exp(x)−1

1

fabs

floor.3m

絶対値

0

floor

floor.3m

以下の整数

0

hypot

hypot.3m

ユークリッド距離

1

ilogb

ieee.3m

指数抽出

0

j0

j0.3m

ベッセル関数

???

j1

j0.3m

ベッセル関数

???

jn

j0.3m

ベッセル関数

???

lgamma

lgamma.3m

対数ガンマ関数(旧 gamma.3m)

log

exp.3m

自然対数

1

log10

exp.3m

底が10の対数

3

log1p

exp.3m

log(1+x)

1

pow

exp.3m

指数 x**y

60−500

remainder

ieee.3m

余り

0

rint

floor.3m

最も近い整数に丸める

0

scalbn

ieee.3m

指数調整

0

sin

sin.3m

三角関数

1

sinh

sinh.3m

双曲線関数

3

sqrt

sqrt.3m

平方根

1

tan

sin.3m

三角関数

3

tanh

sinh.3m

双曲線関数

3

y0

j0.3m

ベッセル関数

???

y1

j0.3m

ベッセル関数

???

yn

j0.3m

ベッセル関数

???

カリフォルニア大学から 1985 年末に配布された 4.3 BSD では、上記の関数に 2 つのバージョンが用意されています。1 つは DEC VAX−11 ファミリ の コ ン ピュー タの倍精度 "D" フォーマットで、もう 1 つは IEEE 754 標準のバイナ リ浮動小数点数演算に準拠した倍精度演算です。 UNIX が誕生したときのも の よ りも、他のプログラムが正確で頑丈になったことから予想できるように、こ の 2 つのバージョンも、非常に似たような挙動を示します。たとえば、これら の プログラムは、上の表の ulps の数以内で正確です。ulp は、最終桁の単位 です。これらのプログラムは、かつての数学ライブラリ libm で発生してい た 異 常状態を修正しました。この異常状態が発生すると、以下のような結果が得 られていました。

sqrt(−1.0) = 0.0 であり log(−1.0)= −1.7e38
cos(1.0e−11) > cos(0.0) > 1.0
x = 2.0, 3.0, 4.0, ..., 9.0 のときに、pow(x,1.0) != x
pow(−1.0,1.0e10) が整数オーバ フ ロー 例 外 を 起 こ し て い た sqrt(1.0e30) と sqrt(1.0e−30) の計算が非常に遅かった

しかしこの 2 つのバージョンには、説明しなければならない違いがあります。 以下に注意してください。

DEC VAX−11 D_floating−point:

これは、オリジナルの数学ライブラリ libm が開発されたフォーマッ ト で あ り、このマニュアルが原則的に依って立つフォーマットです。これは、 PDP−11 および初期の VAX−11 マシンの倍精度フォーマットです。 1983 年 以 降 の VAX−11 は、IEEE 倍精度フォーマットにより近い、オプションの "G" フォー マットを備えています。初期の DEC MiroVAX には D フォーマットが な く、G 倍精度のみが備わっています。

D_floating−point のプロパティ:

ワードサイズ: 64 ビット、8 バイト。基数: バイナリ。
精度: 56 ビット、10 進数で約 17 桁。

x と x’ が連続した正の D_floating−point 浮動小数点数 (1 ulp ずつ異なる) である場合は、
1.3e−17 < 0.5**56 < (x’−x)/x ≤ 0.5**55 < 2.8e−17
となる。

範囲:  オーバフローしきい値           = 2.0**127= 1.7e38

アンダフローしきい値

= 0.5**128= 2.9e-39

注意: この範囲は相対的に狭くなっています。

オーバフローが発生すると、計算は常に中断されます。
アンダフローが発生すると、0 になります。
警告:

アンダフローにより、x != y で x−y = 0 となる可 能 性 が あ ります。同じように、 x > y > 0 であって も、x∗y = 0 もしくは y/x = 0 ということが警告なし に発生してしまいます。

0 は曖昧に表現されます。

ハードウェアでは、2**55 通りの 0 の表現を受理しますが、生 成されるものは明らかな表現だけです。 VAX には -0 がありま せん。

VAX アーキテクチャには無限 ( Infinity ) もありません。
予約オペランド:

ハードウェアは 2**55 通りを認識しますが、生成されるものは 1 つだけです。予約オペランドに対する浮動小数演算 は、MOVF や MOVD でさえ計算を中断するので、あまり使用されません。

例外:

0 除 算とオーバフローするオペレーションは不正なオペレー ションで、計算は中断されます。古いマシンでは予約オペラ ン ドが作成され、計算が中断されます。

丸め:

(PDP−11 では必ずしもそうとは言えませんが) VAXでのすべて の有理数演算 (+, −, ∗, /) は、オーバフ ロー、 ア ン ダ フ ロー、0 除算が発生しない場合、ulp の半分以内に丸められま す。丸め誤差が ulp のちょうど半分である場合、丸めは 0 か ら離れます。

範 囲が狭いことを別にすると、D_floating−point は、1960 年代に設計された ものとしては優れたコンピュータ算術演算のひとつで す。 D_floating−point のプロパティは、4.3 BSD で配布された VAX の基本的な関数に忠実に反映され ています。オーバフローやアンダフローが発生するのは、結果が範囲外にな る か 範囲外に非常に近い場合のみで、その場合はオーバフローやアンダフローを 起こす有理算術オペレーションと同じような動 作 を し ま す。 同 じ よ う に、log(0) と atanh(1) のような式は 1/0 のように動作し、sqrt(−3) と acos(3) は 0/0 のように動作します。これらはすべて予約オペランドを 作 成 し、 計算は中断されます。この状況については、それぞれのマニュアルページ で詳しく説明します。

このレスポンスは厳し過ぎるので、すべての浮動小数演算の例 外 を 適切に処理するように開発された、柔軟で統一されたスキーム で近い将来に置き換えられる予定です。

UNIX 用 4.3 BSD の新しい libm の関数を DEC の VAX/VMS ライブラリと比 較 すると、一部の VMS 関数は多少速く、一部は多少精度が高く、一部はかなり例 外に厳しくなって(pow(0.0,0.0) と atan2(0.0,0.0) など)おり、 ほ と ん ど は、libm よりも多くのメモリを使用します。 VMS コードは、大きいテーブル で補間して速さと精度を達成しています。 libm コードは、将来はす べ て が ROM に収まりそうなほど小型で相当工夫を凝らした式を使用しています。

それより重要なことは、DEC が VMS コードを特許とみなし、許可なく使用する ことを厳しく禁じているのに対し、4.3 BSD の libm コードは public domain を 意図しているということです。つまり、出所を常に明らかにして、ユーザが コードを使用した感想を報告してコードの作者の研究に協力する限り に お い て、libm コー ド は 自 由 に コピーできます。それゆえ、算術処理が VAX D_floating−point に似ているマシンの UNIX ユーザであっても、新しい libm よりも劣るものを使用する必要がないということです。

IEEE STANDARD 754 浮動小数点数算術演算:

こ の 標 準は、他のコンピュータ算術演算方式よりも広く採用されてきていま す。以下のような多くのメーカが、この標準の一部のバージョンに準 拠 し た VLSI チップを製造しています。

     Intel i8087, i80287      National Semiconductor  32081

Motorola 68881

Weitek WTL-1032, ... , -1165

Zilog Z8070

Western Electric (AT&T) WE32106.

こ の 他 にも、ソフトウェアによる実装 (Apple Macintosh は完璧にそれ) か ら、 Hewlett-Packard 9000 シリーズの VLSI による実装や、果ては ECL を駆 使 し、3 Megaflop を達成した ELXSI 6400 まで極めて広い実装の幅を誇りま す。他の企業の中には、丸め、オーバフロー、アンダフローなどの例外の処 理 で こ の標準の方法に従わずに、IEEE 754 のフォーマットだけを採用していま す。 DEC VAX G_floating−point フォーマットは、IEEE 754 Double フォー マッ ト に 非常に似ているので、上に挙げたほとんどの初等関数の IEEE バー ジョンの C プログラムは、MicroVAX で実行するように簡単に変換で き ま す が、わざわざ変換してやろうという人はいないようです。

4.3 BSD libm のコードのうち、IEEE 754 準拠マシン用のものは、 National Semi. 32081 と WTL 1164/65 用になっています。 Intel チップや Zilog チッ プ、または Apple Macintosh か ELXSI 6400 でこのコードを使用する場合は、 これらの企業が提供するよりよいコード (無料になると思われる) か、 上 の コー ド の 作 者 が設計したコードを使用することになります。 atan, cabs, cbrt, erf, erfc, hypot, j0−jn, lgamma, pow, y0−yn を除け ば、 Motorola 68881 は libm のすべての関数をチップに搭載し、しかもより速く正確になっ ています。 Motorola, Apple, i8087, z8070, WE32106 では、有 効 桁 数 64 ビットを使用しています。4.3 BSD の libm コードには、 public domain を意 図しているという長所があります。出所を常に明らかにし、ユーザがコード を 使用した感想を報告して作者の研究に協力すれば、libm コードは自由にコピー できます。算術演算が IEEE 754 に似ているマシンの UNIX ユーザは、新し い libm に劣るものを使用する必要がないということです。

IEEE 754 Double-Precision のプロパティ:

ワードサイズ: 64 ビット、8 バイト、基数: バイナリ
精度: 53 ビット、(10 進数で) 約 16 桁。

x と x’ が連続した正の倍精度浮動小数点数である (1 ulp ず つ異なる) 場合は、
1.1e−16 < 0.5**53 < (x’−x)/x ≤ 0.5**52 < 2.3e−16 と な る。

範囲: オーバフローしきい値         = 2.0**1024= 1.8e308

アンダフローしきい値

= 0.5**1022= 2.2e-308

オー バフローが発生すると、デフォルトで符号付きの無限にな ります。アンダフローは 0.5**1074 = 4.9e-324 の倍数の整 数 のうち最も近いものに段階的に丸められます。

0 は、+0 か −0 のようにあいまいに表現されます。

符号は乗算か除算で正しく変換され、符号付き 0 の加算で維持 されます。しかし x が有限の 場合、x−x は +0 になります。0 の 符号が問題になる演算は、0 除算と copysign(x,±0) のみで す。とくに比較 (x > y, x ≥ y など) は 0 の符号から影響 を 受 け ま せんが、有限な値 x = y である場合は、 Infinity = 1/(x−y) != −1/(y−x) = −Infinity になります。

Infinity は符号付きです。

Infinity をそれ自身か有限の数値に加算しても無限は維持され ま す。 Infinity の符号は乗算と除算で正しく変換され、 (有 限値)/±Infinity = ±0 (非 0 値))/0 = ±Infinity にな り ま す。 し か し、 無 限− 無 限, 無限∗0, 無限/無限は、0/0 や sqrt(−3) と同様に不正な演算であり、NaN が生成されます。

予約オペランド:

2**53−2 通りの予約オペランドがあり、すべ て NaN (Not a Number) と呼ばれます。一部は発信 NaN と呼ばれ、それに対し て浮動小数点演算を実行するとトラップが発生します。値の 欠 落 か未初期化、または配列の存在しない要素をマークするため に使用されます。その他は無発信 NaN と呼ばれます。これは不 正 演算の結果のデフォルト値で、それ以後の算術演算に影響が 波及します。 x != x である場合、x は NaN です。その他すべ て の論理演算 (x > y, x = y, x < y など) は、NaN が関係す る場合は偽になります。
注意: NaN は三分法を妨害します。

NaN が関係する場合、純粋な (不) 等値演算ではな く 順 序付け比較を伴う論理演算は偽になるだけでなく、 不正演算となります。

丸め:

デフォルトでは、すべての代数演 算 (+, −, ∗, /, sqrt) は、ulp の半分以内に丸められます。丸め誤差が ulp のちょう ど半分である場合、丸めた値の最下位有効ビットは 0 になりま す。 通常はこのような丸めが最適で、たとえば x = 1.0, 2.0, 3.0, 4.0, ..., 2.0**52 である場合、商と積の両方が丸められ る に も 関わらず、 (x/3.0)∗3.0 == x, (x/10.0)∗10.0 == x, ... になります。このような結果となるのは、IEEE 754 の よ うな丸めのみです。しかし 1 つの丸めがどの状況でも最適であ るとは限らないため、 IEEE 754 で は、0 へ の 丸 め、 +Infinity への丸め、 -Infinity への丸めをプログラマが選択 できるようになっています。最低でも 1.0e−10 から 1.0e37 く ら いまででは、バイナリと小数の変換に同じ種類の丸めを指定 できます。

例外:

IEEE 754 では、5 種類の浮動小数例外が認識されます。以下で は、重要な順にその 5 種類を挙げます。

例外                   デフォルトの結果

不正演算

NaN, または偽

オーバフロー

±無限

0 除算

±無限

アンダフロー

段階的アンダフロー

不精密

丸めた値

注意: 例外は、適切に処理されればエラーではありません。例 外を例外的なものにしてしまうのは、すべての例を満 た す デ フォ ルトレスポンスがない場合です。逆に、デフォルトレスポ ンスがほとんどの例を満たす場合、それでも満たされ な い 例 を、 例外が発生するたびに演算を中断するという手で正当化す ることはできません。

IEEE 754 では、それぞれの浮動小数例外のためにフラグを用意して い ます。このフラグは例外が発生するたびに上がり、プログラムでリセッ トするまで、上がった状態で残ります。プログラムでは、フラグのテス ト、保存、復元もできます。 IEEE 754 には、デフォルトの結果が満足 の行くものでない例外にプログラムで対処する方法として、以下 の 3 つがあります。

1)

後で例外の原因となる可能性がある条件をテストし、そのテ スト結果に従って処理を変えて例外を避ける。

2)

フラグをテストし、プログラムで最後にフラグをリセット し てから例外が発生したかどうかを確認する。

3)

結果をテストし、1 つの例外のみが発生する値であるかどう かを確認する。

警告: アンダフローが発生したかどうかを確実に発見するには、積 や 商がアンダフローしきい値より 0 に近いかどうかをテストする か、アンダフローフラグをテストするしかありません。 (IEEE 754 で は、和と差がアンダフローを起こすことはありません。 x != y である場合、x−y は完全に正確で、値が小さくても 0 ではあり ま せ ん。 ) 段階的にアンダフローを起こす積と商は、徐々に精度を 失っても 0 にはならないため、0 と比較しても (VAX で実行す る 可 能 性 がある)、ロスは明らかになりません。段階的にアンダフ ローを起こした値をアンダフローしきい値より大きい値と加算する 場合、段階的なアンダフローで失われる桁は切り捨てられるので失 われません。このため、通常の場合、段階的なアンダフローは無視 で き る ことがあります。 0 にフラッシュされるアンダフローで は、同じことが正しいとは限りません。

IEEE 754 に準拠するシステムでは、以下のような方法でも例外に対 処 できます。

4)

異 常終了 (ABORT)。このメカニズムでは、「ON ERROR GO TO」のようなエラー処理ステートメントに関連する方法で処理 す る 事がらとして、例外を前もって分類します。言語が異なれば、 このステートメントの形式は変わりますが、次のような共通し た 特徴があります。

—

例外を起こした演算結果の値を代用し、式の途中から計算を再 開する方法がなくなります。例外の結果は放棄されます。

—

エラー処理ステートメントがないサブプログラムでは、どのプ ロ グラムが呼び出してもそのサブプログラムは例外のために異常 終了し、エラー処理ステートメントが見つかるまでサブプログ ラ ム が次々に呼び出されるか、タスク全体が異常終了してメモリの ダンプが取られます。

5)

停止(STOP)。対話型のデバッグ環境が必要になるこのメカ ニ ズ ムは、プログラマ用であってプログラム向きではありません。 例外はプログラマの間違いの兆候として前もって分類されま す。 例 外が発生すると、例外を起こしたオペレーションの近くで実行 が停止されるので、プログラマは、例外がどのように発生した の か を確認できます。多くの場合は、最初の数個の例外は大きな問 題ではなく、プログラマはそれぞれの例外が発生するたびに、 実 行が停止されなかったかのように実行を再開できます。

6)

... その他の方法については、この文書では説明しません。

例外処理には、範囲 (scope) という重大な問題があります。問題の解決策は分 かっていますが、人手不足のため、4.3 BSD の libm で配布するときまでに 完 全 に実現できませんでした。理想的には、それぞれの初等関数を、次のような 意味で分割せずに自動的に動作させる必要があります。

i)

関数に提供したデータが値しない例外を発生させない。

ii)

発生したすべての例外をその関数で識別し、サブルーチンでは 識 別しない。

iii)

関数の定義が例外処理と関連しているにも関わらず、上記 5 つ の例外処理方法によって呼出し側プログラムを変更し、小さな関数の内 部処理を中断しない。

プ ログラマが、デバッグ済みサブプログラムをユーザが気付かないような小さ なものに簡単にできるようになれば理想的です。しかし、小さな関数の 3 つす べ て の 特徴をシミュレートすることは、多くのテスト、保存、復元が伴うた め、厄介な作業になっています。現在、この不便さを改善する作業が進んで い ます。

この作業が進むまで、libm の関数はそれほど小さくなりませんが、以下の場合 を除き、不適切な例外を発生させません。

オーバフロー/アンダフロー

正しく計算した結果が、範囲内に収まっている場合。

cabs, cbrt, hypot, log10, pow の不正。

エラーの偶発的な取り消しによって正確になった場合。

その他の場合

不正演算は、以下のような場合にのみ発生します。

NaN 以外の結果が誤っていると思われる場合。

オーバフローは、以下のような場合にのみ発生します。

正しい結果が有限であるが、オーバフローしきい値を越えて い る場合。

0 除算は、以下のような場合にのみ発生します。

有限演算で、関数が無限の値を受けた場合。

アンダフローは、以下のような場合にのみ発生します。

正しい結果が 0 ではないが、アンダフローしきい値より小さい 場合。

不正確は、以下のような場合にのみ発生します。

正しい結果の表現に、より広い範囲かより高い精度が必要な 場 合。

信 号が適切である場合、その信号はコードの特定オペレーションから発信され ているので、上記の手法 5) を使用している場合は、サブルーチンをトレー ス し て そ の信号を発している関数を特定する必要があります。すべてのコード は、IEEE 754 のデフォルトを使用します。つまり、すべての 0 除算をト ラッ プしようとすると、トラップしない場合に 0 除算でも正しい結果を出すコード が中断されるということです。

fpgetround(3), fpsetround(3), fpgetprec(3), fpsetprec(3), fpgetmask(3), fpsetmask(3), fpgetsticky(3), fpresetsticky(3) - IEEE 浮動小数インタ フェース

IEEE 754 とエクステンション p854 については、 1984 年 8 月に発 行 さ れ た、IEEE の 雑 誌 『MICRO』 を 参 照してください。 W. J. Cody らが「A Proposed Radix− and Word−length−independent Standard for Floating-point Arithmetic」 という記事の中で説明しています。Apple Macintosh の Pascal, C, BASIC のマニュアルには、IEEE 754 の機能に関する分かりやすい説明が あ ります。 IEEE の雑誌『COMPUTER vol. 14 no. 3』(1981 年 3 月) と 1979 年 10 月の『ACM SIGNUM Newsletter Special Issue』の記事も役立ちますが、 こ の記事はのちに変更された標準仕様案に関係しています。