Search code examples
c++mandelbrot

Multiplications with doubles c++


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!


Solution

  • 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)