I'm trying to write a JavaScript fractal generation algorithm from first principles. I'm aware that there are many examples out there but I wanted to incorporate additional functionality to support both Mandelbrot and 'spinning' Julia with variants such as 'Burning Ship' and 'Tricorn'. With this in mind I implemented a lightweight Complex maths library (again, I'm aware there are standard Complex js libraries out there but I wanted to build one from scratch as a learning exercise).
I tested two alternate functions, one fractal using standard maths functions and the other fractalComplex using my Complex library methods. They both work fine, but I was surprised to find that the standard version is almost twice as fast as the Complex version. I was expecting some additional overhead but not that much!
Can anyone explain why? The Complex library is using the same maths constructs 'under the covers'. Is the additional overhead purely down to object creation?
The code is reproduced below (the input parms z and c are objects of the form {re, im}).
function fractal(z, c, maxiter) {
var i, za, re, im, re2, im2;
c = (settype === JULIA ? c : z);
// Iterate until abs(z) exceeds escape radius
for (i = 0; i < maxiter; i += 1) {
if (setvar === BURNING_SHIP) {
re = Math.abs(z.re);
im = -Math.abs(z.im);
}
else if (setvar === TRICORN) {
re = z.re
im = -z.im; // conjugate z
}
else { // Mandelbrot
re = z.re;
im = z.im;
}
re2 = re * re;
im2 = im * im;
z = { // z = z² + c
re: re2 - im2 + c.re,
im: 2 * im * re + c.im
};
za = re2 + im2 // abs(z)²
if (za > 4) { // abs(z)² > radius²
break;
}
}
za = Math.sqrt(za); // abs(z)
return { i, za };
}
function fractalComplex(z, c, maxiter, n, radius) {
var i, za;
c = (settype === JULIA ? c : z);
// Iterate until abs(z) exceeds escape radius
for (i = 0; i < maxiter; i += 1) {
if (setvar === BURNING_SHIP) {
z = new Complex(Math.abs(z.re), -Math.abs(z.im))
}
if (setvar === TRICORN) {
z = z.conjugate()
}
z = z.quad(n, c); // z = zⁿ + c
za = z.abs();
if (za > radius) {
break;
}
}
return { i, za };
}
My "Complex lite" library is as follows:
// ------------------------------------------------------------------------
// A basic complex number library which implements the methods used for
// Mandelbrot and Julia Set generation.
// ------------------------------------------------------------------------
'use strict';
// Instantiate complex number object.
function Complex(re, im) {
this.re = re; // real
this.im = im; // imaginary
}
Complex.prototype = {
're': 0,
'im': 0,
// Set value.
'set': function (re, im) {
this.re = re;
this.im = im;
},
// Get magnitude.
'abs': function () {
return Math.sqrt(this.re * this.re + this.im * this.im);
},
// Get polar representation (r, θ); angle in radians.
'polar': function () {
return { r: this.abs(), θ: Math.atan2(this.im, this.re) };
},
// Get square.
'sqr': function () {
var re2 = this.re * this.re - this.im * this.im;
var im2 = 2 * this.im * this.re;
return new Complex(re2, im2);
},
// Get complex number to the real power n.
'pow': function (n) {
if (n === 0) { return new Complex(1, 0); }
if (n === 1) { return this; }
if (n === 2) { return this.sqr(); }
var pol = this.polar();
var rn = Math.pow(pol.r, n);
var θn = n * pol.θ;
return cart(rn, θn);
},
// Get conjugate.
'conjugate': function () {
return new Complex(this.re, -this.im);
},
// Get quadratic zⁿ + c.
'quad': function (n, c) {
var zn = this.pow(n);
return new Complex(zn.re + c.re, zn.im + c.im);
},
// Rotate by angle in radians.
'rotate': function (angle) {
var pol = this.polar();
angle += pol.θ;
return new Complex(pol.r * Math.cos(angle), pol.r * Math.sin(angle));
},
// String in exponent format to specified significant figures.
'toString': function (sig = 9) {
return this.re.toExponential(sig) + " + " + this.im.toExponential(sig) + "i";
},
}
// Convert polar (r, θ) to cartesian representation (re, im).
function cart(r, θ) {
var re = r * Math.cos(θ);
var im = r * Math.sin(θ);
return new Complex(re, im);
}
Additional edit 22/12/2021 11:52:
For what it's worth, this is what I eventually settled on...
function fractal(p, c, n, maxiter, radius) {
var i, za, zre, zim, tre, cre, cim, r, θ;
var lastre = 0;
var lastim = 0;
var per = 0;
if (setmode === JULIA) {
cre = c.re;
cim = c.im;
zre = p.re;
zim = p.im;
}
else { // Mandelbrot mode
cre = p.re;
cim = p.im;
zre = 0;
zim = 0;
}
// Iterate until abs(z) exceeds escape radius
for (i = 0; i < maxiter; i += 1) {
if (setvar === BURNING_SHIP) {
zre = Math.abs(zre);
zim = -Math.abs(zim);
}
else if (setvar === TRICORN) {
zim = -zim; // conjugate z
}
// z = z² + c
if (n == 2) {
tre = zre * zre - zim * zim + cre;
zim = 2 * zre * zim + cim;
zre = tre;
}
else { // z = zⁿ + c, where n is integer > 2
r = powi(Math.sqrt(zre * zre + zim * zim), n); // radiusⁿ
//r = Math.pow(Math.sqrt(zre * zre + zim * zim), n); // radiusⁿ
θ = n * Math.atan2(zim, zre); // angleⁿ
zre = r * Math.cos(θ) + cre;
zim = r * Math.sin(θ) + cim;
}
// Optimisation - periodicity check speeds
// up processing of points within set
if (PERIODCHECK) {
if (zre === lastre && zim === lastim) {
i = maxiter;
break;
}
per += 1;
if (per > 20) {
per = 0;
lastre = zre;
lastim = zim;
}
}
// ... end of optimisation
za = zre * zre + zim * zim // abs(z)²
if (za > radius) { // abs(z)² > radius²
break;
}
}
return { i, za };
}
// Optimised pow() function for integer exponents
// using 'halving and squaring'.
function powi(base, n) {
var res = 1;
while (n) {
if (n & 1) { // if n is odd
res *= base;
}
n >>= 1; // n * 2
base *= base;
}
return res;
}