複素数の対数領域計算
ComplexPlotというグラフ表示アプリを作った。この手のアプリはDesmosやGeoGebraで良いだろと思うかもしれないが、複素数のプロットがあまり便利でないので、それに特化した表示器を作った次第だ。(リンク)
これを開発するにあたって、問題になったのが計算可能なレンジだ。元々電子回路の周波数vs利得とかの計算を便利にしたくて作ったやつなので、両対数表示とかをしたくなる。しかし対数スケールでズームアウトしていくと、まあまあ簡単に値が計算可能な幅を超えてしまう。
JavaScriptの数値型は倍精度浮動小数で、最大で10^308くらいである。まあ大抵のユースケースで収まることは収まるかもしれないが、5秒くらいズームアウトし続けたらすぐに「端」が見えてしまうのは心許ないといえば心許ない。
巨大な数を安定的に扱う手法といったら対数領域計算である。英語では Logarithmic Number System(LNS) などとも呼ぶらしい。今回初めて実装してみたので、その知見を書き留めておく。
対数領域計算とは
名前の通り、1とか100とかの数値をそのまま扱うのではなく、対数をとった上で扱う計算方法だ。これによって扱える数値の幅は飛躍的に上がる。logを取るので、単純計算で$\pm e^{10^{\pm 308}}$くらいの幅が扱える計算だ。
対数なので、乗除算は加減算に変換される。そのため乗除算との相性は特に良い。しかし加減算との相性は絶望的に悪い。
複素対数
今回は複素数を扱うことが主目的なので、ただの対数ではなく複素対数を扱うことになった。
複素数の対数を取った場合、一般的な定義では、絶対値の対数が実部になり、偏角が虚部に変換される。
有名なオイラーの公式を踏まえればピンとくるのではないだろうか。
$$e^{ix}=\cos(x)+i\sin(x)$$
指数関数に虚数を与えると原点を中心に回転する。対数関数はその逆である。
純虚数ではない複素数を与えた場合は、その組み合わせである。
$$e^{a+bi}=e^a e^{bi}=e^a(\cos(b)+i\sin(b))$$
$$\log\left[e^a(\cos(b)+i\sin(b))\right]=\log(e^{a+bi})=a+bi$$
$$\log z = \log|z| + i \arg z$$
具体的な対数領域計算の中身
複素数の表示
結局コンピュータで計算するときは浮動小数などで表せる普通の実数の計算に帰着しなけらばならない。
複素数は「絶対値の対数」と「偏角」で表す。この方が数学的に筋が良い理由は上に書いた通りだ。「実部の対数」と「虚部の対数」などとすると数学的に面倒な扱いになるのでやらない。
実数の表示
複素数ではなく実数の場合は「絶対値の対数」と「符号」の組で表すことになる。単なる対数では負数が表せないからだ。
後述するが、符号が±の2値に限定されることにより、負の数の計算はより正確になる。
乗算
対数なので加算に変換される。対数領域計算が最も威力を発揮する部分だ。
$$\log(e^z \times e^w)$$ $$ = z+w$$ $$= (\mathrm{Re}[z]+\mathrm{Re}[w])+i(\mathrm{Im}[z]+\mathrm{Im}[w])$$
除算
これも対数だと減算に変換される。
$$\log(e^z / e^w)$$ $$ = z-w$$ $$= (\mathrm{Re}[z]-\mathrm{Re}[w])+i(\mathrm{Im}[z]-\mathrm{Im}[w])$$
符号反転
これは偏角を$+\pi$すればよい。
$$\log(-e^z) $$ $$= \log(-1) + \log(e^z)$$
$$= \log(e^{i\pi}) + \log(e^z)$$
$$= i\pi + z$$
$$= \mathrm{Re}[z] + i(\mathrm{Im}[z]+\pi)$$
虚部(偏角)は$\pm 2\pi n$の違いで同じ数になるので、余りを取って $[0, 2\pi)$ または $[-\pi, +\pi)$ の範囲に収めるとよい。
なお、$\pi$の値を浮動小数で正確に表すことはできないので、ここにはどうしても誤差が出る。もしかすると、虚部を$\pi$の倍数で表すといった手法も考えられるのかもしれない(未検証)。
複素数でなくただの実数であればそもそも±しかないので、符号反転による誤差は出なくなる。 そのため実数であれば実数のための計算にしておくことが吉である。自作のComplexPlotでも、グラフへのプロット時には一部実数を前提とした計算を実装している。
加算
対数領域計算の最も辛いところがここになる。どうしても面倒な数式になる。
とりあえず式を見て欲しい。
$$\log(e^z+e^w)=\log(e^{\mathrm{Re}[z]+i\mathrm{Im}[z]}+e^{\mathrm{Re}[w]+i\mathrm{Im}[w]})$$ $$=\log(\cos(\mathrm{Im}[z])e^{\mathrm{Re}[z]}+i\sin(\mathrm{Im}[z])e^{\mathrm{Re}[z]}+\cos(\mathrm{Im}[w])e^{\mathrm{Re}[w]}+i\sin(\mathrm{Im}[w])e^{\mathrm{Re}[w]})$$ $$=\log(\sqrt{(\cos(\mathrm{Im}[z])e^{\mathrm{Re}[z]}+\cos(\mathrm{Im}[w])e^{\mathrm{Re}[w]})^2+(\sin(\mathrm{Im}[z])e^{\mathrm{Re}[z]}+\sin(\mathrm{Im}[w])e^{\mathrm{Re}[w]})^2})+i\tan^{-1}(\frac{\sin(\mathrm{Im}[z])e^{\mathrm{Re}[z]}+\sin(\mathrm{Im}[w])e^{\mathrm{Re}[w]}}{\cos(\mathrm{Im}[z])e^{\mathrm{Re}[z]}+\cos(\mathrm{Im}[w])e^{\mathrm{Re}[w]}})$$ $$=\mathrm{Re}[z]+\frac{1}{2}\log((\cos(\mathrm{Im}[z])+\cos(\mathrm{Im}[w])e^{\mathrm{Re}[w]-\mathrm{Re}[z]})^2+(\sin(\mathrm{Im}[z])+\sin(\mathrm{Im}[w])e^{\mathrm{Re}[w]-\mathrm{Re}[z]})^2)+i\tan^{-1}(\frac{\sin(\mathrm{Im}[z])+\sin(\mathrm{Im}[w])e^{\mathrm{Re}[w]-\mathrm{Re}[z]}}{\cos(\mathrm{Im}[z])+\cos(\mathrm{Im}[w])e^{\mathrm{Re}[w]-\mathrm{Re}[z]}})$$
どうもゴチャゴチャしているが、とにかく指数関数で一旦対数領域から元の数に戻したりすることを余儀なくされていることが分かると思う。
巨大な数の回避
注目してほしいのが $e^{\mathrm{Re}[w]-\mathrm{Re}[z]}$ の部分だ。対数領域計算の加算においては、どうしても元の数に戻すために指数関数の使用を免れない。しかし、巨大な数を扱うために対数領域で計算しているのに元の数に戻してしまっては何の意味もない。
そこで出てくるのが $e^{\mathrm{Re}[w]-\mathrm{Re}[z]}$ の式である。ここで $\mathrm{Re}[w]-\mathrm{Re}[z]\leq 1$ であれば、指数関数を使っても微小な数しか出さずに済む。これによって対数領域計算のメリットを保ったまま加算を行うのである。
大きさ比較をして入れ替えれば、$\mathrm{Re}[w]-\mathrm{Re}[z]\leq 1$ の前提条件は常に保つことができる。これが安定な加算の計算方法である。
さらなる計算
最後の式のlogの中身は、さらに以下のように計算できる。
$$\frac{1}{2}\log((\cos(\mathrm{Im}[z])+\cos(\mathrm{Im}[w])e^{\mathrm{Re}[w]-\mathrm{Re}[z]})^2+(\sin(\mathrm{Im}[z])+\sin(\mathrm{Im}[w])e^{\mathrm{Re}[w]-\mathrm{Re}[z]})^2)$$ $$=\frac{1}{2}\log(1+2(\cos(\mathrm{Im}[z])\cos(\mathrm{Im}[w])+\sin(\mathrm{Im}[z])\sin(\mathrm{Im}[w]))e^{\mathrm{Re}[w]-\mathrm{Re}[z]}+e^{2(\mathrm{Re}[w]-\mathrm{Re}[z])})$$ $$=\frac{1}{2}\log(1+2\cos(\mathrm{Im}[z]-\mathrm{Im}[w])e^{\mathrm{Re}[w]-\mathrm{Re}[z]}+e^{2(\mathrm{Re}[w]-\mathrm{Re}[z])})$$
だいぶ簡潔な式になった。しかし実は元の式でも、以下のように途中計算を挟むことでそれなり簡潔な式にはなっている。
$$r=e^{\mathrm{Re}[w]-\mathrm{Re}[z]}$$ $$p=\cos(\mathrm{Im}[z])+\cos(\mathrm{Im}[w])r$$ $$q=\sin(\mathrm{Im}[z])+\sin(\mathrm{Im}[w])r$$ $$\log(e^z+e^w)=\mathrm{Re}[z]+\frac{1}{2}\log(p^2+q^2)+i\tan^{-1}(\frac{q}{p})$$
これが計算を進めた後の式だとこうなる。
$$r=e^{\mathrm{Re}[w]-\mathrm{Re}[z]}$$ $$p=\cos(\mathrm{Im}[z])+\cos(\mathrm{Im}[w])r$$ $$q=\sin(\mathrm{Im}[z])+\sin(\mathrm{Im}[w])r$$ $$\log(e^z+e^w)=\mathrm{Re}[z]+\frac{1}{2}\log(1+2\cos(\mathrm{Im}[z]-\mathrm{Im}[w])r+r^2)+i\tan^{-1}(\frac{q}{p})$$
後の式では $\cos(\mathrm{Im}[z]-\mathrm{Im}[w])$ という三角関数の式が新たに1つ増えていることを考えると、個人的にはどちらが良いかは微妙なところではないかと考えている。
大きさ順の計算？
加算の方法は計算を安定させるため、簡単に言えば「大きい方の数をベースとして小さい数を足す」形になっている。両者の差が大きければ大きいほど、小さい側の数は誤差に埋もれやすいことが予想される。
そこで理想的には、多くの項がある足し算は項を絶対値の大きさ順にソートし、小さい数から順に計算するべきであると考えられる。ただこれは単なる2項計算の工夫に留まらず、式の処理順序といった部分にまで絡む処理なので、自分は一旦諦めてしまった。やる気が出たら検証などしてみたいと思う。
減算
これは加算と符号反転の処理が出来れば、以下のように組み合わせてやるだけでよい。
$$A-B=A+(-B)$$
べき乗
加算ほどではないが、まあまあ面倒な計算になる。
$$\log((e^z)^{(e^w)})$$ $$=(e^w)\log((e^z))$$ $$=e^w \times z$$ $$=e^{\mathrm{Re}[w]+i\mathrm{Im}[w]} \times (\mathrm{Re}[z]+i\mathrm{Im}[z])$$ $$=e^{\mathrm{Re}[w]}(\cos(\mathrm{Im}[w])+i\sin(\mathrm{Im}[w])) \times (\mathrm{Re}[z]+i\mathrm{Im}[z])$$ $$=e^{\mathrm{Re}[w]}(\cos(\mathrm{Im}[w])\mathrm{Re}[z]-\sin(\mathrm{Im}[w])\mathrm{Im}[z]) +ie^{\mathrm{Re}[w]}(\cos(\mathrm{Im}[w])\mathrm{Im}[z]+\sin(\mathrm{Im}[w])\mathrm{Re}[z])$$
$\cos(\mathrm{Im}[w]), \sin(\mathrm{Im}[w]), e^{\mathrm{Re}[w]}$はそれぞれ2回出てくるので、メモ化すると良いかもしれない。
その他注意点
対数領域である以上、ゼロの扱いはやや気を付ける必要があるように思われる。
実際のところ、CでもJavaScriptでもlog(0)は-Infinityに評価され、またexp(-Infinity)は0に評価される。これのおかげで殆どの場合は問題を生じないようだ。他の言語のライブラリでも同様の挙動をするかは分からない(少なくともPythonのmathモジュールではlog(0)はmath domain errorを起こす)ので個別に対応する必要があると思うが、このように0と-∞の対応が取れていればそこまで神経質になる必要はないと思われる。
まとめ
まともに計算できるのか不安だったが、全ての内部処理を対数領域で計算しているComplexPlotが動いているのを見れば分かる通り、問題なく計算出来ることが分かった。対数領域なので絶対値が大きい計算ほど誤差は激しくなるものと思われるが、それは浮動小数も同じなので程度問題になる。