C言語をはじめとする多くのプログラミング言語には、exp, log, sqrt などの指数、対数、平方根の関数や、 cos, sin, tan などの三角関数、cosh, sinh, tanh などの双曲線関数が備わっている。これらの数学関数はECMAScript (JavaScript)からは、Mathオブジェクト経由で Math.exp や Math.cos などのようにアクセスできる（が、現状ECMAScriptには双曲線関数はない）。
さて、C99でこれらに加えて、acosh, asinh, atanh などの逆双曲線関数や、expm1 や log1p などの「精度を意識した」関数などが追加された（expm1は\(\exp x-1\), log1p は \(\log(1+x)\) を計算する関数である）。これらの関数をECMAScriptで使うにはどうすればいいか？
ECMAScriptの次期標準であるECMAScript 6では、双曲線関数 Math.(cosh|sinh|tanh) や、C99で追加された関数に対応する Math.(acosh|asinh|atanh) や Math.expm1, Math.log1p が入るらしい（Firefoxなどでは既に実装されている）。しかし、ECMAScript 6に非対応の環境では、これらの関数は自前で実装するしかない。
幸い、双曲線関数は指数関数で書ける。よって、例えば Math.sinh はECMAScriptに既に存在する Math.exp 関数を使って次のように書ける:
[sourcecode lang=”js”]
Math.sinh = Math.sinh || function(x) {
return (Math.exp(x)-Math.exp(-x))/2;
};
[/sourcecode]
Mozillaのサイトにも代替コード(Polyfill)として同じようなコードが載っている。めでたしめでたし。
といいたいところだが、精度を気にする場合はそうはいかない。例えば、\(x\) にとても小さい数、例えば \(x=10^{-20}\) を代入してみるとどうなるだろうか。\(\sinh\) のテーラー展開\[
\sinh x=x+\frac{x^3}{3!}+\cdots
\]を考えると（\(x^3\) は微小なので無視して）値はおおよそ \(x=10^{-20}\) と一致するはずである。この場合、ECMAScriptの「数」はIEEE754の64bit浮動小数点数で、精度は53bitしかないため、\(x\) が十分に小さければ \(x+\frac{x^3}{3!}+\dots\) は完全に \(x\) と一致する。
一方、指数関数のテーラー展開は\[
\exp x=1+x+\frac{x^2}{2!}+\cdots
\]で、今書いたのと同じ理由で Math.exp(x) や Math.exp(-x) は完全に 1 と等しくなる。従って、さっき実装したオレオレ Math.sinh は0を返す。情報落ちだか桁落ちだか忘れたが、なんかそういう現象が起こっている。これは理想的とはいえない。ちなみに、C言語の sinh はこの場合「正しい」 \(10^{-20}\) を返す。
同様に、Math.expm1 を次のように愚直に実装しても、同じ問題が発生する。
[sourcecode lang=”js”]
Math.expm1 = Math.expm1 || function(x) {
return Math.exp(x)-1;
};
[/sourcecode]
では、どうすれば \(x\) が小さい時でも「正しい」答えを出す expm1 関数を作れるか？（双曲線関数よりも指数関数の方が単純なので、以後こちらを議論する）
まあ単純な答えとしては、 \(x\) が小さいときに関数のテイラー展開を使って多項式として計算してやれば良い。さっきのように \(x=10^{-20}\) とかだと1次近似で \(x\) そのものを返してやればよかったが、この近似は \(x\) がどのぐらい小さければ問題ないのか？あるいは、近似の次数をどのぐらい増やしてやれば誤差はどれぐらい小さくなるか？
仮に、\(\exp x\) を \(n\) 次の多項式で近似するとしよう。すると、真の値（無限級数で表される）との差は\[
\left\lvert(\text{差})\right\rvert=\left\lvert\exp x-1-\sum_{k=1}^{n}\frac{x^k}{k!}\right\rvert=\left\lvert\sum_{k=n+1}^{\infty}\frac{x^k}{k!}\right\rvert
\]となる。これを適当に評価してやれば良い。実際にやってみると、例えば次のようになる。\begin{align*}
\left\lvert(\text{差})\right\rvert=\left\lvert\sum_{k=n+1}^{\infty}\frac{x^k}{k!}\right\rvert
&\le\sum_{k=n+1}^{\infty}\frac{\lvert x\rvert^k}{k!} \\
&=\frac{\lvert x\rvert^{n+1}}{(n+1)!}\sum_{k=0}^{\infty}\frac{\lvert x\rvert^k}{(n+2)\cdot\dots\cdot(n+k+1)} \\
&\le\frac{\lvert x\rvert^{n+1}}{(n+1)!}\sum_{k=0}^{\infty}\frac{1}{(n+2)^k} \quad (\left\lvert x\right\rvert <1\text{を仮定}) \\
&\le\frac{\left\lvert x\right\rvert^{n+1}}{(n+1)!}\frac{1}{1-\frac{1}{n+2}} \\
&\le\frac{\left\lvert x\right\rvert^{n+1}}{(n+1)!}\frac{n+2}{n+1}=\frac{(n+2)\left\lvert x\right\rvert^{n+1}}{(n+1)(n+1)!}
\end{align*}
では、どういうときに \(n\) 次の多項式近似の値を使ってよいかというと、近似値 \(\sum_{k=1}^{n}\frac{x^k}{k!}\) の大きさに対してこの \(\left\lvert(\text{差})\right\rvert\) の大きさが（精度が53bitの浮動小数点数の場合は） \(2^{-53}\) よりも小さければ確実に問題ない（十分条件）。つまり\begin{align*}
\left\lvert\frac{(\text{差})}{\sum_{k=1}^{n}\frac{x^k}{k!}}\right\rvert<2^{-53}
\end{align*}だ。この左辺をもうちょっと分かりやすい形で評価すると
\begin{align*}
\left\lvert\frac{(\text{差})}{\sum_{k=1}^{n}\frac{x^k}{k!}}\right\rvert&\le\left\lvert\frac{(\text{差})}{x}\right\rvert \\
&=\left\lvert\frac{1}{x}\frac{(n+2)\left\lvert x\right\rvert^{n+1}}{(n+1)(n+1)!}\right\rvert \\
&=\frac{(n+2)\left\lvert x\right\rvert^n}{(n+1)(n+1)!}
\end{align*}
となる（書いた後に気付いたが、最初の行の評価が \(x<0\) のときにマズい気がする）。つまり、\(n\) 次近似を使っても良い十分条件として\[
\frac{(n+2)\left\lvert x\right\rvert^n}{(n+1)(n+1)!}<2^{-53},
\]すなわち\[
\left\lvert x\right\rvert<\sqrt[n]{\frac{(n+1)(n+1)!}{n+2}2^{-53}}
\]が得られる。
例えば、\(n=3\) として右辺を計算すると \(1.286976\ldots\times 10^{-5}\) が得られるので、さっきのオレオレ Math.expm1 の実装は
[sourcecode lang=”js”]
Math.expm1 = Math.expm1 || function(x) {
if (Math.abs(x) <= 1.28e-5) {
return x*(1+x/2*(1+x/3));
} else {
return Math.exp(x)-1;
}
};
[/sourcecode]
のように改良できる。
まあ、改良したと言っても所詮お遊びで、実際のC言語のライブラリの実装はこんな適当なコードよりもちゃんとやっているので、真面目にやりたい人はC言語のライブラリの実装を移植しよう。
Math.expm1 の他に、Math.log1p, Math.sinh, Math.tanh, Math.asinh, Math.atanh あたりの関数を実装する時にも同じような（小手先の）テクニックが使える。
ところで、MDNに載っている Math.tanh のPolyfillに問題があるのを見つけた。MDNのPolyfillでは x が Infinity かどうかを関数の最初でチェックしているが、x が Infinity でなくても Math.exp(x) が Infinity になる場合を考慮していない。その結果、Math.tanh(x) に 1 を返してほしいのに実際には NaN が返ってくる場合がある。
最後に全く関係ない話だが、プログラミング言語Luaは標準ライブラリに math.sinh などの双曲線関数を計算する関数(中ではC言語のライブラリ関数を呼び出す)を持っている。しかし、どうやら将来的にはこれらの双曲線関数を標準ライブラリから外すようだ（次期バージョン5.3でdeprecated扱い）。理由として作者は「あまり使われていない」「必要なら mathx などの外部ライブラリを使えば良い」としている（この記事に書いた理由により、「math.exp で実装できる」は理由にならない）。ECMAScript と Lua はよく似ていると言われる事が多いが、一方がWebの世界のアセンブラとなるべく数学関数・浮動小数点数関係の関数を追加しようとしているのに対し、もう一方はそれを削減しようとしているのは対照的だと思った。（まあ筆者自身はLuaを活発に触っていたのは5.1時代の事で、5.2が出て以降はあまり触っていないので、ぶっちゃけどうでもいいのだが）
