std::complex and clang++

[element90]

As well as switching from Gtkmm to Qt for Saturn and Titan development I also changed compilers: clang++ (3.0.2) instead of g++ 4.7.

As part of testing the new version of Saturn and Titan I load seed files produced by previous versions of Saturn and there were some unexpected results i.e. the seed image was not correctly reproduced. I tracked the problem down to the result given by std::norm, the fractal in question had an unnecessarily large bailout value and iteration value reached inf - nan i, the bailout condition used the value from std::norm(inf - nan i) which clang++ returns as nan (not a number). Switching back to g++ the seed file images were correctly reproduced because g++ returns inf (infinity) for std::norm(inf - nan i).

The problem also affects std::abs(inf - nan i).

I managed to work around this problem but there were other fractals that weren't correctly rendered and after a lot a poking around in the guts of the GNU implementation of std::complex I discovered that clang++ returns -nan - nan i for any complex number that has been divided by a zero complex value, g++ on the hand will return inf - nan i, -inf - nan i, -nan + inf i, -nan - inf i &c. and only returns -nan - nan i if both the complex quotient and the complex divisor are zero. The problem is with the clang++ implementation of its built in complex division used by the specialisation code of std::complex<long complex>.

I've been releasing Windows and Linux versions of my programs but I haven't tested the Windows version to anywhere near the extent that I've tested the Linux version so I need to verify that this problem doesn't exist using the Microsoft C++ compiler. For version 4.0.0 of Saturn and Titan I intend releasing a version for OS X, which has clang++ 3.2 and an old version of g++ that doesn't support c++11 features, Saturn on OS X exhibits the same problem as the clang++ built Linux version and doesn't build using g++.

I've also found that the non-specialised code in the GNU std::complex template class contains errors, it is fortunate that I found those as the multi-precision complex class was converted from the std::complex template class as mpfr::mpreal was incompatible with the Windows implementation of std::complex. I suspect the errors haven't been fixed because there are specialisations for float, double and long double so the generic versions aren't used.

The errors in the std::complex generic code are in the implementation of complex number to a complex power and complex division, the complex power will always return 0 if the value to be raised to the specified complex power is zero though zero to a negative power (imaginary part zero) should be inf - nan i and zero to power zero should be nan + nan i, complex division does not take into account zero and will always return -nan - nan i (comments say "This is a grammar school implementation").

[element90]

There's more on this topic on my blog see http://element90.wordpress.com/2013/07/02/an-unexpected-difference/ and http://element90.wordpress.com/2013/07/05/more-on-an-unexpected-difference/.

[hsmyers]

If I remember correctly (always in doubt), IEEE defines the standard with regards to inf and nan and for that matter the rest of the math co-processor (time was). Since all IEEE docs cost a non-small fortune, the fallback is any edition of the appropriate Intel or AMD processor manuals. Then minimally I'd do the following:

• rework the return values
• crib the rest from GCC
• regress
• rinse
• repeat

Given GCC as a source you should at least get matching errors which can be fixed incrementally. This would also give you a cross platform solution as well.

--hsm

[element90]

I've sorted out all the cases involving inf and nan by using two new classes ldcomplex (long double) and mpcomplex (mpfr::mpreal) both adapted from the GNU implementation of std::complex<T>. The only thing that remains to be tracked down is the loss of precision in the power function where the power is long double or ldcomplex,




The first image is produced by Saturn version 4.0.0 (beta) using ldcomplex and the second by Saturn 3.0.1 using std::complex<long double>. Formula Zcpac i.e. z = z^alpha + c, alpha = 4.0 + 0.0i, c is the location in the complex plane and the initial value of z = 0.

Version 3.0.x is implemented using C++/Gtkmm-2.4 and Version 4.0.0 (beta) is implemented using C++/Qt 4.8 on Linux, C++/Qt 5.1 on Windows. The OS X version is implemented using the same code but can't be finished as my iMac has a GPU fault and won't boot into OS X, it'll only boot into command line Linux as X windows does work due to the GPU fault.

Location:

Image centre: 0.629930414528152718527373 + 1.09107091125381725320439i
Image width: 7e-16

Note: there as similar problems with nan and inf when using Microsoft's C++ and std::complex<long double>, additionally Microsoft have defined long double as being the same as double. So using ldcomplex and mpcomplex also fix the problems with the Windows version although it does require the use of mpfr::mpreal and mpcomplex earlier than the Linux version.

[hsmyers]

Good work! Once you get it resolved, it would make a very good blog article. I can then post a link to Hacker News so you can avoid a self post.

--hsm

[element90]

The loss of precision when using power with a long double or ld::complex (previously referred to as ldcomplex) has been tracked down, the arg function used an unqualified atan2 which calls the double version and not the long double version thus losing precision.

One of the pitfalls of C++ is in how it handles numeric functions, to get the right versions the correct namespace MUST be used either via a "using" statement or by qualifying the call (e.g. atan2 would become std::atan2). The use of std namespace can easily be omitted and with g++ and clang++ with the default warnings option you will be none the wiser as the compilers will silently compile the code and use the double versions in all cases (except abs). The the treatment of the abs function is even worse, if the namespace is not used it first casts your float, double or long double into an integer calls the abs function for integers and then casts it back to the required type, in C the abs functions for floats, doubles and long doubles are fabsf, fabs and fabsl.

One of the problems with the std namespace is that often it can be omitted and everything works as expected, and if you use floating values and don't use long double or std::abs nothing amiss will happen. Note: on Windows only the std::abs function is likely to cause problems as long double is the same as double.

[hsmyers]

Post here and today's blog both excellent!   

--hsm