Newton Fractal writen in C

Question

The code has no errors, it is just slow. It takes 1.19s to start.

How can I reduce runtime without using threads or multiple processes?

I have tried compiler optimization

(-o1, -o2 and -o3)

These did not do much.

typedef struct s_float2{
    double  x;
    double  y;
}   t_float2;

static t_float2 next(t_float2 z)
{
    t_float2    new;
    double      temp_deno;

    temp_deno = 3 * (z.x * z.x + z.y * z.y) * (z.x * z.x + z.y * z.y);
    new.x = ((z.x * z.x * z.x * z.x * z.x + 2 * z.x * z.x * z.x * z.y * z.y
                - z.x * z.x + z.x * z.y * z.y * z.y * z.y + z.y * z.y)
            / temp_deno);
    new.y = ((z.y * (z.x * z.x * z.x * z.x + 2 * z.x * z.x
                    * z.y * z.y + 2 * z.x + z.y * z.y * z.y * z.y))
            / temp_deno);
    z.x -= new.x;
    z.y -= new.y;
    return (z);
}

static t_float2 find_diff(t_float2 z, t_float2 root)
{
    t_float2    new;

    new.x = z.x - root.x;
    new.y = z.y - root.y;
    return (new);
}

void    init_roots(t_float2 *roots)
{
    static unsigned int flag;

    if (!flag)
    {
        roots[0] = (t_float2){1, 0};
        roots[1] = (t_float2){-0.5, 0.866025403784438};
        roots[2] = (t_float2){-0.5, -0.866025403784438};
        flag = 1;
    }
}

unsigned int    checker_tolerance(t_fractal *fractal, t_float2 z,
    unsigned int iterations)
{
    t_float2        difference;
    static t_float2 roots[3];

    init_roots(&roots);
    difference = find_diff(z, roots[0]);
    if (fabs(difference.x) < fractal->tolerance
        && fabs(difference.y) < fractal->tolerance)
        return (1);
    difference = find_diff(z, roots[1]);
    if (fabs(difference.x) < fractal->tolerance
        && fabs(difference.y) < fractal->tolerance)
        return (1);
    difference = find_diff(z, roots[2]);
    if (fabs(difference.x) < fractal->tolerance
        && fabs(difference.y) < fractal->tolerance)
        return (1);
    return (0);
}

void    calculate_newton(t_fractal *fractal)
{
    t_float2        z;
    unsigned int    iterations;

    iterations = 0;
    z.x = (fractal->x / fractal->zoom) + fractal->offset_x;
    z.y = (fractal->y / fractal->zoom) + fractal->offset_y;
    while (iterations < fractal->max_iterations)
    {
        z = next(z);
        if (checker_tolerance(fractal, z, iterations) == 1)
            return (put_color_to_pixel(fractal, fractal->x, fractal->y,
                    fractal->color * iterations / 2));
        ++iterations;
    }
    put_color_to_pixel(fractal, fractal->x, fractal->y, 0x000000);
}

The github link is https://github.com/MehdiMirzaie2/fractol.

Also, if you just want the main here it is:

# include <mlx.h>
# include <stdlib.h>
# include <unistd.h>

typedef struct s_fractal
{
    void    *mlx;
    void    *window;
    void    *image;
    void    *pointer_to_image;
    int     bits_per_pixel;
    int     size_line;
    int     endian;
    int     x;
    int     y;
    double  zx;
    double  zy;
    double  cx;
    double  cy;
    int     color;
    double  offset_x;
    double  offset_y;
    double  zoom;
    int     name;
    int     max_iterations;
    float   tolerance;
}           t_fractal;

static void free_and_explain(t_fractal *fractal)
{
    ft_putendl_fd("Available fractals: a = mandel, b = julia, c = newton\
    \n\033[1;31mASK FOR EXTRA, check your program!\033[0m", 1);
    exit_fractal(fractal);
}

int draw_fractal(t_fractal *fractal, double cx, double cy)
{
    fractal->x = 0;
    fractal->y = 0;
    while (fractal->x < SIZE)
    {
        while (fractal->y < SIZE)
        {
            if (fractal->name == 'a')
                calculate_mandelbrot(fractal);
            else if (fractal->name == 'b')
                calculate_julia(fractal, cx, cy);
            else if (fractal->name == 'c')
                calculate_newton(fractal);
            else
                free_and_explain(fractal);
            fractal->y++;
        }
        fractal->x++;
        fractal->y = 0;
    }
    mlx_put_image_to_window(fractal->mlx, fractal->window, fractal->image, 0,
        0);
    return (0);
}

int main(int argc, char **argv)
{
    t_fractal   *fractal;

    if (argc != 2 || argv[1][1] != '\0')
    {
        ft_putendl_fd("Available fractals: a = mandel, b = julia, c = newton\
        \n\033[1;31mASK FOR EXTRA!\033[0m", 1);
        return (0);
    }
    fractal = malloc(sizeof(t_fractal));
    init_fractal(fractal);
    init_mlx(fractal);
    fractal->name = argv[1][0];
    mlx_key_hook(fractal->window, key_hook, fractal);
    mlx_mouse_hook(fractal->window, mouse_hook, fractal);
    mlx_hook(fractal->window, 17, 0L, exit_fractal, fractal);
    draw_fractal(fractal, -0.79, 0.15);
    mlx_loop(fractal->mlx);
    return (0);
}
```

Answer 1

I don't like parentheses around expressions returned.

I do not understand init_roots(). Why not just initialise checker_tolerance()'s roots, or make it a compilation unit static const?
(Do the numbers have a meaning? ±sqrt(¾), sin(±60°)?)

Do not repeat yourself:

unsigned int
checker_tolerance(t_fractal const *fractal, t_float2 z, unsigned int iterations)
{
    static const t_float2 roots[3] = {
            (t_float2){1, 0},
            (t_float2){-0.5, 0.866025403784439},
            (t_float2){-0.5, -0.866025403784439}
        };
    const int nroots = sizeof roots / sizeof *root;

    for (int root = 0 ; root < nroots ; root++) {
        t_float2 difference = find_diff(z, roots[root]);
        if (fabs(difference.x) < fractal->tolerance
            && fabs(difference.y) < fractal->tolerance)
        return 1;

    return 0;
}

– still looks weird:
make roots a vector parameter? (Darn - C, not C++, and array/pointer & length sucks.)

In next(), I need to count "multiplicities" - I'd rather not:

static t_float2 next(t_float2 z) // const t_float2 *z?
{
    const double
        x2 = z.x * z.x, // x4 = x2 * x2, // x3 = x2 * z.x, x5 = x2 * x3,
        y2 = z.y * z.y, // y4 = y2 * y2,
        sum_of_squares = x2 + y2,
        sum_of_squares2 = sum_of_squares * sum_of_squares, // x4 + 2*x2*y2 + y4
        denominator = 3 * sum_of_squares,
        //   z.x * z.x * z.x * z.x * z.x + 2 * z.x * z.x * z.x * z.y * z.y
        // - z.x * z.x                   + z.x * z.y * z.y * z.y * z.y + z.y * z.y
        x = // (       x5 + 2 *  x3 * y2 -     x2  + z.x * y4  + y2) / denominator,
            // (z.x * (x4 + 2 *  x2 * y2 -     z.x +       y4) + y2) / denominator,
               (z.x * (sum_of_squares2 - z.x) + y2) / denominator,
        // z.y * (z.x * z.x * z.x * z.x + 2 * z.x * z.x
        //        * z.y * z.y + 2 * z.x + z.y * z.y * z.y * z.y)
        y = // (z.y * (x4 + 2 *  x2 * y2 + 2 * z.x +       y4) / denominator;
                z.y * (sum_of_squares2 + 2 * z.x) / denominator;
    z.x -= x;
    z.y -= y;
    return z;
}

Oh, look: I think I recognized how similar the modifications to x and y have been.

Seeing how fractal in calculate_newton(t_fractal *fractal) gets used, consider a const t_fractal *fractal).
And return ‹expression› looks weird in a void function even with a void expression.

Answer 2

Use complex numbers

Since C99, C supports complex numbers, and has many mathematical functions that work on complex numbers. I would rewrite all your code to use that. For example, for a Newton fractal with \$p(z) = z^3 - 1\$:

static complex double next(complex double z)
{
    complex double denominator = 3 * z * z;
    complex double numerator = z * z * z - 1;
    return z - numerator / denominator;
}

This will make your code look much nicer. It may or may not be more performant, but that brings me to:

Use -ffast-math and -march=native

While in pure mathematics, multiplication of real and complex numbers is associative, it's not for floating point numbers stored in computers. The reason is that after every operation, the calculated number is rounded back to the precision of a float or double. The order of operations therefore matters for the final result. Compilers won't reorder those operations by default. However, if you use the -ffast-math compiler flag (or /fp:fast for MSVC), then it will be able to reorder operations, even if it might change the precision of the result.

By default, compilers also compile binaries that are guaranteed to run on a wide range of processors; for example if you compile for the amd64 architecture, it will run on all amd64 machines, even old ones. If you have a newer CPU however, it might have some extra instructions that could have sped up your calculation, but would break backwards compatibility with older CPUs. If you use the -march=native compiler flag, then it will use newer instructions if possible.

Answer 3

Suspected subtle bug

Try

  printf("%a\n", 0.866025403784438);
  printf("%a\n", sqrt(3)/2);
  printf("%a\n", 0.86602540378443865);
  printf("%a\n", 0.86602540378443864676372317075294);
  printf("%a\n", nextafter(sqrt(3)/2,0));

My output

0x1.bb67ae8584ca5p-1  // 5 ULP away from the best root3/2
0x1.bb67ae8584caap-1
0x1.bb67ae8584caap-1
0x1.bb67ae8584caap-1
0x1.bb67ae8584ca9p-1

Rather than code a less than precise value, use enough significant decimal digits to form the best floating point value.
For double, use DBL_DECIMAL_DIG or more significant digits.
DBL_DECIMAL_DIG is commonly 17.

typedef struct s_float2{
    double  x;
    double  y;
}   t_float2;

// (t_float2){-0.5, 0.866025403784438};
//                 12345678901234567
(t_float2){-0.5, 0.86602540378443865};

---

OTOH, perhaps OP wanted to start further away than closest double to sqrt(3)/2?