Newton Fractal Optimization

Question

I have implemented the Newton fractal in C. I was wondering if anyone knows how I can improve the efficiency of the code without using threads or multiple processes. I have tried unrolling the while loops, but it did not do much. I also tried using -o1, -o2, and -o3 when compiling this also did not do much. If you guys have any suggestions please let me know, thanks.

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

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

    temp_deno = 3 * pow((pow(z.x, 2) + pow(z.y, 2)), 2);
    new.x = ((pow(z.x, 5) + 2 * pow(z.x, 3) * pow(z.y, 2)
                - pow(z.x, 2) + z.x * pow(z.y, 4) + pow(z.y, 2)) / temp_deno);
    new.y = ((z.y * (pow(z.x, 4) + 2 * pow(z.x, 2)
                    * pow(z.y, 2) + 2 * z.x + pow(z.y, 4))) / 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 int  flag;

    flag = 0;
    if (flag == 0)
    {
        roots[0] = (t_float2){ 1, 0 };
        roots[1] = (t_float2){ -0.5, sqrt(3) / 2 };
        roots[2] = (t_float2){ -0.5, -sqrt(3) / 2 };
        flag = 1;
    }
    return ;
}

void    calculate_newton(t_fractal *fractal)
{
    t_float2        z;
    int             iterations;
    int             i;
    t_float2        difference;
    static t_float2 roots[3];

    fractal->name = "newton";
    init_roots(&roots);
    iterations = -1;
    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);
        i = -1;
        while (++i < 3)
        {
            difference = find_diff(z, roots[i]);
            if (fabs(difference.x) < fractal->tolerance
                && fabs(difference.y) < fractal->tolerance)
                return (put_color_to_pixel(fractal, fractal->x, fractal->y,
                        fractal->color * iterations / 2));
        }
    }
    put_color_to_pixel(fractal, fractal->x, fractal->y, 0x000000);
}

Would it be better to manually do the math operations, rather than using the <math.h> library.

Answer 1

My suggestion:

A key factor in the total cost is the number of iterations performed from each starting point. What about memoizing these ? I mean, for every pixel processed, keep the number of iterations. But when iterating, when you land on a pixel already processed, assume the number of remaining iterations to be the same.

This is an approximation, but you assume anyway a single color per pixel.

Answer 2

You should start by computing the powers with multiplications, not using the pow function:

static t_float2 next(t_float2 z)
{
    t_float2    new;
    float       temp_deno;
    float       x2 = z.x * z.x;
    float       y2 = z.y * z.y;
    float       x4 = zx2 * zx2;
    float       y4 = zy2 * zy2;

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

Then note that you can simplify the algebra by factoring (x2 + y2) * (x2 + y2) which is also x4 + 2 * x2 * y2 + y4:

static t_float2 next(t_float2 z)
{
    float x2 = z.x * z.x;
    float y2 = z.y * z.y;
    float temp_deno = 3 * (x2 + y2) * (x2 + y2);
    return (t_float2){
        z.x * 2 / 3 - (y2 - x2) / temp_deno,
        z.y * 2 / 3 - (2 * z.x * z.y) / temp_deno
    };
}

roots should be precomputed at compile time. In C you should just use constants:

    #define COS_PI_6  0.866025403784439  // sqrt(3) / 2
    static t_float2 roots[3] = {{ 1, 0 }, { -0.5, COS_PI_6 }, { -0.5, -COS_PI_6 }};

Answer 3

yes using pow for small integer exponents is a bad idea (as its relatively complex operation involving log,exp,*,/) computing it simply with * is much better. On top of this you use multiple powers of the same variable and compute them over again you can reuse lower powers instead for example change:

x3 = pow(x,3);
x5 = pow(x,5);

into:

x3 = x*x*x;
x5 = x3*x*x;

this will lower the number of operations.

On top of this you compute everything over and over for each pixel ...

I assume your code looks like this:

t_fractal fractal;
for (fractal.y=0;fractal.y<ys;fractal.y++)
 for (fractal.x=0;fractal.x<xs;fractal.x++)
  calculate_newton(&fractal);

while inside calculate_newton(t_fractal *fractal) you recalculate all the stuff related to x,y position over and over... if you have y loop the outer one you can compute all the stuff just once per each y change instead of doing it xs times more ... for that you should move those variables into fractal context/struct

in some cases you can avoid using * for powers of incrementing values (however depends on used computing architecture if it has any meaning from performance aspect) see

• How to get a square root for 32 bit input in one clock cycle only?

for details on how to do it.

On top of that if you continuously zooming/panning you can reuse already computed values from previous frame like I did in here for Mandelbrot

However without parallelism that is about all you can do...

With parallelism unless you got many cores CPU/computer it will not boost the speed much its much better to use GPU for this see:

• GLSL RT Mandelbrot

Where you render single QUAD/rectangle covering your view and compute the fractal inside fragment shader (which for GLSL is more or less C like language so it does not differ much from what you already have).