floating point tests with NTL lib

test.cpp

/*
To comppile it on Debian, first install the package libntl-dev:
  sudo apt-get install libntl-dev
then you can compile and run it like this:
  g++ test.cpp -lntl -o test ; ./test

output on my system:

mandelbrot calculation at point (-0.025, 0.643265)

float:
height: 430, mag: 12.3007

double:
height: 1233, mag: 8.95292

RR, 53 bits:
height: 1233, mag: 8.952915868

RR, 80 bits:
height: 1423, mag: 7.389421261

RR, 1000 bits:
height: 1423, mag: 7.249571922057591375674421545391787210883200832900173534966107998908019153055426981005408394157085535
*/

#include <NTL/RR.h>
#include <iostream>

using namespace std;
using namespace NTL;

template<class T>
uint32_t mbrot(T c_r, T c_i, uint32_t bail_out, T& mag)
{

    /* point c belongs to the Mandelbrot set if and only if
     * the magnitude of the f(c) <= 2.0 */
    uint32_t height = 0;
    T zr(0.0);
    T zi(0.0);
    T tmp_r, tmp_i;
    mag = 0.0;

    while ( ( height < bail_out ) && ( mag < 4.0 ) ) {
        tmp_r = ( zr * zr ) - ( zi * zi );
        tmp_i = ( zr * zi ) + ( zr * zi );
        zr = tmp_r + c_r;
        zi = tmp_i + c_i;
        mag = zr * zr + zi * zi;
        height += 1;
    }

    return height;

}

template<class T>
void test(double c_r, double c_i)
{
    T mag;
    uint32_t height = mbrot<T>(T(c_r), T(c_i), 2048, mag);
    cout << "height: " << height << ", mag: "  << mag << endl << endl;
}

int main()
{
    double c_r(-0.025);
    double c_i(0.64326485);

    cout << "mandelbrot calculation at point (" << c_r << ", " << c_i << ")" << endl << endl;

    cout << "float:" << endl;
    test<float>(c_r, c_i);

    cout << "double:" << endl;
    test<double>(c_r, c_i);

    cout << "RR, 53 bits:" << endl;
    RR::SetPrecision(53);
    test<RR>(c_r, c_i);

    cout << "RR, 80 bits:" << endl;
    RR::SetPrecision(80);
    test<RR>(c_r, c_i);

    cout << "RR, 1000 bits:" << endl;
    RR::SetPrecision(1000);
    RR::SetOutputPrecision(100);
    test<RR>(c_r, c_i);

    return 0;
}