Three years ago, I published a YouTube video all about creating a program capable of rendering Newton-Raphson Fractals in the form of NetBPM images. In this video, I put particular emphasis on my learning process where I started with zero knowledge of the fractal or its brilliant algorithm, and ended with a concrete creation based off my research and experimentation.
Since that video, I have been obsessed with figuring out how I can improve this program, making it more efficient and scalable. In this article, I’m going to share how I’ve optimized and extended my Newton Fractal Generator, the main addition being multithreading.
Re-Write in C
I began by taking my original Python program and re-writing it in C, re-implementing algorithms in it as follows:
Abstrating a Fractal into a Struct
In order to facilitate multithreading later in the program, it would be advantageous to start representing all of the relevant information of our fractal in a typedef struct. The essential information for any newton fractal is as follows:
- The roots
- The order (number of roots)
- The center to render from
- The amount of default iterations
- The radius of the render
- The zoom
With this in mind, I settled on this implementation:
typedef struct {
complex * roots;
int numRoots;
complex center;
int iterations;
int radius;
double zoom;
} fractal;
Newton-Raphson Root Function
I began by re-writing the star of the show, the newton() function, and adding a new iterations value:
...
// Code snippet showing the Newton algorithm
complex newton(complex x, int iterations) {
for (int i = 0; i < iterations; i++) {
if (cabs(p(x)) < H)
return x;
x -= p(x) / dydx(x);
}
}
...
You can think of the iterations value as a slider which controls the level of detail in the fractal.
As we increase this value from 0 to 7, we can visualize this unique phenomenon for the expression x^5 + x^3 + 1:
Root Finder
Outside of the basic mathematical functions (p(x), dydx(x), newton()`)
we need some systemic way of determining the roots of any given polynomial p(x).
In my previous implementation, the program cheated by already being provided with the number of roots.
This time I wanted a function that would guess up to an arbitrary number of roots (let’s say, 100) and populated an array with the given roots of our polynomial:
void determineRoots(fractal * f) {
for (int i = 0; i < ROOT_GUESSES; i++) {
complex x = i*cexp(i*I);
x = newton(x, 100);
// If a unique root is found, add it!
if (!inRoots(x *f) && isnormal(creal(x)) && isnormal(cimag(x))) {
fprintf(stderr, "Found root %d: %.2f + %.2fi\n", f->numRoots, x);
f->roots[f->numRoots] = x;
f->numRoots++;
}
}
}
Guessing occurs using the following function:
complex x = i * cexp(i*I);
The purpose of this function is to generate a sufficiently “different” complex number to perform newton root guesses on, in search of a new unique root. This specific implementation uses Euler’s formula, which we can visualize as a point advancing through spiral:
