大抵のプログラミング言語には指数関数 exp や三角関数 sin, cos, tan などの初等関数が用意されている。これらの関数はよく使うので、標準に用意されているというのは合理的だろう。しかし、初等関数の組み合わせで書けるような関数がわざわざ専用の関数として用意されている場合がある。

指数部の不必要なオーバーフロー・アンダーフローを回避する

浮動小数点数は固定幅の整数や固定小数点数よりは広い範囲の数を表現できるが、倍精度でも指数部はせいぜい11ビット（うち1ビットが符号に相当するので、絶対値として表せるのは \(2^{10}\) 程度まで）なので、 \(2^{1024}\) や \(10^{309}\) というような数を表現できない（inf になる）。

これはまあ仕方ないのだが、計算の途中にこういう数が出現してしまうと、本来の計算結果自体が範囲に収まる場合でも浮動小数点数で計算した結果が意図しない物になる可能性がある。

hypot

hypot は直角三角形の斜辺(hypotenuse)の長さを計算する関数である。数学的にはこうだ：
\[\mathrm{hypot}(x,y)=\sqrt{x^2+y^2}\]
この式に \(x=y=2^{600}\) を代入して、倍精度浮動小数点数として計算してみよう。
\begin{align*}
\sqrt{x^2+y^2}&=\sqrt{(2^{600})^2+(2^{600})^2} \\
&=\sqrt{\mathtt{inf}+\mathtt{inf}} \\
&=\sqrt{\mathtt{inf}} \\
&=\mathtt{inf}
\end{align*}
途中に出てくる \((2^{600})^2=2^{1200}\) の指数部が符号付き11ビットに収まらないので、結果が inf になってしまった。しかし、数学的には \(\sqrt{x^2+y^2}=\sqrt{2}\cdot2^{600}\) であり、この値は倍精度浮動小数点数で表せる範囲に収まっているはずである。

そこで、\(\sqrt{x^2+y^2}\) の計算前に \(x\) と \(y\) を定数倍して、指数部をオーバーフローしない範囲に収めておく。具体的には、 \(x=ax’, y=ay’\), \(a=2^{600}, x’=1, y’=1\) とおくと、
\begin{align*}
\sqrt{x^2+y^2}&=\sqrt{(ax’)^2+(ay’)^2} \\
&=a\sqrt{x’^2+y’^2} \\
&=2^{600}\sqrt{1^2+1^2} \\
&=1.4142135623730951\cdot 2^{600}
\end{align*}
となり、浮動小数点数の範囲で正しい答えが得られる。

複素数の逆数

複素数の逆数は、数学的には次のように計算できる：
\[\frac{1}{x+iy}=\frac{x-iy}{x^2+y^2}\]
しかし、\(x\) や \(y\) がとても大きい（または、絶対値がとても小さい）場合、\(\frac{x-iy}{x^2+y^2}\) の計算結果自体は浮動小数点数の範囲で表せるにもかかわらず、途中の \(x^2\) や \(y^2\) がオーバーフロー（アンダーフロー）する可能性がある。

これも、分母がオーバーフロー（アンダーフロー）しないように、分母分子に定数をかけてやればよい。

ガンマ関数

C99 には特殊関数のガンマ関数が tgamma 関数として用意されている。一方、ガンマ関数の対数を取ったものも lgamma 関数として用意されている。
\[\mathrm{lgamma}(x)=\log(\Gamma(x))\]

一見すると、 lgamma は必要なくて、 tgamma と初等関数の log があれば足りるように思うかもしれない。しかし、御存知の通り、ガンマ関数（階乗）というのは急速に大きくなるので、割と簡単に浮動小数点数で表せる範囲（絶対値が \(2^{1024}\) 未満）を超えてしまう。具体的には、 \(171!\) が \(2^{1026}\) 以上 \(2^{1027}\) 未満くらいの大きさとなる。

そんな場合でも、ガンマ関数の対数であれば浮動小数点数の範囲に収まるので、 \(\mathrm{tgamma}(x)\) が inf となっても \(\mathrm{lgamma}(x)\) が意味を持つ可能性がある。

【2019年3月31日 追記】ガンマ関数は負の値をとることもある。上記の説明は説明不足で、lgamma 関数は、実際にはガンマ関数の「絶対値の」自然対数を返す。処理系によっては、絶対値をとる前のガンマ関数の符号を何らかの方法で返せるようになっている場合がある。

ベータ関数

ベータ関数は、ガンマ関数を使うと次のように書ける。
\[B(x,y) = \frac{\Gamma(x)\Gamma(y)}{\Gamma(x+y)}\]
C言語でこの式の通りに

double beta(double x, double y) {
    return tgamma(x) * tgamma(y) / tgamma(x + y);
}
// printf("%g\n", beta(200, 100)); -> nan

と実装してしまうと、 tgamma がオーバーフローした場合に結果が nan となってしまう。

そこで、対数を取って計算してやる。こうすると、不要なオーバーフローを回避できる。

double beta(double x, double y) {
    return exp(lgamma(x) + lgamma(y) - lgamma(x + y));
}
// printf("%g\n", beta(200, 100)); -> 3.60729e-84

【2019年3月31日 追記】上記の例はガンマ関数が負の値をとった場合の考慮が足りない。

誤差回避

前にちょろっと書いた。

ECMAScriptでのいくつかの数学関数の実装

expm1, log1p

\begin{align*}
\mathrm{expm1}(x)&=\exp(x)-1=x+\frac{x^2}{2!}+\frac{x^3}{3!}+\cdots \\
\mathrm{log1p}(x)&=\log(1+x)
\end{align*}
\(x\approx 0\) のときに \(\exp(x)-1\) を \(\exp(x)\) を計算してから1を引くというやり方で計算すると桁落ちとか情報落ちとかの誤差が発生する。

sinh, cosh, asinh, atanh

同様。双曲線関数や逆双曲線関数は指数関数・対数関数や平方根で書けるが、素朴に実装すると 0 の付近での挙動がよろしくない。

おまけ

【2019年3月31日 追記】以下の内容は最新のNode.jsには当てはまらない。Node.js v10.15.0で試したところ、 Math.asinh(1e-100) も Math.atanh(1e-100) も 1e-100 を返す。

以前の記事に書いたように、 ECMAScript 2015 (ECMAScript 6) では Math.asinh や Math.atanh などの関数が追加された。しかし、現在もっとも広く使われていると考えられる JavaScript の実装である V8 では、誤差を回避する努力が行われていないように見える。

$ node --version
v6.2.2
$ node
> Math.asinh(1e-100)
0
> Math.atanh(1e-100)
0
> Math.asin(1e-100)
1e-100
> Math.atan(1e-100)
1e-100

Math.asin や Math.atan はおそらく C とかの標準ライブラリに投げているので小さい値についても良い結果を返すが、 Math.asinh や Math.atanh はおそらく自前で実装しているので、不正確な答えが帰ってくるのだろう。ちなみに、Math.expm1, Math.log1p, Math.hypot, Math.sinh, Math.tanh などはこの記事に書いたような注意点を考慮して実装されているようだった。（2年前に調べた時は、これらも「雑な」実装だった気がする）

もちろん、 V8 (Chrome) 以外の主要な JavaScript 実装 (Firefox, Safari) ではこれらの関数はまともに実装されているようである。

TODO: V8 でのこのへんの実装を調べる。某 Visual なんとかが最近まで C99 の標準ライブラリ（特に数学関数）を実装してこなかったのが関係有るのかどうか。