More on an Unexpected Difference   2 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

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.

Advertisements

Posted 5 July 2013 by element90 in Fractal, Programming, Software

Tagged with , , , ,

2 responses to “More on an Unexpected Difference

Subscribe to comments with RSS.

  1. Hi Mark.

    My name is Chris M. Thomasson and am having a problem with getting a
    couple of fractal equations to properly compute using native `std::complex’
    functionality:

    F(Z) = Z^(1.0 – (1 / (LOG(Z) + 0.1))) + 0.283079

    F(Z) = Z^(1.3 – (1 / (LOG(Z) + 0.11))) + 0.42457

    The fractals should produce a hyperbolic 3d polyhedral-like entity in which
    Cantor-like bifurcations and Glynn spirals are embedded within each
    polygonal surface. The following renderings, created with Xaos, create a
    fairly interesting fractal:

    However, with current MSVC++, I am getting extreme voids where there
    should be vast detail. Before, I have to create my own complex class, I was
    wondering if you could render these equations with the custom/correct
    complex class you describe here:

    https://element90.wordpress.com/2013/07/02/an-unexpected-difference

    and get something like Xaos produces.

    If these render correctly for you, I am going to have to follow the exact
    same
    road you followed wrt buggy `std::complex’… Damn!

    ;^/

    Here is a link to my site wrt Fractals:

    https://plus.google.com/101799841244447089430/posts

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: