manページ  — MATH

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

解説 関数リスト 注 バグ 関連項目 注

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

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

名称    ページ解説エラー拘束 (ULP)

acossin.3m 逆三角関数3
acosh      asinh.3m逆双曲線関数3
asinsin.3m 逆三角関数3
asinh      asinh.3m逆双曲線関数3
atansin.3m 逆三角関数1
atanh      asinh.3m逆双曲線関数3
atan2      sin.3m逆三角関数2
cabshypot.3m複素絶対値1
cbrtsqrt.3m立方根1
ceilfloor.3m以上の整数0
copysign   ieee.3m符号ビットのコピー0
cos sin.3m 三角関数1
coshsinh.3m双曲線関数3
erf erf.3m 誤差関数???
erfcerf.3m 補足誤差関数???
exp exp.3m 指数1
expm1      exp.3mexp(x)-11
fabsfloor.3m絶対値0
floor      floor.3m以下の整数0
hypot      hypot.3mユークリッド距離1
ilogb      ieee.3m指数抽出0
j0j0.3m    ベッセル関数???
j1j0.3m    ベッセル関数???
jnj0.3m    ベッセル関数???
lgamma     lgamma.3m対数ガンマ関数(旧 gamma.3m)
log exp.3m 自然対数1
log10      exp.3m底が10の対数3
log1p      exp.3mlog(1+x)1
pow exp.3m 指数 x**y60-500
remainder  ieee.3m余り0
rintfloor.3m最も近い整数に丸める0
scalbn     ieee.3m指数調整0
sin sin.3m 三角関数1
sinhsinh.3m双曲線関数3
sqrtsqrt.3m平方根1
tan sin.3m 三角関数3
tanhsinh.3m双曲線関数3
y0j0.3m    ベッセル関数???
y1j0.3m    ベッセル関数???
ynj0.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, i80287National Semiconductor  32081
     Motorola 68881Weitek WTL-1032, ... , -1165
     Zilog Z8070Western 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,&#177;0) のみです。 とくに比較 (x > y, x ≥ y など) は 0 の符号から 影響を受けませんが、有限な値 x = y である場合は、 Infinity

amp;= 1/(x-y) !=

!= -1/(y-x) = -Infinity

- になります。

Infinity

は符号付きです。

Infinity

をそれ自身か有限の数値に加算しても 無限は維持されます。 Infinity

の符号は乗算と除算で正しく変換され、 (有限値)/&amp;#177;Infinity = &amp;#177;0

(有限値)/&amp;#177; = &amp;#177;0 (非 0 値))/0 = &amp;#177;Infinity

&amp;#177; になります。しかし、 無限-無限, 無限*0, 無限/無限

∞-∞, ∞*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 除算 ±無限 

オーバフロー±∞
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)	amp;... その他の方法については、この文書では説明しません。

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

i)	関数に提供したデータが値しない例外を発生させない。
ii)	発生したすべての例外をその関数で識別し、サブルーチンでは 識別しない。
iii)	関数の定義が例外処理と関連しているにも関わらず、上記 5 つの 例外処理方法によって呼出し側プログラムを変更し、小さな関数の 内部処理を中断しない。

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

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

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

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

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

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

その他の場合

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

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

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

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

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

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

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

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

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

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

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

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』の記事も役立ちますが、この 記事はのちに変更された標準仕様案に関係しています。