0
\$\begingroup\$

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.

\$\endgroup\$
5
  • \$\begingroup\$ You have for y / for x / for iter nested loops. Please extract one or two helper functions. Also, dealing with a getchar() keystroke should clearly be broken out as a helper. We strive to make each function do one thing well. It should be possible to read and absorb the meaning of a function in a single screenful of code, with no scrolling around. No one has posted an Answer, so there's still time to revise the Question. Also, you might comment on the decision to not spawn N threads when running on an N-core machine. \$\endgroup\$
    – J_H
    Sep 7 at 23:36
  • \$\begingroup\$ I don't know what you mean by "dealing with a getchar() keystroke should clearly be broken out as a helper". I will separate the loops and the other function into separate code blocks. \$\endgroup\$
    – trxgnyp1
    Sep 7 at 23:41
  • \$\begingroup\$ Prior to the "separated code into multiple blocks" edit, there was a big switch (getchar()) { chunk of code in the middle of the main loop. Currently it is broken out as a handle_input() helper. Thank you, that is exactly the refactor I was encouraging you to take on. \$\endgroup\$
    – J_H
    Sep 8 at 5:53
  • \$\begingroup\$ Ah, okay. And yeah, you were right. \$\endgroup\$
    – trxgnyp1
    Sep 8 at 14:35
  • \$\begingroup\$ @trxgnyp1 Where is scaled_h declared as used in double real_y = (y_px / scaled_h) - 1.0;? Have you tried compiling what is posted here? \$\endgroup\$ Sep 14 at 15:57

2 Answers 2

2
\$\begingroup\$

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]);
}
\$\endgroup\$
9
  • \$\begingroup\$ Wow. Increasing the "raw" index of the framebuffer from the Y and X for loops is actually really smart. I might use that trick in other places of my OS, since I actually saw it was using multiplications for indexing in the disassembly. \$\endgroup\$
    – trxgnyp1
    Sep 9 at 13:57
  • \$\begingroup\$ Your goto trick doesn't work here because what I am checking is: Iterate N times, if the condition is false inside the loop, exit because we are not in the set; otherwise, continue checking. If we checked all iterations, and inside_set is true, we draw black and continue. Otherwise, draw HUE depending on iterations. What you are doing is: Iterate N times, if the condition is false (we are not in the set), set color to black and exit the loop. Otherwise, if we checked all iterations, we are inside the set, use HUE depending on iter... Except it's always going to be max_iter. \$\endgroup\$
    – trxgnyp1
    Sep 9 at 15:11
  • \$\begingroup\$ well - it seems I got the inside handling downside up in rev. 1. \$\endgroup\$
    – greybeard
    Sep 10 at 0:18
  • 1
    \$\begingroup\$ (Sorry for the edit rejection - it was too hot to think yesterday, and is the hour past midnight when I missed the corrected ,;.) \$\endgroup\$
    – greybeard
    Sep 10 at 0:26
  • 1
    \$\begingroup\$ (The differences in use of whitespace are deliberate, even where my spacing is inconsistent. E.g. one rule I usually follow is use blank lines to separate things only where they aren't separated visually by other means. In the inner loop, that's the reason there's no blank line between declarations and comment. The blank line for pronounced separation "before the x/y update" is inconsistent - works for me, with both decisions close to works both ways. Doubling zoom without spaces happens to be closer to how I' jot down the expression/formula.) \$\endgroup\$
    – greybeard
    Sep 10 at 1:16
2
\$\begingroup\$

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.

\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.