# Lambdabulb

# 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
    };