3

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!

kamelfanger83
  • 97
  • 1
  • 7
  • 2
    Do you want [imaginary numbers](https://en.cppreference.com/w/c/numeric/complex/imaginary)? Also consider [`std::complex`](https://en.cppreference.com/w/cpp/numeric/complex). – tadman Jan 08 '21 at 09:49
  • 5
    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 Jan 08 '21 at 10:00
  • @Damien the [Mandelbrot set iterative function](https://en.wikipedia.org/wiki/Mandelbrot_set) is fairly well known. – Patrick Roberts Jan 08 '21 at 10:02
  • 2
    @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 Jan 08 '21 at 10:10
  • @Damien indeed, using a program I wrote several years ago. Nostalgic memories :) – Patrick Roberts Jan 08 '21 at 10:12
  • 2
    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 Jan 08 '21 at 10:30

2 Answers2

6

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)
Patrick Roberts
  • 49,224
  • 10
  • 102
  • 153
  • 1
    Thanks a lot! But I still have to major questions: 1. How can you call the function if you want the imaginary part to be a variable. You can't just add an "i" because that would change the name of the variable. 2. What does "std::norm(z) < 4" do? – kamelfanger83 Jan 08 '21 at 10:21
  • @kamelfanger83 good question, the [`std::complex`](https://en.cppreference.com/w/cpp/numeric/complex/complex) constructor accepts the real and imaginary components as arguments, i.e. `std::complex c(0.2, 1);` will work. – Patrick Roberts Jan 08 '21 at 10:23
  • @kamelfanger83 `std::abs(z)` is basically `std::sqrt(std::norm(z))` so obviously it's far more costly – phuclv Jan 08 '21 at 10:23
  • @kamelfanger83 and `std::norm(z)` computes the absolute square of `z`, which is the absolute value squared. It's cheaper to compute because it does not require use of [`std::sqrt()`](https://en.cppreference.com/w/cpp/numeric/math/sqrt). – Patrick Roberts Jan 08 '21 at 10:25
  • Ooh thank you very much! It works just fine now! But it seems to work a lot slower than it did before. Well of course it didn't even work then but it nearly worked ;). So if we would be able to achieve the multiplication of a double by 2 without chanching the sign of it I think that would be a big computing time benefit. Did you understand what I was trying to do with my original code? – kamelfanger83 Jan 08 '21 at 10:27
  • 1
    @kamelfanger83 _"But it seems to work a lot slower than it did before."_ Have you enabled optimizations on your compiler? Add `-O3` to your compilation flags and in theory you should have the same performance as without using `std::complex`. And yes, your original code is just an explicit complex multiply-add of the separated real and imaginary components. It's exactly the same as what `std::complex` does. – Patrick Roberts Jan 08 '21 at 15:06
0

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.

kamelfanger83
  • 97
  • 1
  • 7