## More on an Unexpected Difference2 comments

In the post “An Unexpected Difference” I fixed a problem related to std::complex<long double> and the clang++ compiler, at the end of the post I mentioned that there are other problems with std::complex using that compiler. The following image was produced from a seed file using Titan 3.0.0. Successfully reproduced by Titan 3.0.0

I won’t show you image produced by the new version of Saturn because it is a featureless rectangle. I discovered that failure to reproduce this fractal was again due to clang++.

The fractal type is “Combination 3” and it has the following formula:

z = transform(z)
z = (αz
β + γ)δ + (εzζ + η)θ

z0 = 0
α = the location in the complex plane
β = 2
γ = -1
δ = 3
ε = the location in the complex plane
ζ = -7
η = 1
θ = -1
There are no transforms defined so z = transform(z) does nothing.

Using clang++ the first iteration for all locations results in -nan – nan i so calculation stops after one iteration resulting in a featureless rectangle. So why does this happen?

The error occurs in the second expression of the formula, zero is raised to the power -7 which is infinity, it is then multiplied by the location in complex plane and has 1 added to it so it is still infinity, then infinity is raised to the power -1 which is the same as dividing 1 by infinity which sould result in zero, using g++ it does but using clang++ the result is nan + nan i. Any calculation involving nan (or -nan) results in nan. When printing the value of an std::complex or long double number inf, -inf, nan and -nan can result, nan means not a number and for some strange reason can be printed with a sign which is nonsensical as it is NOT A NUMBER.

Regarding infinity, it can be signed and for complex numbers it can be represented as inf in one of the components of the complex number and the other may contain nan, complex numbers can also have both components inf and they may or may not be signed. In all cases if either of the components is inf or -inf the result of using it to divide another number is zeo unless the number itself is zero in which case the result is nan + nan i as zero divided by infinity is undefined.

Here is a comparision of the results produced by g++ and clang++.

g++

(0,0)^(-7,0) = (inf, nan)
(0,0)^-7 = (inf, -nan)
(0,0)^(-1,0) = (inf, nan)
(inf,nan)^(-1,0) = (0,0)
1/0 = inf
(0,1)/(0,0) = (-nan, -inf)
(1,0)/(0,0) = (inf, -nan)

clang++

(0,0)^(-7,0) = (inf, nan)
(0,0)^-7 = (-nan, -nan)
(0,0)^(-1,0) = (inf, nan)
(inf,nan)^(-1,0) = (nan,nan)
1/0 = inf
(0,1)/(0,0) = (-nan, -inf)
(1,0)/(0,0) = (-nan, -nan)

The differences affect calculations that involve infinity which in some circumstances is not a problem because an operation on infinity can result in zero but once both components of a complex number are nan then any operation with that value will result in not a number.

Delving into the guts of GNU complex template library, the code for a complex number to a complex power is:

```  inline __complex__ long double
__complex_pow(const __complex__ long double& __x,
const __complex__ long double& __y)
{ return __builtin_cpowl(__x, __y); }
```

This code is called for both g++ and clang++, it looks a bit strange as it uses C99 types to handle complex numbers using the compiler’s internal routines which means that the built routine for complex power in the clang++ compiler is at fault.

For division of a complex number the following is called:

```      template<class _Tp>
complex&
operator/=(const complex<_Tp>& __z)
{
_ComplexT __t;
__real__ __t = __z.real();
__imag__ __t = __z.imag();
_M_value /= __t;
return *this;
}
```

Again C99 types are used so the /= is compiled to call the C99 complex code, so again the fault is with clang++.

I could ignore clang++ and just use g++, however, the supported OS X compiler for C++ is clang++. For multi-precision I found that although std::complex<mpfr::mpreal> worked on Linux it failed dismally on Windows, so I produced a non-template complex class specific to mpfr::mpreal based on the non-specialised code in the GNU template library. So all I needed to do was produce a similar class for long double and then I found that I had opened a can worms.

The generic code for a complex number to a complex power used if the std::complex<T> is used with a T other than float, double or long double is:

```  template<typename _Tp>
inline complex<_Tp>
__complex_pow(const complex<_Tp>& __x, const complex<_Tp>& __y)
{ return __x == _Tp() ? _Tp() : std::exp(__y * std::log(__x)); }
```

This code is seriously flawed, if __x is zero the result is zero regardless of the value of __y, if __y has a negative real part and the imaginary part is zero the answer should be infinity (inf, nan) and if both parts are zero the answer should be not a number (nan, nan).

Elsewhere there is this:

```  // 26.2.5/13
// XXX: This is a grammar school implementation.
template<typename _Tp>
template<typename _Up>
complex<_Tp>&
complex<_Tp>::operator*=(const complex<_Up>& __z)
{
const _Tp __r = _M_real * __z.real() - _M_imag * __z.imag();
_M_imag = _M_real * __z.imag() + _M_imag * __z.real();
_M_real = __r;
return *this;
}

// 26.2.5/15
// XXX: This is a grammar school implementation.
std::cout << "1/0 = " << z_1 << std::endl;

template<typename _Tp>
template<typename _Up>
complex<_Tp>&
complex<_Tp>::operator/=(const complex<_Up>& __z)
{
const _Tp __r =  _M_real * __z.real() + _M_imag * __z.imag();
const _Tp __n = std::norm(__z);
_M_imag = (_M_imag * __z.real() - _M_real * __z.imag()) / __n;
_M_real = __r / __n;
return *this;
}
```

For multiplication the code fails if both values are infinity where at least one of the components is nan as the result is not a number instead of infinity. For division by infinity the result is non a number even if both compnents are +inf or -inf as the values of _M_real and _M_imag will be set to infinity divided by infinity which is not a number or not a number divided by not a number, the answer should be zero. It should be that much of a surprise that this code is faulty as it does say “This is a grammar school implementation.” as a comment for both routines.

This means that for “Combination 3” Saturn would fail to correctly produce the fractal image once multi-precision was required. As I had to correct the multi-precision complex class I also created a new class for long doubles. I’ve corrected all the methods that produced faulty results, the routines are:

real abs(complex)
real norm(complex)
complex pow(complex, complex)
complex pow(complex, real)
complex operator/=(complex)
complex operator*=(complex)
complex sqrt(complex)

where real is long double or mpfr::real and complex is mpfr::mpcomplex or ld::complex.

Since extra code had been added to correctly handle infinity and zero to the power zero I had expected an adverse affect on performance, happily I was wrong for most of the seed files I’ve tried so far the time to produce the image is shorter, there are some exceptions where it takes longer but the differences are small.

The fact that the code is quicker is not entirely unexpected, several years ago I had compared the performance of different versions of the g++ compiler (I’ve since lost figeres) the performance degraded with the introduction of a newer version I can’t remember which one but examining the assembler code produced by the newer slower version showed that it specifically called subroutines for simple addition, subtraction and multiplication of complex numbers, the older version operated on registers and did not call subroutines, the call overhead is the cause of the degraded performance. The new complex class specific to long doubles when compiled probably doesn’t call built in subroutines either.

Posted 5 July 2013 by in Fractal, Programming, Software

Tagged with , , , ,

While testing the version of Saturn I found some seed files that were not reproduced correctly. The following picture is an example: Quincunx

which the new version of Saturn reproduced as … Incorrect Quincunx

At least with fractal applications bugs manifest themselves visually. So, where did the black areas come from?

After checking the usual things like the formula and all the settings used were the same and that the image was correctly reproduced by version 3.0.1 it was clear that a debug session was required. The cause of the problem was unexpected, the fractal has an unnecessarily large bailout value and the iterations for the new black areas reach inf – nan i (inf is infinity and nan is not a number), the bailout condition uses the std::norm of the iteration value which is returned by clang++ as nan. Odd.

After further investigation I discovered that if I built Saturn using the g++ compiler instead of the clang++ compiler the seed file is reproduced correctly. The difference being that g++ returns inf for the std::norm of inf – nan i and not nan.

A side effect of nan being returned is that the bailout condition is not met and calculation continues using nan which always results in nan and is very slow, the locations where the problem occurs are coloured using the inner colouring method which is to colour the location black. To work around the problem the bailout condition is deemed to have been met if either the real or imaginary components are nan, inf or -inf. The change wasn’t sufficient, the time to complete the picture was much reduced but the black areas remained, this was because nan was added to the colour statistics and when the colour was determined nan was used to determine how close the colour was to one colour and the next in the colour map resulting incorrectly with black. To stop that happening the iteration value is only added to the colour statistics if is not nan.

I have since found that there are more problems with clang++ and std::complex that the fixes for this fractal do not address but that is for a future post.

Posted 2 July 2013 by in Fractal, Programming, Software

Tagged with , , , ,