Dの複素数型 vs C++のstd::complex

Dの複素数は、C++のstd::complexクラスと比較するとどうなのでしょう？

構文的な美しさ

C++ では、複素数型は：

complex<float>
complex<double>
complex<long double>

です。C++ では虚数型というものは区別されていません。D には 3 つの複素数型と 3 つの虚数型があります：

cfloat
cdouble
creal
ifloat
idouble
ireal

C++の複素数は数値リテラルと相互利用できますが、虚数型を持たないため、虚数はコンストラクタを呼び出す構文でのみ作れます：

complex<long double> a = 5;		// a = 5 + 0i
complex<long double> b(0,7);		// b = 0 + 7i
c = a + b + complex<long double>(0,7);	// c = 5 + 14i

Dでは、虚数レテラルは接尾辞 'i' で表します。対応するコードはより自然なものに見えるでしょう：

creal a = 5;		// a = 5 + 0i
ireal b = 7i;		// b = 7i
c = a + b + 7i;		// c = 5 + 14i

定数の混ざったもっと複雑な式でも：

c = (6 + 2i - 1 + 3i) / 3i;

C++では、これはこんな風になります：

c = (complex<double>(6,2) + complex<double>(-1,3)) / complex<double>(0,3);

あるいはC++に虚数クラスが追加されたら、こうなるでしょう：

c = (6 + imaginary<double>(2) - 1 + imaginary<double>(3)) / imaginary<double>(3);

要するに、Dでは虚数nnは、complex<long double>(0,nn) のようなコンストラクタ呼び出しではなく、 nni と書くだけで表現できます。

効率性

C++に虚数型が欠如していることは、純虚数に対する操作の際に、 0の実部に関する余計な計算が発生することにつながっています。例えば、 D で二つの純虚数を足すには１回の加算で済みます：

ireal a, b, c;
c = a + b;

C++ では、実部も加算されるため２回の加算になります：

c.re = a.re + b.re;
c.im = a.im + b.im;

乗算はもっと深刻で、１回の乗算で済むところが４回の乗算と２回の加算になってしまいます。

c.re = a.re * b.re - a.im * b.im;
c.im = a.im * b.re + a.re * b.im;

除算に至っては最悪で - Dでは１回の除算ですが、C++ での典型的な実装では１回の比較、３回の除算、３回の乗算に３回の加算が行われます：

if (fabs(b.re) < fabs(b.im))
{
    r = b.re / b.im;
    den = b.im + r * b.re;
    c.re = (a.re * r + a.im) / den;
    c.im = (a.im * r - a.re) / den;
}
else
{
    r = b.im / b.re;
    den = b.re + r * b.im;
    c.re = (a.re + r * a.im) / den;
    c.im = (a.im - r * a.re) / den;
}

この非効率性をC++で避けるには、虚部の計算を double を使ってシミュレートする方法があります。例えば、次のDのコード：

cdouble c;
idouble im;
c *= im;

は、C++で次のように実装されます：

complex<double> c;
double im;
c = complex<double>(-c.imag() * im, c.real() * im);

しかしこうすると、std::complex が算術演算子を備えたライブラリである、という利点が失われてしまいます。

意味論

最大の欠点は、虚数型が無いことで、気付かないうちに誤った計算結果が得られる可能性があることです。 Kahan教授の著述から引用すると：

"Fortranや、現在C/C++コンパイラと共に配布されているライブラリに必要な SQRTやLOGといった複素関数の実装では、この流線はおかしな方向に流れることになります。 IEEE 754 で規定された 0.0 の符号を無視した結果、負の実数値を持つ複素数Zに関して SQRT( CONJ( Z ) ) = CONJ( SQRT( Z ) ) や LOG( CONJ( Z ) ) = CONJ( LOG( Z ) ) といった等式が成立しなくなっているのです。この異常現象は、複素数演算が実数と虚数変数の和 x + i*y ではなくペア(x, y)に対して行われる限り、避けられません。ペアを使った言語は複素数演算に関して不正確です。虚数型が必要なのです。"

この意味論的問題は、まとめると：

式 (1 - infinity*i) * i を考えると、この結果は (infinity + i) となるべきです。しかし、もし第二オペランドが単なる i でなく (0 + i) だったとすると、結果は (infinity + NaN*i) となり、誤ったNaNが作られます。
虚数型を区別することで、分枝切断を持つ計算に必要な 0 の符号は保持されます。

C99標準の Appendix G にはこの問題を扱うために推奨される方法が書かれていますが、この付記はC++98標準の一部ではないため、環境を問わず信頼することはできません。

参考文献

How Java's Floating-Point Hurts Everyone Everywhere Prof. W. Kahan and Joseph D. Darcy

The Numerical Analyst as Computer Science Curmudgeon by Prof. W. Kahan

"Branch Cuts for Complex Elementary Functions, or Much Ado About Nothing's Sign Bit" by W. Kahan, ch.
7 in The State of the Art in Numerical Analysis (1987) ed. by M. Powell and A. Iserles for Oxford U.P.

参考書籍

Dの複素数型 vs C++のstd::complex

構文的な美しさ

効率性

意味論

参考文献