www.digitalmars.com

D Programming Language 1.0


Last update Sun Dec 30 20:34:42 2012

D Complex Types and C++ std::complex

How do D's complex numbers compare with C++'s std::complex class?

Syntactical Aesthetics

In C++, the complex types are:
complex<float>
complex<double>
complex<long double>
C++ has no distinct imaginary type. D has 3 complex types and 3 imaginary types:
cfloat
cdouble
creal
ifloat
idouble
ireal
A C++ complex number can interact with an arithmetic literal, but since there is no imaginary type, imaginary numbers can only be created with the constructor syntax:
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
In D, an imaginary numeric literal has the 'i' suffix. The corresponding code would be the more natural:
creal a = 5;		// a = 5 + 0i
ireal b = 7i;		// b = 7i
c = a + b + 7i;		// c = 5 + 14i
For more involved expressions involving constants:
c = (6 + 2i - 1 + 3i) / 3i;
In C++, this would be:
c = (complex<double>(6,2) + complex<double>(-1,3)) / complex<double>(0,3);
or if an imaginary class were added to C++ it might be:
c = (6 + imaginary<double>(2) - 1 + imaginary<double>(3)) / imaginary<double>(3);
In other words, an imaginary number nn can be represented with just nni rather than writing a constructor call complex<long double>(0,nn).

Efficiency

The lack of an imaginary type in C++ means that operations on imaginary numbers wind up with a lot of extra computations done on the 0 real part. For example, adding two imaginary numbers in D is one add:
ireal a, b, c;
c = a + b;
In C++, it is two adds, as the real parts get added too:
c.re = a.re + b.re;
c.im = a.im + b.im;
Multiply is worse, as 4 multiplies and two adds are done instead of one multiply:
c.re = a.re * b.re - a.im * b.im;
c.im = a.im * b.re + a.re * b.im;
Divide is the worst - D has one divide, whereas C++ implements complex division with typically one comparison, 3 divides, 3 multiplies and 3 additions:
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;
}
To avoid these efficiency concerns in C++, one could simulate an imaginary number using a double. For example, given the D:
cdouble c;
idouble im;
c *= im;
it could be written in C++ as:
complex<double> c;
double im;
c = complex<double>(-c.imag() * im, c.real() * im);
but then the advantages of complex being a library type integrated in with the arithmetic operators have been lost.

Semantics

Worst of all, the lack of an imaginary type can cause the wrong answer to be inadvertently produced. To quote Prof. Kahan:

"A streamline goes astray when the complex functions SQRT and LOG are implemented, as is necessary in Fortran and in libraries currently distributed with C/C++ compilers, in a way that disregards the sign of 0.0 in IEEE 754 arithmetic and consequently violates identities like SQRT( CONJ( Z ) ) = CONJ( SQRT( Z ) ) and LOG( CONJ( Z ) ) = CONJ( LOG( Z ) ) whenever the COMPLEX variable Z takes negative real values. Such anomalies are unavoidable if Complex Arithmetic operates on pairs (x, y) instead of notional sums x + i*y of real and imaginary variables. The language of pairs is incorrect for Complex Arithmetic; it needs the Imaginary type."

The semantic problems are: Appendix G of the C99 standard has recommendations for dealing with this problem. However, those recommendations are not part of the C++98 standard, and so cannot be portably relied upon.

References

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.



Forums | Comments |  D  | Search | Downloads | Home