# # Lambdabulb

A 3D fractal formula related to the Mandelbulb.

History/credits:

- Fractal Chats Discord discussion August 2024
- original definition at fractalforums.com/mandelbulb-renderings/lambdabulb,
- modified by machina-infinitum
- fragmentarium version assisted by dark-beam
- translated to lycium et al’s FractalTracer by bezo97
- polynomial optimization by me

## # 1 Formula

\[ z \to c (z - z^p)\]

### # 1.1 Trigonometric

The Lambdabulb formula is specified using spherical trigonmetry with a “triplex” vector multiplication and power similar to the Mandelbulb. Pseudo-code:

```
formula(z, c):
z = triplexMult(c, z - triplexPow(z, power, phase));
triplexPow(z, power, phase):
r = length(z);
theta = atan2(z.y, z.x);
phi = acos(z.z / r);
r = pow(r, power)
theta = theta * power;
phi = phi * power + phase;
return
{ r * sin(phi) * cos(theta)
, r * sin(phi) * sin(theta)
, r * cos(phi)
};
triplexMul(u, v):
ru = length(u);
rv = length(v);
a = 1 - (u.z * v.z) / (ru * rv);
return
{ a * (u.x * v.x - u.y * v.y)
, a * (v.x * u.y + u.x * v.y)
, rv * u.z + ru * v.z
};
```

Image rendering can be done via ray marching with automatic differentiation for distance estimates, for example code see github.com/lycium/FractalTracer.

### # 1.2 Polynomial (Power 4)

Profiling with `valgrind --tool=cachegrind`

showed that
50% of the program runtime was in `libm`

for the trigonmetry and powers.

Similar to Chebyshev polynomials,
`trig(f(inverse trig))`

can often be replaced by algebra,
which is usually less computationally expensive
and sometimes more accurate.

I used wxMaxima computer algebra system to work out the maths for power 4 triplexPow. A similar approach should work for other integer powers.

```
(%i1) power : 4;
(%o1) 4
(%i2) factor(trigexpand(cos(power*atan2(y,x))));
(%o2) ((y^2-2*x*y-x^2)*(y^2+2*x*y-x^2))/(y^2+x^2)^2
(%i3) factor(trigexpand(sin(power*atan2(y,x))));
(%o3) -(4*x*y*(y-x)*(y+x))/(y^2+x^2)^2
(%i4) at(factor(expand(trigexpand(factor(trigexpand(cos(phase + power * acos(z / sqrt(x^2+y^2+z^2)))))))), [cos(phase) = c, sin(phase) = s]);
(%o4) (c*z^4-4*s*sqrt(y^2+x^2)*z^3-6*c*y^2*z^2-6*c*x^2*z^2+4*s*y^2*sqrt(y^2+x^2)*z+4*s*x^2*sqrt(y^2+x^2)*z+c*y^4+2*c*x^2*y^2+c*x^4)/(z^2+y^2+x^2)^2
(%i5) at(factor(expand(trigexpand(factor(trigexpand(sin(phase + power * acos(z / sqrt(x^2+y^2+z^2)))))))), [cos(phase) = c, sin(phase) = s]);
(%o5) (s*z^4+4*c*sqrt(y^2+x^2)*z^3-6*s*y^2*z^2-6*s*x^2*z^2-4*c*y^2*sqrt(y^2+x^2)*z-4*c*x^2*sqrt(y^2+x^2)*z+s*y^4+2*s*x^2*y^2+s*x^4)/(z^2+y^2+x^2)^2
```

With some manual common subexpression elimination,
and noting that dividing `cos_phi`

and `sin_phi`

by \(r^4\) is unnecessary
as the end result is multiplied by \(r^4\),
this translates into pseudo-code like so:

```
triplexPow4(z, cos_phase, sin_phase):
x2 = z.x * z.x;
y2 = z.y * z.y;
z2 = z.z * z.z;
xy2 = 2 * z.x * z.y;
y2mx2 = y2 - x2;
y2px2 = y2 + x2;
y2px22 = y2px2 * y2px2;
cos_theta = (y2mx2 - xy2) * (y2mx2 + xy2) / y2px22;
sin_theta = -2 * xy2 * y2mx2 / y2px22;
z22 = z2 * z2;
a = z22 + y2px22 - 6 * y2px2 * z2;
b = sqrt(y2px2) * z.z * (z2 - y2px2);
cos_phi = cos_phase * a - 4 * sin_phase * b;
sin_phi = sin_phase * a + 4 * cos_phase * b;
return
{ sin_phi * cos_theta,
, sin_phi * sin_theta,
, cos_phi
};
```

## # 2 Benchmarks

FractalTracer compiled with g++ 12 on AMD Ryzen 2700X CPU (16 threads) running Debian Bookworm, at 1280x720 resolution. Lower timings are better.

- trigonometric, double precision, double epsilon: 126 seconds per pass
- trigonometric, double precision, single epsilon: 58 seconds per pass
- trigonometric, single precision, double epsilon: 73 seconds per pass
- trigonometric, single precision, single epsilon: 36 seconds per pass
- polynomial, double precision, double epsilon: 78 seconds per pass
- polynomial, double precision, single epsilon: 40 seocnds per pass
- polynomial, single precision, double epsilon: 42 seconds per pass
- polynomial, single precision, single epsilon: 22 seconds per pass

## # 3 Images

Comparison of trigonometric (4 samples per pixel) vs polynomial (32 samples per pixel), double vs single floating point precision, double vs single tuned ray marching epsilons.

Using single-tuned epsilons (which are larger) has the effect of thickening the fractal shape.

Polynomial seems improved over trigonometric, the fractal is much more solid. Not sure whether there is a bug somewhere or if it’s just numerical precision issues.

Note the different appearance near the horizon of the penultimate polynomial-single-double image, this is probably due to self-intersections caused by to too-low epsilons for the precision.

Parameters:

```
power = 4
phase = 1.815142
c = { 1.035, -0.317, 0.013 }
```

The fixed c makes it a “Juliabulb” variant.

- trigonometric, double precision, double epsilon

- trigonometric, double precision, single epsilon

- trigonometric, single precision, double epsilon

- trigonometric, single precision, single epsilon

- polynomial, double precision, double epsilon

- polynomial, double precision, single epsilon

- polynomial, single precision, double epsilon

- polynomial, single precision, single epsilon