浮動小数点数による複素数の演算に関する注意点


コンピューターで複素数を表す時は、通常は実部と虚部の組をそれぞれ浮動小数点数として持つ。

複素数の演算は、数学的には\begin{align*}
(a+bi) \pm (c+di) &= (a\pm c) + (b\pm d)i \\(a+bi)(c+di) &= (ac-bd) + (ad+bc)i \\
\frac{a+bi}{c+di} &= \frac{(ac+bd)+(-ad+bc)i}{c^2+d^2} \\
\left\lvert a+bi\right\rvert &= \sqrt{a^2+b^2}
\end{align*}

となるので、これに対応する浮動小数点演算を使えば良い。…という単純な話とはならない。

この記事では主にC言語を使って話をするが、浮動小数点数の組として複素数を実装する限りこの記事に書く問題とは無縁ではいられない。むしろ、これらは複素数の演算を実装する時の問題なので、組み込みの複素数型がない環境でこそこの記事の内容が関係してくると言うこともできる。

不必要なオーバーフロー・アンダーフロー (undue overflow and underflow)

除算・逆数と絶対値(実数の言葉で言えば hypot 関数)が不必要なオーバーフロー・アンダーフローを引き起こす可能性およびその対処法については、以前の記事で述べた。よって、ここでは省略する。

浮動小数点数の関数とオーバーフロー

無限大と NaN

掛け算に関しては、不必要なオーバーフローとかそういうのはないので、数式通りに実装すれば良い。…というわけでもない。

複素数の成分のどちらか、または両方が無限大の場合を考えよう。

例:a = inf, b = inf, c = 1, d = 0 の場合:
実部:a * c – b * d = inf * 1 – inf * 0 = inf – nan = nan
虚部:a * d + b * c = inf * 0 + inf * 1 = nan + inf = nan

この例では、無限大 (inf + inf i) と有限の値 (1 + 0 i) の積を計算したところ、実部と虚部の両方が NaN になってしまった。

(無限大)×(無限大)、(無限大)×(有限の値) は NaN ではなく何らかの無限大になるべきだと考えると、この挙動には問題がある。なので、素朴に計算した結果が NaN + NaN i となった場合は、引数が (無限大)×(無限大) または (無限大)×(有限の値) のパターンでないかチェックして、必要なら結果が無限大になるようにうまく計らう必要がある。(C言語の規格の annex には実装例が載っている)

ちなみに、実部か虚部の少なくとも一方が無限大であれば、たとえもう一方が NaN であっても、その複素数は無限大であるとみなす。

複素数の無限大をこういう風に扱うと約束した場合、絶対値の実装時にも注意が必要である。絶対値を(不必要なオーバーフロー・アンダーフローの話は置いておいて)素朴に sqrt(x*x + y*y) で計算すると、片方が無限大であってももう片方が NaN の時に計算結果が NaN となってしまう。ここは適宜チェックして無限大を返すようにする必要がある。

不必要なオーバーフローだとか無限大の取り扱いなんて知るか!俺は数学的な式の通りに計算してくれれば満足なんだ!

不必要なオーバーフローとか無限大の取り扱いとかそういうのはどうでもいいので、とにかく数学の式通りに、変なチェックなしで高速に計算してくれればいい、という人もいるかもしれない。

C言語では、そういう人のために

#pragma STDC CX_LIMITED_RANGE ON
#pragma STDC CX_LIMITED_RANGE OFF
#pragma STDC CX_LIMITED_RANGE DEFAULT

というプラグマが用意されている。

…のだが、2大 C コンパイラーであるところの GCC と Clang はこれを実装していない。ただし、 GCC では -fcx-limited-range コンパイルオプションにより、翻訳単位ごとにこのプラグマの効果を指定できる。Clang は知らん。

複素数の構築時の注意

実部が \(a\), 虚部が \(b\) な複素数の値を作りたいとき、数学なら \(a+bi\) と書く。C言語では #include <complex.h> すれば虚数単位が I という名前のマクロとして定義されるので、 a + b * I と書けば良い。…という単純な話ではない。

C の処理系が _Imaginary 型に対応していない場合、 I_Complex_I, つまり実部が 0 で虚部が 1 の複素数として定義される。b が有限の浮動小数点数ならば問題ないのだが、 b が無限大や NaN の場合は、 b * I = b * (0.0 + 1.0i) = (b * 0.0) + (b * 1.0)i = NaN + bi というように、実部に意図しない NaN が出現してしまう。

C の処理系が _Imaginary 型に対応していて I_Imaginary_I として定義されている場合は、この問題は起こらない。ちなみに、 _Imaginary 型を実装している C 処理系としては、筆者の知る限り Digital Mars C/C++ Compiler と Sun CC があるようである(実際に使って試したわけではない)。

この問題を回避するために、 C11 では CMPLX / CMPLXF / CMPLXL マクロが導入された。これらのマクロを使えば、実部と虚部を別々に指定できる。

ちなみに、 Linux で clang を使う場合や、 OS X で GCC を使う場合は、<complex.h> でこれらのマクロが定義されない場合がある。その場合は、各コンパイラーの拡張機能を使って、次のような定義をどこかに書いておけば良い:

#if !defined(CMPLX)
#if defined(__clang__)
// Clang では _Complex 型を構造体っぽく初期化できる: http://clang.llvm.org/docs/LanguageExtensions.html#initializer-lists-for-complex-numbers-in-c
#define CMPLX(x,y) ((double _Complex){(x), (y)})
#elif defined(__GNUC__)
// GCC では __builtin_complex で複素数を作れる: https://gcc.gnu.org/onlinedocs/gcc-6.1.0/gcc/Other-Builtins.html#index-g_t_005f_005fbuiltin_005fcomplex-4209
#define CMPLX(x,y) __builtin_complex((double)(x), (double)(y))
#endif
#endif

言いたいこと

C 言語の規格(の最終ドラフト)を読もう!

おまけ:cproj 関数

この記事の趣旨とは特に関係ないが、C言語での複素数の扱いについてこの記事でいろいろ書いているので、ついでに触れておく。

C言語では cproj という関数が定められている。この関数でググって出てくる説明文を読むと、「複素数のリーマン球面への射影を計算する」とある。しかし、 cproj 関数の返り値の型も直交座標による複素数を表す double _Complex であり、どこがどうリーマン球面なのか意味がわからない。

先に触れたように、浮動小数点数を使った複素数の演算において、実部と虚部の少なくとも一方が無限大であればその複素数は無限大であるとみなされる。しかし、リーマン球面の無限遠点 ∞ に対応する浮動小数点数の複素数値がたくさんあるというのは、あまり嬉しくない。

そこで、たくさんある無限大を、実軸上の正の無限大 +inf + 0 i に「正規化」することが考えられる。これを行うのが cproj 関数である。細かいことを言うと、虚部の 0 の符号は、元の複素数の虚部の符号が反映される。

浮動小数点数による複素数の演算に関する注意点」への1件のフィードバック

  1. ピンバック: プログラミング言語における複素数の実装 格付けチェック | 雑記帳

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です