Improving performance of mandelbrot set calculation

Question

I am making a hobby OS, and I thought about adding a command for interactively rendering the Mandelbrot set. The "interactive" part is not really important, but I wanted to check if the Mandelbrot algorithm is as good as it can be.

The most performance-sensitive part is the code inside the inner for loop, the one that iterates max_iter.

I tried making it as optimized as possible, but I might have left something (also note that my OS doesn't support multicore yet).

/* Smaller -> More zoom (1.0 is defaul) */
#define DEFAULT_ZOOM 1.0

/* Offsets for moving in the set. */
#define DEFAULT_X_OFF 0.0 /* -2.0..+2.0 (Left..Right) */
#define DEFAULT_Y_OFF 0.0 /* -1.0..+1.0 (Up..Down) */

/* Default steps when moving around */
#define ZOOM_STEP         0.5
#define DEFAULT_MOVE_STEP 0.1

int mandelbrot(int argc, char** argv) {
    const uint32_t w   = screen_width();
    const uint32_t h   = screen_height();
    const int max_iter = atoi(argv[1]);

    /* Check arguments are valid... */

    /* Framebuffer. Each 32 bit entry is a color */
    volatile uint32_t* fb = fb_get_ptr();

    /* Calculate some values here for performance */
    const double scaled_h = h / 2.0;
    const double scaled_w = w / 3.0;

    /* Will become smaller when zooming */
    double zoom     = DEFAULT_ZOOM;
    double x_offset = DEFAULT_X_OFF;
    double y_offset = DEFAULT_X_OFF;

    /* Will be scaled when zooming, so we don't move too much */
    double move_step = DEFAULT_MOVE_STEP;

    bool main_loop = true;
    while (main_loop) {
        do_mandelbrot();
        handle_input();
    }
}

My do_mandelbrot() function:

static void do_mandelbrot(void) {
    for (uint32_t y_px = 0; y_px < h; y_px++) {
        double real_y = (y_px / scaled_h) - 1.0;
        real_y *= zoom;
        real_y += y_offset;

        for (uint32_t x_px = 0; x_px < w; x_px++) {
            double real_x = (x_px / scaled_w) - 2.0;
            real_x *= zoom;
            real_x += x_offset;

            double x = real_x;
            double y = real_y;

            bool inside_set = true;

            int iter;
            for (iter = 0; iter < max_iter; iter++) {
                /* Calulate squares once */
                double sqr_x = x * x;
                double sqr_y = y * y;

                /* Absolute value of a complex number is the distance from
                 * origin: sqrt(x^2 + y^2) > 2 */
                if ((sqr_x + sqr_y) > 2 * 2) {
                    inside_set = false;
                    break;
                }

                y = (2.0 * x * y) + real_y;
                x = (sqr_x - sqr_y) + real_x;
            }

            /* If it's inside the set, draw black */
            if (inside_set) {
                fb[y_px * w + x_px] = 0x000000;
                continue;
            }

            /* Get 0..360 hue for color based on iter..max_iter */
            int scaled_hue      = iter * MAX_H / max_iter;
            fb[y_px * w + x_px] = hue2rgb(scaled_hue);
        }
    }
}

My hue2rgb function:

static uint32_t hue2rgb(float h) {
    float prime = fmod(h / 60.f, 6);
    float x     = 1 - fabs(fmod(prime, 2) - 1);

    uint32_t ret = 0x000000;

    if (prime >= 0 && prime < 1) {
        ret |= 0xFF0000;                 /* r = 255 */
        ret |= (uint8_t)(x * 0xFF) << 8; /* g = x */
    } else if (prime < 2) {
        ret |= 0x00FF00;                  /* g = 255 */
        ret |= (uint8_t)(x * 0xFF) << 16; /* r = x */
    } else if (prime < 3) {
        ret |= 0x00FF00;                 /* g = 255 */
        ret |= (uint8_t)(x * 0xFF) << 0; /* b = x */
    } else if (prime < 4) {
        ret |= 0x0000FF;                 /* b = 255 */
        ret |= (uint8_t)(x * 0xFF) << 8; /* g = x */
    } else if (prime < 5) {
        ret |= 0x0000FF;                  /* b = 255 */
        ret |= (uint8_t)(x * 0xFF) << 16; /* r = x */
    } else if (prime < 6) {
        ret |= 0xFF0000;                 /* r = 255 */
        ret |= (uint8_t)(x * 0xFF) << 0; /* b = x */
    }

    return ret;
}

My handle_input function should not be important, but here it is:

static void handle_input() {
    /* User input, not really important for performance */
    switch (getchar()) {
        case KEY_ZOOM_IN:
            zoom *= ZOOM_STEP;
            move_step *= ZOOM_STEP;
            break;
        case KEY_ZOOM_OUT:
            zoom /= ZOOM_STEP;
            move_step /= ZOOM_STEP;
            break;
        case KEY_UP:
            if (y_offset > -1.1)
                y_offset -= move_step;
            break;
        case KEY_DOWN:
            if (y_offset < 1.1)
                y_offset += move_step;
            break;
        case KEY_LEFT:
            if (x_offset > -2.1)
                x_offset -= move_step;
            break;
        case KEY_RIGHT:
            if (x_offset < 2.1)
                x_offset += move_step;
            break;
        case KEY_RESET:
            zoom      = DEFAULT_ZOOM;
            move_step = DEFAULT_MOVE_STEP;
            x_offset  = DEFAULT_X_OFF;
            y_offset  = DEFAULT_X_OFF;
            break;
        default:
            break;
    }
}

Edit: I changed the hue2rgb return value.

Edit: Separated code into multiple blocks. Didn't bother with variables or returns, contents of do_mandelbrot and handle_input should be replaced where the functions are called.

Answer 1

There is a lot of literature on performant 2D visualisations of the Mandelbrot set, including wikipedia - I intend not to go into algorithmic improvements.

I expect most of what I think of to improve do_mandelbrot()'s speed is covered by contemporary "optimising" compilers, and may be at odds with code the way you think about the solution.
An inlining compiler would even cover using return for lack of continue a loop enclosing the current one - I suggest using goto instead of introducing a flag to re-establish why the innermost loop was terminated.
I consider goto a red flag rather than a problem: don't use it to code in an unstructured way.

I prefer declaring identifiers I want to have the same type for the same reason in a single declaration.
Switching the outer loops allows just walking the frame buffer with an index or pointer.

static const uint32_t BLACK = 0;

static void do_mandelbrot(void) {
    const double scale_h = zoom / scaled_h,
                 scale_w = zoom / scaled_w,
                 hue_scale = MAX_H / max_iter;
    for (uint32_t x_px = 0, fb_index = 0; x_px < w; x_px++) {
        const double real_x = x_px * scale_w - 2*zoom;
    pixel_loop:
        for (uint32_t y_px = 0 ; y_px < h; y_px++, fb_index++) {
            const double real_y = y_px * scale_h - zoom,

            double x = real_x,
                   y = real_y;

            for (int iter = 0; iter < max_iter; iter++) {
                /* Calculate squares once */
                const double sqr_x = x * x,
                             sqr_y = y * y;
                /* Absolute value of a complex number is the
                 * distance from origin: sqrt(x^2 + y^2) > 2 */
                if (sqr_x + sqr_y > 2 * 2) {  /* outside set */
                   /* scale to 0..360 hue for color 
                    * based on ratio of iter to max_iter */
                    fb[fb_index] = hue2rgb(iter * hue_scale);  
                    goto next_pixel;  /* continue pixel_loop; */
                }

                y = (2.0 * x * y) + real_y;
                x = (sqr_x - sqr_y) + real_x;
            }
            fb[fb_index] = BLACK;  /* inside set */
        next_pixel:
        }
    }
}

hue2rgb() shows avoidable repetition in handling x:

static uint32_t hue_to_rgb(float h) {
    float prime = fmod(h / 60.f, 6);

    uint32_t rgb;
    int shift;

    switch (floor(prime)) {
    default: return 0xffffff;  /* not BLACK (inside) */
    case 0:
        rgb = 0xFF0000;  /* r = 255 *//* introduce symbolic names? */
        shift = 8;       /* g */
        break;
    case 1:
        rgb = 0x00FF00;  /* g = 255 */
        shift = 16;      /* r */
        break;
    case 2:
        rgb = 0x00FF00;  /* g = 255 */
        shift = 0;       /* b */
        break;
    case 3:
        rgb = 0x0000FF;  /* b = 255 */
        shift = 8;       /* g */
        break;
    case 4:
        rgb = 0x0000FF;  /* b = 255 */
        shift = 16;      /* r */
        break;
    case 5:
        rgb = 0xFF0000;  /* r = 255 */
        shift = 0;       /* b */
        break;
    }

    float x = 1 - fabs(fmod(prime, 2) - 1);
    return rgb | ((uint8_t)(x * 0xFF) << shift);
}
/* can be replaced by array access - less readable/mode difficult to comment */
static const uint32_t HUE2RGB[] = {
        0xFF0000,  /* r = 255 */
        0x00FF00,  /* g = 255 */
        0x00FF00,  /* g = 255 */
        0x0000FF,  /* b = 255 */
        0x0000FF,  /* b = 255 */
        0xFF0000,  /* r = 255 */
    };
static const int HUE_SHIFT[] = {
        shift = 8;  /* g */
        shift = 16; /* r */
        shift = 0;  /* b */
        shift = 8;  /* g */
        shift = 16; /* r */
        shift = 0;  /* b */
    };
static uint32_t hue2rgb(float h) {
    float prime = fmod(h / 60.f, 6);
    if (prime < 0 || 6 <= prime)
        return 0xffffff;
    float x     = 1 - fabs(fmod(prime, 2) - 1);
    int i_prime = floor(prime);

    return HUE2RGB[i_prime] | ((uint8_t)(x * 0xFF) << HUE_SHIFT[i_prime]);
}

Answer 2

Avoid double functions for float problems

OP's code promotes the float to a double, calls the double function and then narrows the double result to a float. Better to call a float function and skip the promotion and narrowing.

Review function calls. Examples:

// float prime = fmod(h / 60.f, 6);
// float x     = 1 - fabs(fmod(prime, 2) - 1);
float prime = fmodf(h / 60.f, 6);
float x     = 1 - fabsf(fmodf(prime, 2) - 1);

Tip: Save time and enable all warnings as a well enabled compiler would advise about this.

Consider integer math

The if (prime) tree does a lot of FP compares against integers values. Consider converting prime to an int to perform the multiple compares -perhaps even a switch(). Profile to see if this improved performance.