Julia: one rule, infinite starting points in C

/* Arithmos: draw the filled Julia set using only standard C and libm.
 *
 * Compile: cc -std=c11 -O2 julia.c -lm -o julia
 * Run:     ./julia
 * Options: ./julia WIDTH HEIGHT ITERATIONS OUTPUT.ppm
 * Example: ./julia 1600 1000 300 julia.ppm
 *
 * The output is a binary PPM (P6) image. It needs no graphics library.
 * Dark pixels have not escaped within the iteration limit; that alone
 * does not prove they belong to the filled Julia set.
 * The constant below gives the Douady rabbit; each pixel starts its own orbit.
 */
#include <errno.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>

/* Read a positive integer, rejecting malformed or excessive values. */
static int read_int(const char *text, int minimum, int maximum) {
    char *end;
    errno = 0;
    long value = strtol(text, &end, 10);
    if (errno || end == text || *end || value < minimum || value > maximum) {
        fprintf(stderr, "Expected an integer from %d to %d: %s\n",
                minimum, maximum, text);
        exit(EXIT_FAILURE);
    }
    return (int)value;
}

/* Interpolate through the website's mint, blue, rose and gold palette. */
static void gradient(double smooth, unsigned char rgb[3]) {
    const double stops[4][3] = {
        {91, 227, 201}, {128, 146, 240},
        {218, 138, 220}, {244, 200, 127}
    };
    double t = fmod(log1p(fmax(0.0, smooth)) * 0.35, 1.0) * 3.0;
    int band = (int)t;
    double blend = t - band;
    double brightness = fmin(1.0, 0.11 + smooth / 12.0);
    for (int channel = 0; channel < 3; ++channel) {
        double color = stops[band][channel] * (1.0 - blend)
                     + stops[band + 1][channel] * blend;
        rgb[channel] = (unsigned char)lround(color * brightness);
    }
}

int main(int argc, char **argv) {
    if (argc != 1 && argc != 5) {
        fprintf(stderr, "Usage: %s [WIDTH HEIGHT ITERATIONS OUTPUT.ppm]\n", argv[0]);
        return EXIT_FAILURE;
    }
    int width = argc == 5 ? read_int(argv[1], 2, 4096) : 1200;
    int height = argc == 5 ? read_int(argv[2], 2, 4096) : 800;
    int limit = argc == 5 ? read_int(argv[3], 1, 2000) : 250;
    const char *filename = argc == 5 ? argv[4] : "julia.ppm";
    FILE *image = fopen(filename, "wb");
    if (!image) { perror(filename); return EXIT_FAILURE; }
    if (fprintf(image, "P6\n%d %d\n255\n", width, height) < 0) {
        fclose(image); return EXIT_FAILURE;
    }

    /* Choose a view of the complex plane; keep both axes at one scale. */
    const double center_real = 0.0, center_imag = 0.0;
    const double parameter_real = -0.1225611668766536;
    const double parameter_imag = 0.7448617666197442;
    double span_x = fmax(3.4, 3.0 * width / height);
    double span_y = span_x * height / width;

    for (int py = 0; py < height; ++py) {
        for (int px = 0; px < width; ++px) {
            /* Fix c; each pixel supplies a different starting value of z. */
            double zr = center_real + ((px + 0.5) / width - 0.5) * span_x;
            double zi = center_imag + (0.5 - (py + 0.5) / height) * span_y;
            double cr = parameter_real, ci = parameter_imag;
            int step = 0;

            /* z <- z*z + c, written in real and imaginary parts.
             * (a + bi)^2 = (a*a - b*b) + (2*a*b)i
             * Preserve the old real part until both updates are ready.
             */
            while (step < limit && zr * zr + zi * zi <= 4.0) {
                double next_real = zr * zr - zi * zi + cr;
                zi = 2.0 * zr * zi + ci;
                zr = next_real;
                ++step;
            }

            unsigned char rgb[3] = {9, 16, 27};
            if (zr * zr + zi * zi > 4.0) {
                /* Smooth the escape count to avoid hard color bands. */
                double magnitude = sqrt(zr * zr + zi * zi);
                double smooth = fmax(0.0, step + 1.0 - log2(log2(magnitude)));
                gradient(smooth, rgb);
            }
            if (fwrite(rgb, 1, 3, image) != 3) {
                perror("Writing image"); fclose(image); return EXIT_FAILURE;
            }
        }
    }
    if (fclose(image) != 0) { perror("Closing image"); return EXIT_FAILURE; }
    printf("Wrote %s (%d x %d, %d iterations)\n", filename, width, height, limit);
    return EXIT_SUCCESS;
}