8
\$\begingroup\$

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);
}
```
\$\endgroup\$
5
  • 5
    \$\begingroup\$ Welcome to Code Review! You'll receive better reviews if you show a complete example. For example, I recommend that you edit to show the necessary #include lines, and a main() that shows how to call your function. It can really help reviewers if they are able to compile and run your program. \$\endgroup\$ Jul 7, 2023 at 4:51
  • 7
    \$\begingroup\$ -o3 sets the output file name to 3, still no optimization. Hopefully you actually tried -O3. (Also try -O3 -march=native -ffast-math with GCC or clang, and also -fprofile-generate / run it / -fprofile-use. See also Which gcc optimization flags should I use? and this for some interesting compiler options.) \$\endgroup\$ Jul 7, 2023 at 21:44
  • 3
    \$\begingroup\$ @mehdi mirzaie "It takes 1.19s to start." ---> There is no main(). Post your code (with main() so others may attempt to compare times. \$\endgroup\$ Jul 8, 2023 at 2:46
  • \$\begingroup\$ @chux - Reinstate Monica, the link to the program is link. Here you will find more fractals that I have implemented. \$\endgroup\$ Jul 8, 2023 at 4:09
  • 5
    \$\begingroup\$ @mehdimirzaie Links too often become broken over time or code is deep is a high level link. Best to put a sample main() here. No need to add difficulties for reviewers who provide quality feedback at zero cost. \$\endgroup\$ Jul 8, 2023 at 5:24

3 Answers 3

9
\$\begingroup\$

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.

\$\endgroup\$
3
  • \$\begingroup\$ A pity there is no main() to give it a spin. \$\endgroup\$
    – greybeard
    Jul 7, 2023 at 7:03
  • \$\begingroup\$ thanks for the help, here is the program link: github.com/MehdiMirzaie2/fractol \$\endgroup\$ Jul 7, 2023 at 7:36
  • 2
    \$\begingroup\$ Interesting this code used a different constant with no explanation: 0.866025403784439, which is not much closer to sqrt(¾) than OP's 0.866025403784438. 0.86602540378443865 is closest double to sqrt(¾). \$\endgroup\$ Jul 8, 2023 at 14:19
17
\$\begingroup\$

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.

\$\endgroup\$
4
  • 4
    \$\begingroup\$ Aside: z * z * z - 1 in not as numerically stable as (z*z + z + 1)*(z - 1), especially for z near 1.0`. \$\endgroup\$ Jul 8, 2023 at 2:57
  • \$\begingroup\$ Hmmm, Need to apply the idea further: (z*z + z + 1)*(z - 1) --> (z*(z + 1) + 1)*(z - 1) to also improve z near -1. \$\endgroup\$ Jul 10, 2023 at 0:09
  • \$\begingroup\$ @chux-ReinstateMonica I would be more worried about the denominator than the numerator. How are you going to make the term \$-1 / 3 z^2\$ in the result stable near \$z = 0\$? \$\endgroup\$
    – G. Sliepen
    Jul 10, 2023 at 6:27
  • 1
    \$\begingroup\$ G. Sliepen The above suggestions help the precision for values of z near 1.0 and -1.0. Near z == 0, code relies on the exponential range of the floating point type and is not a stability problem, but an overflow one. Near z== 0, next() is stable, even if it approaches infinity. \$\endgroup\$ Jul 10, 2023 at 6:36
5
\$\begingroup\$

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?

\$\endgroup\$

Your Answer

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

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