Archive for the ‘Programming’ Category

Another Unexpected Difference   3 comments


I received a query regarding the implementation of a fractal formula that didn’t produce the expected pictures. Here is such a picture. The implicated culprit was std::complex<T>, I’ve had trouble with this in the past leading to the use of my own version of the complex class in my software, it looks like it is entirely innocent.

I checked the formulae by trying them out in Saturn which also failed to produce the expected pictures. The program used for producing the fractals is called XoaS, the fractals in question aren’t part of XaoS but can be produced by setting the user formula.

An example of these formulae is shown below:

zn+1 = zn^(1 – (1/(log(z) + 0.1))) + 0.3

where z0 = the location in the complex plane

This is what is produced by the Window’s version of Xaos (the image is centred on (0,0) and has a width of 3):

XaoS.Windows

The Linux version of XaoS would be expected to produce the same image, it doesn’t (the image is centred on (0,0) and has a width 8):

XaoS.Linux

Saturn produces this (the image is centred on (0,0) and has width 32):

Saturn

Gnofract4d agrees with Saturn:

Gnofracta4d

Ultra Fractal also agrees with Saturn and Gnofract4d, the much smaller circle produced by XaoS for Linux must be to do with its bailout function, the bailout value is 4 which is the same as Gnofract4d and Saturn where the absolute value of z is used (abs(z) in Saturn, cabs in Gnofract4d).

The picture for the Linux version of XaoS took several attempts, initially there was just featureless black, zooming out I got a simple circle, but it was only when I went back to XaoS on Linux to check the location and size of the image that the blobs appeared.

I’ve now found that Saturn produces the same image as XaoS on Linux when the bailout function is changed to norm(z).

I can only conclude that the fractal images produced by the Windows version of XaoS are in error, as XaoS, Gnofract4d and Saturn on Linux and Ultra Fractal on Windows all produce the same image.

The images produced by Saturn for fractals with similar formulae are not all interesting so fractals of this sort won’t be added to Saturn. The erroreous images produced by XaoS on Windows on the the other hand are definitely worth adding to any fractal program provided the error can be found and replicated, indeed Saturn already has several fractal types found due to programming errors.

Advertisements

C++ Numeric Pitfalls   Leave a comment


C++ is a general perpose object orientated language and is by no means perfect. There are many pitfalls in its use, a “pitfall” is where the code in question compiles and links but when it is run it does something unexpected. For this post I’m restricting the scope to numerical functions, I fell into two of these traps through the course of developing the new version of Saturn and Titan.

The problem is best illustrated with a couple of pictures:

Saturn Version 4.0.0 (beta)

Saturn Version 4.0.0 (beta)


Saturn Version 3.0.1

Saturn Version 3.0.1

As part of the testing my software I load a great many seed files into the new software. A seed file is a PNG image with fractal parameters embedded saved by Saturn and they are useful for determing whether the new version still works as expected, in the case shown about the resulting image was certainly not expected as it shows the symptoms of not enough precision. The seed file is correctly reproduced using the previous version, so what’s changed?

In earlier testing I found problems with the std::complex standard template library regarding the handling of nan (not a number) and inf (infinity) when using the clang++ compiler. Those problems were fixed by implementing a new complex class for handling only long doubles which is when I found a pitfall. Somewhere in the class precision was being lost, I had noted an oddity in power of -1 determined using the new class verses the old class when complex values where used for both the value and the power, the result included a very small value for the imaginary part for both classes except that the std::complex value was the smaller of the two.

Part of the power function is to take the log of the input value, -1 results in PI in the imaginary part, so this test program was used to investigate:

#include <iostream>
#include <iomanip>
#include <complex>
#include "ldcomplex.h"

using namespace std;

std::string longdoubleToString(long double value)
{
    std::stringstream str;
    str << std::setprecision(std::numeric_limits<long double>::digits10) << value;
    return std::string(str.str());
}

int main(int argc, char** argv)
{
    std::complex<long double> a(-1, 0);
    ld::complex b(-1, 0);
    
    a = std::log(a);
    b = ld::log(b);
    
    std::cout << "std::complex<long double> log(-1 + 0i)" << std::endl;
    std::cout << "a.r = " << longdoubleToString(real(a)) << std::endl;
    std::cout << "a.i = " << longdoubleToString(imag(a)) << std::endl;
    std::cout << "ld::complex<long double> log(-1 + 0i)" << std::endl;
    std::cout << "b.r = " << longdoubleToString(real(b)) << std::endl;
    std::cout << "b.i = " << longdoubleToString(imag(b)) << std::endl;
    std::cout << "PIl = " << longdoubleToString(M_PIl) << std::endl;
    std::cout << "PI  = " << longdoubleToString(M_PI) << std::endl;

    return 0;
}

Running the test program:

std::complex<long double> log(-1 + 0i)
a.r = 0
a.i = 3.14159265358979324
ld::complex<long double> log(-1 + 0i)
b.r = 0
b.i = 3.14159265358979312
PIl = 3.14159265358979324
PI  = 3.14159265358979312

The value of PI for ld:complex differs from that of std::complex, the value produced by ld::complex matches the double value (M_PI) and not the long double value (M_PIl). So the loss of precision indicates that a double value is being introduced in the log function.

    inline complex log(const complex& z)
    {
        return complex(std::log(ld::abs(z)), arg(z));
    }

Note: the qualification of log and abs, arg is not qualified as it already in the ld namespace. I return to why the abs function is qualified, later, even tough it also in the ld namespace.

    inline long double arg(const complex& z) { return atan2(z.imag(), z.real()); }

The arg function is the cause of the problem, atan2 isn’t qualified and because it isn’t qualified the C function atan2 is called which is applied to doubles, the inputs are cast from long double to double and the output is cast from double to long double. The atan2 function in C has two sister versions atan2f for floats and atan2l for long doubles, in C++ when the std namespace is used all versions have the same name and the version called is determined by type.

Changing the code to:

    inline long double arg(const complex& z) { return std::atan2(z.imag(), z.real()); }

solves the problem. Checking the rest of class showed that this was the only occurance of an unqualified call to a numeric function.

Returning to the use of ld::abs function, it is qualified to make it obvious that function isn’t the C abs function. The C abs function has an integer input and an integer output, to perform the same function on float, double and long double there are functions fabsf, fabs and fabsl, C++ uses abs for float, double, long double and std::complex, the std:: standard namespace in not used and float, double or long double is used the call to abs casts the input value to integer and the output value to the required type. This NOT what is wanted and is the second pitfall.

When using numeric functions in C++ to get the right versions the std namespace MUST be used either by qualifying the name or via a using namespace statement.

When using the maths functions values in your code you may be ignorant of the fact that when the std namespace is not used all the functions (except abs) are casting to and from double and that the double versions of the functions are being called.

Posted 19 August 2013 by element90 in Development, Programming

Tagged with , ,

When is C++ long double not long double?   3 comments


When the compiler is Microsoft’s. The data type long double exists and can be used but there is no difference between long double and double.

Saturn and Titan have been developed on Linux using the g++ and clang++ compilers, long double occupies 128 bits (only 80 bits are used) and double occupies 64 bits. I used the long double datatype because it allowed deeper zooms than double before the generated picture degrades due to lack of precision. I introduced code to check the required precision to correctly generate fractals as the level of zoom increased, the initial shift from long double works most of the time, for some reason when long double is in fact double it doesn’t. Update: the example fractal illustrating the problem doesn’t shift to multi-precision because it is one of those fractals that the automatic mechanism doesn’t work, zooming in further with the Linux version exhibits the same problem it just does it later because of the greater precision of long doubles.

For Windows versions of Saturn and Titan the testing was cursory so I missed this problem which is present in all released versions of Saturn and Titan to date.

The problem is best illustrated using two screen shots of the yet to be released version 4.0.0 of Saturn, one running on Windows and the other on Linux.

Saturn 4.0.0 on Windows 7

Saturn 4.0.0 on Windows 7

Saturn 4.0.0 on Linux

Saturn 4.0.0 on Linux

So I now have to add code to Saturn and Titan to account for Microsoft defining long double as double.

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.

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

Tagged with , , , ,

An Unexpected Difference   Leave a comment


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

Quincunx

Quincunx

which the new version of Saturn reproduced as …

Incorrect Quincunx

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 element90 in Fractal, Programming, Software

Tagged with , , , ,

I Knew it was too easy – part 2   Leave a comment


At the end of the the post “I knew it was too easy” I mentioned that easiest method to get the parameter data out of Saturn seed files was to read the file directly and extract the keys and their values from the tEXt chunks. It was simple to implement and so I continued with development of version 4.0.0 of Saturn, today I noticed that when loading a seed file the wrong colour map was used, a colour map called Greyscale was unexpectedly used. The Greyscale colour map is a simple “auto” generated colour map which is defined as a single cycle sine wave for the red, green and blue components, the expected colour map was also an auto generated colour map. The problem turns out to be that the keys, “colour0”, “colour1” and “colour2” and their data were not extracted from the seed file, so that when the data was requested empty strings were returned and the default values were used instead which corresponds to the Greyscale colour map.

So why weren’t the keys and their data extracted from the seed file? Using a Linux utility called pngmeta all the keys and their data can be printed for any PNG file and that showed that the seed file did indeed contain the keys and their data for the “auto” generated colour map. But what exactly is in the file? Using a hex editor the answer is that Qt had embedded the colour map data using zTXt chunks and not tEXt chunks. The difference between zTXt and tEXt is that zTXt compresses the data but not the key, the data for the “colour0” was 46 characters, the length of the compressed data in the zEXt chunk is 120 bytes, the length of the data for “colour1” and “colour2” are 48 bytes each and the compressed size is again 120 bytes. This makes no sense, the compressed data is larger than the data to be compressed!

PNG files only support zlib compression (as far as I know). I initially looked for documentation on zlib which of course would have a C interface. Saturn is written in C++ so some method of accessing zlib’s functionality from C++ would be preferred, my initial thought was “Does Boost handle zlib decompression?”, after a bit of searching the answer is, “Yes, it does”.

I loaded the compressed data into a std::vector and then passed it to a utility routine to convert it into std::string:

std::string uncompress(std::vector<char>& buffer)
{
    std::string str;
    io::filtering_ostream os;
    os.push(io::zlib_decompressor());
    os.push(io::back_inserter(str));
    int size = buffer.size();
    for (int i = 0; i < size; i++)
    {
        os.put(buffer[i]);
    }
    // the above loop could be replaced using
    //
    // os.write(reinterpret_cast<char*>(&buffer[0]), buffer.size());
    //
    // or
    //
    // os.write(buffer.data(), buffet.size()); // C++11
    //
    return str;
}

A fragment of code on its own, usually results in questions such as “Where is io defined?” and “What headers should I include?”, so here are the missing pieces required to compile the code:

#include <boost/iostreams/device/back_inserter.hpp>
#include <boost/iostreams/filtering_stream.hpp>
#include <boost/iostreams/filter/zlib.hpp>

namespace io = boost::iostreams;

At one time the different path delimiters used on Windows and Unix/Linux meant that different include statements had be used and were selected by condition compilation, fortunately the Microsoft C++ compiler I use happily accepts Unix style delimiters.

I’ve used some C++11 features in my code, I would use more except that compiler support is not yet complete and as I release Saturn and Titan for Linux and Windows I have to take account of the compiler that has the least support and that is Microsoft’s. The sooner better C++11 support is more widespread the better as much cleaner code can written.

So now all the embedded data can be extracted from Saturn seed files, so the next release is one step closer.

Posted 8 May 2013 by element90 in Programming, Software

Tagged with , , , ,

%d bloggers like this: