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.

Advertisements

Posted 19 August 2013 by element90 in Development, Programming

Tagged with , ,

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: