Multiplications with doubles c++

3
votes

I'm new to C++ and I wanted to code the mandelbrot set in it to get some practise and compare the speed to Python. My function to iterate to see if the complex number blows up to infinity looks like this:

int iterate(double ir, double ii, int iterations){
    double sr = 0.0;
    double si = 0.0;
    unsigned int steps = 0;
    while(sr < 2.0 && sr > -2.0 && si < 2.0 && si > -2.0 && steps < iterations){
        sr = sr*sr-si*si+ir;
        si = si*sr*2 + ii;
        steps ++;
    }
    return steps;
}

I compared the outputs to my python code and realised the the line where it updates the imaginary part doesn't work properly. The *2 doesn't work properly. Insted of doubling it makes the results negative. For example when I insert the number 0.2 + 1i it does the following:

0.2 + 1i, -0.76 - 0.52i, 0.5072 + 0.472512i, 0.233984 + 1.22112i, -1.23639 - 2.01956i

What it should do instead (what I got from my Python programm) is:

0.2+ 1i, -0.76 + 1.4i, -1.1823999999999997 - 1.1279999999999997i, 0.3256857599999999 + 3.6674943999999985i

My guess is that it writes to the sign instead of either doubling the mantissa or increasing the exponent. I can't just make it unsigned because it must also be able to store negative values. What can you do about it? Thanks for your awnsers!

2
c++mandelbrot
Do you want imaginary numbers? Also consider std::complex. - tadman
You did not recall what formula you want to implement. Note that your code uses the newly calculated sr to calculate si. Is it what you want? - Damien
@Damien the Mandelbrot set iterative function is fairly well known. - Patrick Roberts
@PatrickRoberts It is well known by people who know it ! :) I was sure that I could find it by a google search, but a general advice to OP is to provide all relevant information in the post. That said, thank you for the link! Did you use this formula for your avatar? - Damien
The problem is that when you compute si, you have modified sr. For instance, your first iteration uses -0.76 instead of 0.2 as the real part when computing the imaginary part. - molbdnilo

2 Answers

6
votes

I would write your function like this, using std::complex:

template <class T>
std::size_t iterate(std::complex<T> c, std::size_t max_iterations) {
    std::complex<T> z = 0;
    std::size_t steps = 0;

    while (std::norm(z) <= 4 && steps < max_iterations) {
        z = z * z + c;
        ++steps;
    }

    return steps;
}

Then you can call your function like this:

using std::complex_literals::operator""i;
std::complex<double> c = 0.2 + 1i; // or c(0.2, 1);
const auto steps = iterate(c, 256);

The reason you should use std::norm(z) <= 4 instead of std::abs(z) <= 2 is because it's less costly to compute.

You can verify that the following program outputs your expected sequence if you stream z to std::cout after each iteration:

(0.2,1)(-0.76,1.4)(-1.1824,-1.128)(0.325686,3.66749)
0
votes

Uhm I know I'm a little late but I think it's worth it: I fixed my problem without the library. It actually was quite a dumb mistake: I took the new real value for the calculation of the new imaginary value instead of the previous one. I fixxed it by adding a "buffer" variable that stores the previous real value. My current code is:

int iterate(double r,double  i, int max_iterations) {
    double ir = 0;
    double nir;
    double ii = 0;
    int steps = 0;
    while (ir*ir+ii*ii <= 4 && steps < max_iterations) {
        nir = ir*ir-ii*ii+r;
        ii = 2*ir*ii + i;
        ir = nir;
        ++steps;
    }

    return steps;
}

It works perfectly fine and is a lot faster than the solution with the library. I estimate it's about 5-10 times faster which is the reason why I decided to share this. If I compile it in a release build it pretty much renders instantly using the windows.h header for my graphics.