| Contents
 
 | Hypercomplex Fractals Hypercomplex numbers are similar to the usual 2D complex numbers, except they can be extended to 3 dimensions or more. When you use hypercomplex numbers to generate fractals, you can create some interesting looking 3D fractals.
 
 Quadratic 3D Mandelbulb Set - C++, 7/3/09
 
 |  This hypercomplex fractal is based on Daniel White's creative formula for squaring a 3D hypercomplex number ("triplex") by applying two consecutive rotations using a spherical coordinate system: {x,y,z}2 = r2{cos(θ)cos(φ),sin(θ)cos(φ),-sin(φ)}
 r2=x2+y2+z2, θ=2atan2(y,x), φ=2sin-1(z/r)
 This formula can also be expressed in non-trigonometric form:
 
  One drawback about this formula is that it doesn't allow you to define a very consistent set of operators (click here for a more consistent formula). However, I think it creates the best looking quadratic 3D Mandelbrot set so far. Here is some Mathematica code demonstrating a simple way to render this 3D fractal as a depth map by slowly marching tiny cubes (voxels) toward the boundary:
 (* runtime: 1.7 minutes, increase n for higher resolution *)
 n = 100; norm[x_] := x.x; TriplexSquare[{x_, y_, z_}] := If[x == y == 0, {-z^2, 0, 0}, Module[{a = 1 - z^2/(x^2 + y^2)}, {(x^2 - y^2)a, 2 x y a, -2 z Sqrt[x^2 + y^2]}]];
 Mandelbulb[c_] := Module[{p = {0,0, 0}, i = 0}, While[i < 24 && norm[p] < 4, p = TriplexSquare[p] + c; i++]; i];
 image = Table[z = 1.5; While[z >= -0.1 && Mandelbulb[{x, y, z}] < 24, z -= 3/n]; z, {y, -1.5, 1.5, 3/n}, {x, -2, 1, 3/n}];
 ListDensityPlot[image, Mesh -> False, Frame -> False, PlotRange -> {-0.1, 1.5}]
 
 Here is some inefficient code demonstrating how you can render this isosurface in POV-Ray using recursive functions:
 // runtime: 4.5 minutes
 camera{location <-2.5,5,5> look_at <-0.5,0,0.25> up z sky z angle 25}
 light_source{20*z,1}
 #declare f=function(i,x,y,z,xc,yc,zc) {select(i>0 & x*x+y*y+z*z<4, 0, sqrt(x*x+y*y+z*z), f(i-1,(x*x-y*y)*(1-z*z/(x*x+y*y))+xc,2*x*y*(1-z*z/(x*x+y*y))+yc,-2*z*sqrt(x*x+y*y)+zc,xc,yc,zc))};
 isosurface{function{f(24,x,y,z,x,y,z)} threshold 2 max_gradient 10 contained_by{sphere{<-0.5,0,0>,2}} pigment{rgb 1}}
 
 When I first set out to render these fractals, I couldn't find any software that could ray-trace iterated isosurfaces so I had to write my own C++ program for it. At first, this seemed like a daunting task, but with some help from Thomas Ludwig, and some useful links, it turned out not to be as difficult as I thought. I made the top image using a volumetric technique described by Krzysztof Marczak. This image was precomputed using a volume of 1.4 billion voxels. I made the other images using James Kajiya's path tracing method for Global Illumination (GI). This method uses Monte Carlo integration to solve the rendering equation.
 
 Links
 Mandelbrot 3D, cut away view - by Krzysztof Marczak
 metallic rendering - by Thomas Ludwig
 3D Mandelbrot, Coral Rift - rendering with global illumination, by Thomas Ludwig, see also the discussion on Fractal Forums
 inside view - by Krzysztof Marczak
 animation - by Bib
 animations - uses different powers for each rotation, by Bill Rood
 
 Global Illumination and Participating Media Links
 smallpt - Global Illumination in 99 lines of C++, by Kevin Beason
 MiniLight - simple global illumination renderer
 Subsurface Scattering - Wikipedia
 Subsurface Scattering - papers by Henrik Jensen, explains how to simulate subsurface scattering using photon mapping
 
 | 
 3D Power Mandelbulb Set - C++, 8/20/09
 
 |    Here is my first rendering of an 8th order Mandelbulb set, based on the following generalized variation of Daniel White's original squarring formula: {x,y,z}n = rn{cos(θ)cos(φ),sin(θ)cos(φ),-sin(φ)}
 r=sqrt(x2+y2+z2), θ=n atan2(y,x), φ=n sin-1(z/r)
 I originally posted this image on 8/30/09. Daniel White was very excited when he first saw this rendering because it has significantly finer 3D details than the quadratic version. He wrote a beautiful article about it, calling it the "Mandelbulb" (I believe it was Rudy Rucker who coined the name). The Mandelbulb rapidly grew in popularity, and it was presented in many popular articles including Slashdot, New Scientist magazine, Wikipedia, and Science & Vie magazine. Here is some Mathematica code:
 (* runtime: 1.4 minutes, increase n for higher resolution *)
 n = 100; norm[x_] := x.x; TriplexPow[{x_, y_, z_}, n_] := If[x == y == 0.0, 0.0, Module[{r = Sqrt[x^2 + y^2 + z^2], theta = n ArcTan[x, y], phi}, phi = n ArcSin[z/r]; r^n{Cos[theta]Cos[phi], Sin[theta]Cos[phi], -Sin[phi]}]];
 Mandelbulb[c_] := Module[{p = {0, 0, 0}, i = 0}, While[i < 24 && norm[p] < 4, p = TriplexPow[p, 8] + c; i++]; i];
 image = Table[z = 1.1; While[z >= -0.1 && Mandelbulb[{x, y, z}] < 24, z -= 2.2/n]; z, {y, -1.1, 1.1, 2.2/n}, {x, -1.1, 1.1, 2.2/n}];
 ListDensityPlot[image, Mesh -> False, Frame -> False, PlotRange -> {-0.1, 1.1}]
 
 We can increase the speed / accuracy and emphasize the boundaries using distance estimation:
 (* runtime: 1.7 minute *)
 abs[x_] := Sqrt[x.x]; TriplexMult[{x1_, y1_, z1_}, {x2_, y2_, z2_}] := Module[{r1 = Sqrt[x1^2 + y1^2], r2 = Sqrt[x2^2 + y2^2], a}, a = 1 - z1 z2/(r1 r2); {a(x1 x2 - y1 y2), a(x2 y1 + x1 y2), -(r2 z1 + r1 z2)}];
 DistanceEstimate[c_, n_] :=Module[{p = {0, 0, 0}, dp = {1, 0, 0}, r, dr = 1, i = 0, loop = True}, While[loop, p = TriplexPow[p, n] + c; r = abs[p]; i++; loop = (i < 24 && r < 2); If[loop, dp = n TriplexMult[dp, TriplexPow[p, n - 1]] + {1, 0, 0}; dr = abs[dp]]]; 2 r Log[r]/dr/4];
 n = 100; image = Table[z = d = 1.1; While[d > 0.001 && z >= -0.1, d = DistanceEstimate[{x, y, z}, 8]; z -=d]; z, {y, -1.1, 1.1, 2.2/n}, {x, -1.1, 1.1, 2.2/n}];
 ListDensityPlot[image,Mesh -> False, Frame -> False, PlotRange -> {-0.1, 1.1}]
 
 Here is some inefficient code demonstrating how you can render this isosurface in POV-Ray using recursive functions:
 // runtime: 28 minutes
 camera{location <1.75,4,5.5> look_at 0 up z sky z angle 25}
 light_source{20*z,1}
 #declare n=8;
 #declare f=function(i,x,y,z,xc,yc,zc) {select(i>0 & x*x+y*y+z*z<4, 0, sqrt(x*x+y*y+z*z),f(i-1,pow(x*x+y*y+z*z,n/2)*cos(n*atan2(y,x))*cos(n*atan2(z,sqrt(x*x+y*y)))+xc,pow(x*x+y*y+z*z,n/2)*sin(n*atan2(y,x))*cos(n*atan2(z,sqrt(x*x+y*y)))+yc,-pow(x*x+y*y+z*z,n/2)*sin(n*atan2(z,sqrt(x*x+y*y)))+zc,xc,yc,zc))};
 isosurface{function{f(24,0,0,0,x,y,z)} threshold 2 max_gradient 1000 contained_by{sphere{0,1.5}} pigment{rgb 1}}
 
 I posted some non-trigonometric expansions for this formula here and here. The above animation was created by adding a phase shift to φ.
 
 | 
 | Links Mystery of the Real 3D Mandelbrot Fractal - beautiful article by Daniel White
 Mandelbulber - C++ program by Krzysztof Marczak
 Mandel-Maya Madness! - Maya article by Duncan Brinsmead
 Mandelbulb - French article in Images of Mathematics, by Jos Leys
 Mandelbulb Flight - "hot" glowing animation, by Krzysztof Marczak, see also Minibulbs, Flight through BulbBox, Dangerous Spaceship
 3D Mandelbulb crater transformation - animation by continuously varying the power, by Bib
 The Lava Dome - Mandelbulb flight, by Snakehand
 Mandelbulb 2 - animation by Iņigo Quilez
 Siebenfach, remix - 8th order 3D Mandelbrot set (has 7-fold symmetry), by Thomas Ludwig
 7th order 3D Mandelbrot set - rendered from a volume of 27 billion voxels, by Krzysztof Marczak
 Mandelbulb Spirals - cross-section zoom, by Krzysztof Marczak
 In Search of the Mandelbulb - article by Rudy Rucker
 Mandelbulb gallery - by Jos Leys
 Mandelier - negative power 3D Mandelbrot set, by Christian Buchner
 6th order animation - by Bib
 7th order Mandelbrot set - nice metallic look, by Christopher Kulla
 8th order Mandelbrot set rotation animation - by Iņigo Quilez
 fake ambient occlusion using orbit traps - rendered in 2 seconds using the Graphics Processing Unit (GPU), by Iņigo Quilez
 Sliced Mandelbulb - spectacular deep zoom renderings
 3D Mandelbulb power 8 growth - iteration depth animation
 realtime animation continuously varying the power - by Christian Buchner
 Mandelbulb Radiolaria, Ring Structure, The Formula - interesting variations, by Tom Beddard
 Into the heart of the Mandelbulb, phase animation, zoom - animations by Daniel White
 Exploding Chicken - degree pi Mandelbulb, by Kraftwerk
 Lost in the Mandelbulb - Mandelbulb with lights, by Daniel White
 Monopolar Coordinate System - proposed by Michael, see these renderings by Jesse
 3D printed Mandelbulb - it's not actually the first, but one of the better ones I've seen, by Gary Quinn
 
 | 
 Positive Z-Component Variation of the 3D Mandelbulb Set - C++, 10/3/09
 
 |  One strange thing about the previous formula is that it gives {x, y, z}1 = {x, y, -z}. Here is a more consistent formula without the negative sign in front of the z-component: {x,y,z}n = rn{cos(θ)cos(φ),sin(θ)cos(φ),sin(φ)}, r=sqrt(x2+y2+z2), θ=n atan2(y,x), φ=n sin-1(z/r)
 
 Using this formula, we have been able to come up with definitions for multiplication, division, powers, roots, exponentials, logarithms, and trigonometric functions. It should be noted that the multiplication operator is commutative but not associative. I also posted some non-trigonometric expansions for the power formula here. One drawback about this formula is that it gives a flat looking quadratic Mandelbrot set. Here is the non-trigonometric form of the multiplication operator:
 
  Here is the non-trigonometric form of the division operator:
 
  
 | 
 3D Juliabulb Set - C++, 7/28/09
 
 Lambdabulb 3D Juliabulb Set - C++, 11/25/09
 
 |  An endless variety of 3D Mandelbulb and Juliabulb sets are possible using triplex numbers. For example, here is a "Lambdabulb" fractal. This is actually a Juliabulb set that is based on a unique formula suggested by Ron Barnett. Here is some Mathematica code: (* runtime: 2 minutes, increase n for higher resolution *)
 TriplexPow[{x_, y_, z_}, n_] := Module[{r = Sqrt[x^2 + y^2 + z^2], theta = n ArcTan[x, y], phi}, phi = n ArcSin[z/r]; r^n{Cos[theta]Cos[phi], Sin[theta]Cos[phi], Sin[phi]}];
 TriplexMult[{x1_, y1_, z1_}, {x2_, y2_, z2_}] := Module[{r1 = Sqrt[x1^2 + y1^2], r2 = Sqrt[x2^2 + y2^2], a}, a = 1 -z1 z2/(r1 r2); {a(x1 x2 - y1 y2), a(x2 y1 + x1 y2), r2 z1 + r1 z2}];
 c = {1.0625, 0.2375, 0.0}; n = 100; norm[x_] := x.x;
 Lambdabulb[c_] := Module[{p = c, i = 0}, While[i < 40 && norm[p] < 4, p = TriplexMult[c, p - TriplexPow[p, 4]]; i++]; i];
 image = Table[z = 1.2; While[z >= 0 && Lambdabulb[{x, y, z}] < 40, z -= 2.4/n]; z, {y, -1.2, 1.2, 2.4/n}, {x, -1.2, 1.2, 2.4/n}];
 ListDensityPlot[image, Mesh -> False, Frame -> False, PlotRange -> {0, 1.2}]
 
 | 
 3D Mandelbrot / Mandelbar Set - C++, 7/5/09
 
 |    {x,y,z}2 = {x2-y2-z2, 2xy, -2xz} This method combines the Mandelbrot set with the Mandelbar set (Tricorn). I think Eric Baird was the first to suggest this formula. The right picture is a zoom of the first miniature Mandelbrot set. I think it looks like a space ship with some sort of radar gun. Unfortunately, Julia sets created using this method look the same as quaternion Julia sets. Here is some Mathematica code:
 (* runtime: 3.2 minutes, increase n for higher resolution *)
 n = 100; norm[x_] := x.x; TricornSquare[{x_, y_, z_}] := {x^2 - y^2 - z^2, 2x y, -2x z};
 Mandelbrot[c_] := Module[{p = {0, 0, 0}, i = 0}, While[i < 20 && norm[p] < 4, p = TricornSquare[p] + c; i++]; i];
 image = Table[y = 0.04; While[y > -0.02 && Mandelbrot[{x, y, z}] < 20, y -= 0.08/n]; y, {z, -0.04, 0.04, 0.08/n}, {x, -1.8, -1.72, 0.08/n}];
 ListDensityPlot[image, Mesh -> False, Frame -> False, PlotRange -> {-0.02, 0.04}]
 
 Here is some inefficient code demonstrating how you can render this in POV-Ray:
 // runtime: 4 minutes
 camera{location <-3.5,5,4> look_at <-0.5,0,0.4> up z sky z angle 25}
 light_source{20*z,1}
 #declare f=function(i,x,y,z,xc,yc,zc) {select(i>0 & x*x+y*y+z*z<4, 0, sqrt(x*x+y*y+z*z), f(i-1,x*x-y*y-z*z+xc,2*x*y+yc,-2*x*z+zc,xc,yc,zc))};
 isosurface{function{f(24,x,y,z,x,y,z)} threshold 2 max_gradient 25 contained_by{sphere{<-0.5,0,0>,2}} pigment{rgb 1}}
 
 | 
 3D Mandelbrot Set Pickover Stalks - C++, 8/18/09
 
 |    Pickover stalks are 2D fractals, composed of cross-shaped orbit traps invented by Clifford Pickover. These images were created by extending this concept into 3 dimensions. The right image has helical orbit traps wrapped around the Pickover stalks for added effect. See also my 2D Pickover stalks fractal. Here is some Mathematica code: (* runtime: 1.5 minutes, increase n for higher resolution *)
 n = 100; norm[x_] := x.x; TricornSquare[{x_, y_, z_}] := {x^2 - y^2 - z^2, 2x y, -2x z};
 OrbitTrap[{x_, y_, z_}] := 1.0 - Min[x^2 + y^2, x^2 + z^2, y^2 + z^2]/0.3^2;
 Mandelbrot[c_] := Module[{p = {0, 0, 0},i = 0}, While[i < 100 && norm[p] < 4 && (i < 2 || OrbitTrap[p] < 0.0), p = TricornSquare[p] + c; i++]; OrbitTrap[p]];
 image = Table[z = 1.5; While[z > -1.5 && Mandelbrot[{x, y, z}] < 0.0, z -= 3/n]; z, {y, -1.5, 1.5, 3/n}, {x, -2.0, 1.0, 3/n}];
 ListDensityPlot[image, Mesh -> False, Frame -> False, PlotRange -> All]
 
 Links
 3D Pickover Stalks - by David Makin
 3D Orbit Traps - by Iņigo Quilez
 
 | 
 |      Here are some Pickover stalks using other 3D Mandelbrot formulas. 
 | 
 4D Quaternion Mandelbrot ("Mandelquat") Set - C++, 7/7/09
 Quaternion Mandelbrot Set Orbits - POV-Ray 3.1, 1/6/10
 
 |  Here are the orbits of the quaternion Mandelbrot set. This image was created by simply revolving the 2D Mandelbrot orbits. 
 | 
 4D Quaternion Julia Set - Mathematica 4.2: 7/9/04, POV-Ray 3.1: 8/5/04
 
 |    This used to be the most popular method for creating hypercomplex fractals. Here is some Mathematica code: (* runtime: 1.5 minutes *)
 n = 100; dx = 3/n; c = {-0.2, 0.8, 0, 0}; norm[q_] := q.q;
 QuatSquare[{x_, y_, z_, w_}] := {x^2 - y^2 - z^2 - w^2, 2 x y, 2 x z, 2 x w};
 Julia[q0_] := Module[{q = q0, i = 0}, While[i < 12 && norm[q] < 4, q = QuatSquare[q] + c; i++]; i];
 image = Table[z = 1.5; While[z > -1.5 && Julia[{x, y, z, 0}] < 12, z -= 3/n]; z, {y, -1.5, 1.5, dx}, {x, -1.5, 1.5, dx}];
 ListDensityPlot[image, Mesh -> False, Frame -> False]
 
 We can use fast Phong shading to make it look 3D:
 (* runtime: 0.3 second *)
 Normalize[x_] := x/Sqrt[x.x];
 ListDensityPlot[Table[x = (j - 1) dx - 1.5; y = (i - 1) dx - 1.5; z = image[[i, j]]; normal = Normalize[{(image[[i, j + 1]] - image[[i, j - 1]])z, (image[[i + 1, j]] - image[[i - 1, j]]) z, 2dx}]; L = {1, 1, 0} - {x, y, z}; reflect = Normalize[2(normal.L)normal - L]; -(normal.Normalize[L] + reflect[[3]]), {i, 2, n}, {j, 2, n}],Mesh -> False, Frame -> False]
 
 POV-Ray has an efficient algorithm for ray tracing quaternion Julia sets:
 // runtime: 8 seconds
 camera{location <1,1,-2> look_at 0}
 light_source{<1.5,1.5,-1.5>,1}
 julia_fractal{<-0.2,0.8,0,0> quaternion sqr max_iteration 12 precision 400 pigment{rgb 1}}
 
 Links
 Ray Tracing Quaternion Julia Sets on the GPU - by Keenan Crane, uses Graphics Processing Unit (GPU) for fast calculations, includes source code, see also his explanation
 Hypercomplex Iterations - nice pictures
 Quaternion Julia Sets - multiple isosurfaces at varying iteration depths, by Gert Buschmann
 YouTube Video - by John Hart
 real time ray tracer - ported from Keenan Crane's program by Tom Beddard
 
 | 
 Inverse 4D Quaternion Julia Set - C++ version: 9/15/09; POV-Ray, version: 7/4/06; Mathematica version: 7/18/04
 
 |        Here is the 4D quaternion Julia set using the inverse Julia set method. The 4th dimension has been color-coded. This method is much faster, but some regions are faint because they attract much slower. The left image contains 2 billion points and the right image contains 500 million tiny spheres. See also my POV-Ray code and rotatable 3D version. Here is some Mathematica code: (* runtime: 5 seconds *)
 Sqrt2[q_] := Module[{r = Sqrt[Plus @@ Map[#^2 &, q]], a, b}, a = Sqrt[(q[[1]] + r)/2]; b = (r - q[[1]]) a/(q[[2]]^2 + q[[3]]^2 + q[[4]]^2); {a, b q[[2]], b q[[3]], b q[[4]]}];
 QInverse[qlist_] := Flatten[Map[Module[{q = Sqrt2[# - c]}, {q, -q}] &, qlist], 1];
 c = {-0.2, 0.8, 0, 0}; qlist = {{1.0, 1.0, 1.0, 1.0}};
 Do[qlist = QInverse[qlist], {12}];
 Show[Graphics3D[{PointSize[0.005], {Hue[#[[4]]], Point[Delete[#, 4]]} & /@ qlist}, Boxed -> False, Background -> RGBColor[0, 0, 0]]]
 
 Links
 QJulia - inverse quaternion Julia set program by Chris Laurel
 Quaternion Julia Crystal - laser-etched cube by Bathsheba Grossman
 
 | 
 Mandelbox - 2/22/11
 
 |  The Mandelbox fractal is based on an unusual formula originally proposed by Tom Lowe, inspired from the Mandelbulb. Technically, I don't think you can call this a "hypercomplex" fractal, but it seems to exhibit a greater variety of character than the Mandelbulb and has rapidly grown in popularity among fractal artists. Here is some Mathematica code: (* runtime: 3.5 minutes *)
 n = 100; scale = 2; r0 = 0.5; f = 1; norm[x_] := x.x;
 BoxFold[p_] := Map[If[# > 1, 2 - #, If[# < -1, -2 - #, #]] &, p];
 BallFold[r0_, p_] := Module[{r = Sqrt[norm[p]]}, p/If[r < r0,r0, If[r < 1, r, 1]]^2];
 MandelBox[c_] := Module[{p = {0, 0, 0}, i = 0}, While[i < 24 && norm[p] < 4,p = scale BallFold[r0, f BoxFold[p]] + c; i++]; i];
 image = Table[z = 1.0; While[z >= 0 && MandelBox[{x, y, z}] < 3,z -= 1/n]; z, {y, -1, 1, 2/n}, {x, -1, 1, 2/n}];
 ListDensityPlot[image, Mesh -> False, Frame -> False, PlotRange -> {-0.1, 1}]
 
 Links
 Mandelbox Trip - extremely popular animation by Krzysztof Marczak
 Happy Birthday Mandelbox! - showcase of works by various artists on Fractal Forums
 Mandelbox - French article in Images of Mathematics, by Jos Leys
 Could Be Trouble - impressive images by Hal Tenny
 Fractal Wasp Troll - impressive image by Johan Andersson
 Boxplorer - simple program to zoom into the Mandelbox in real time, by Andy Deck
 The Last Tree - an interesting example of a lone tree, by stoni
 Untitled - by Florian Sanwald
 
 | 
 Inverse 3D Juliabulb Set - C++, 9/21/09
 
 |        These 3D Juliabulb sets were calculated using the inverse Julia set method. In order to implement this method, I had to find the square root of Daniel White's squarring formula. The left image contains 2 billion tiny spheres and the right image contains points. This square root formula has 4 roots, so the total number of points increases by a factor of 4 after each iteration. The following Mathematica code is written for ease of understanding and is not intended to be efficient: (* runtime: 3 seconds *)
 TriplexSqrt[{x_, y_, z_}, sign1_, sign2_] := Module[{r = Sqrt[x^2 + y^2 + z^2], a, b, c,d}, a = sign1 Sqrt[y^2 + x^2 (2 + z^2/(x^2 + y^2)) - 2 x  r]; b = sign2 Sqrt[r - x - a]/2; c = y^4 + x^2(y^2 - z^2); d = a (x (r + x) + y^2); {b ((x y^2 - d)(x^2 + y^2) - x^3 z^2)/(y c), b, -z/Sqrt[2(r - d(x^2 + y^2)/c)]}];
 c = {-0.2, 0.8, 0}; plist = {{1.0, 1.0, 1.0}, {-1.0, -1.0, -1.0}};
 Do[plist = Flatten[Map[{TriplexSqrt[# - c, 1, 1], TriplexSqrt[# - c, -1, 1], TriplexSqrt[# - c, 1, -1], TriplexSqrt[# - c, -1, -1]} &, plist], 1], {6}];
 Show[Graphics3D[{PointSize[0.005], Point /@ plist}, Boxed -> False, ViewPoint -> {0, 10, 0}]]
 
 | 
 |  Here is a 4th order inverse Juliabulb set with 600 million spheres. The triplex raised to the 4th power has 16 unique valid roots. Garth Thornton was the first person to point out that, in general, the triplex raised to the nth power has n2 unique valid roots. I posted a formula for all integer powers here. 
 There are a variety of Modified Inverse Iteration Methods (MIIM) which can be used to speed up the calculation by pruning dense branches of the tree. One such method is to prune branches when the map becomes contractive (the cumulative derivative becomes large). Here is some Mathematica code:
 (* runtime: 3 seconds *)
 norm[x_] := x.x;
 TriplexPow[{x_, y_, z_}, n_] := Module[{r = Sqrt[x^2 + y^2 + z^2], theta = n ArcTan[x, y], phi},phi = n ArcSin[z/r]; r^n{Cos[theta]Cos[phi], Sin[theta]Cos[phi], Sin[phi]}];
 TriplexRoot[{x_, y_, z_}, n_, ktheta_, kphi_] := Module[{dk = If[Mod[n, 2] == 0 && Abs[n] < 4kphi + If[z < 0, 0, 1] <= 3Abs[n], Sign[z]If[Mod[n, 4] == 0, -1, 1], 0], r = Sqrt[x^2 + y^2 + z^2],theta, phi}, theta = (ArcTan[x, y] + (2ktheta + dk)Pi)/n; phi = (ArcSin[z/r] + (2kphi - dk)Pi)/n; r^(1/n){Cos[theta] Cos[phi], Sin[theta]Cos[phi], Sin[phi]}];
 c = {-0.2, 0.8, 0}; plist = {}; dpmax = 1.0*^6; imax = 100; p = Table[{1.0, 1.0, 1.0}, {imax}]; dp = Table[1, {imax}]; roots = Table[0, {imax}]; i = 2; power = 4; nroot = power^2;
 While[i > 1 && Length[plist] < 5000, p[[i]] = TriplexRoot[p[[i - 1]] - c, power, Mod[roots[[i]], power], Floor[roots[[i]]/power]]; dp[[i]] = Abs[power]Sqrt[norm[TriplexPow[p[[i]], power - 1]]]dp[[i - 1]]; plist = Append[plist, p[[i]]]; prune = (i == imax || dp[[i]] > dpmax); If[prune, While[i > 1 && roots[[i]] == nroot - 1, roots[[i]] = 0; i--]; roots[[i]]++, i++; roots[[i]] = 0]];
 Show[Graphics3D[Point /@ plist]]
 
 | 
 3D Glynn Juliabulb Set - C++, 11/17/09
 
 |  Here is a 3D variation of the Glynn Julia set, based on my modified version of Daniel White's squarring formula for raising a 3D hypercomplex number ("triplex") to the power of 1.5. This was rendered as a point cloud with many tiny spheres using the Modified Inverse Iteration Method (MIIM). Finding all the unique valid roots was a challenge. In general, I found that there are at least 2 roots when x < 0 and at least one root when x > 0. But there is also a strange conic section region containing 3 roots when z2 > x2 + y2. It should also be noted that convergence is very sensitive to the initial seed value. I found the best results when the seed value is inside one of the cones of the conic section. But this only finds one half of the overall fractal. You can find the other half by placing the seed in the other cone, but I preferred to leave this part out, because it gives a nice cross-section so you can see inside the fractal. See also my 2D Glynn Julia sets. Here is some Mathematica code: (* runtime: 15 seconds *)
 power = 1.5; norm[x_] := x.x;
 TriplexPow[{x_, y_, z_}, n_] := Module[{r = Sqrt[x^2 + y^2 + z^2], theta = n ArcTan[x, y], phi},phi = n ArcSin[z/r]; r^n{Cos[theta] Cos[phi], Sin[theta]Cos[phi], Sin[phi]}];
 NumRootsGlynn[{x_, y_, z_}] := If[z^2 > x^2 + y^2, 3, If[x < 0, 2, 1]];
 TriplexRootGlynn[{x_, y_, z_}, k_] := Module[{r = Sqrt[x^2 + y^2 + z^2], ktheta, kphi, theta, phi}, {ktheta, kphi} = If[k == 0, {0, 0}, If[z^2 > x^2 + y^2, If[k == 1, If[x <0 && y < 0, {2, 0}, {0.5, If[z > 0, 0.5, 2.5]}], If[x < 0 && y > 0, {1, 0}, {2.5, If[z > 0, 0.5, 2.5]}]], {If[y < 0,
 2, 1], 0}]]; theta = (ArcTan[x, y] + ktheta Pi)/power; phi = (ArcSin[z/r] +kphi Pi)/power; r^(1/power){Cos[theta] Cos[phi], Sin[theta]Cos[phi], Sin[phi]}];
 c = {-0.2, 0, 0}; imax = 100; dpmax = 50.0; p = Table[{-0.61, 0.0, -0.42}, {imax}]; dp = Table[1, {imax}]; roots = Table[0, {imax}]; plist = {}; i = 2;
 While[i > 1 && Length[plist] < 10000, p[[i]] = TriplexRootGlynn[p[[i - 1]] - c, roots[[i]]]; dp[[i]] = Abs[power]Sqrt[norm[TriplexPow[p[[i]], power - 1]]]dp[[i - 1]];plist = Append[plist, p[[i]]]; prune = (i == imax || dp[[i]] > dpmax); If[prune, While[i > 1 && roots[[i]] == NumRootsGlynn[p[[i - 1]] - c] - 1, roots[[i]] = 0; i--]; roots[[i]]++, i++; roots[[i]] = 0]];
 Show[Graphics3D[{PointSize[0.005], Point /@ plist}]]
 
 Links
 Glynn Juliabulb - interesting variation based on the cosine power formula, looks like a "weaved" clam shell, by Jos Leys
 Tripcove Fractal - another interesting variation probably based on the cosine power formula, by Dizingof
 Lichenfan Lamp Shade - another interesting variation probably based on the cosine power formula, by unellenu
 
 | 
 4D Bicomplex Mandelbrot Set ("Tetrabrot") - C++, 7/30/09
 
 |    {x,y,z,w}2 = {x2-y2-z2+w2, 2(xy-zw), 2(xz-yw), 2(xw+yz)} This is the 4D bicomplex Mandelbrot set, also called the Tetrabrot set. The shape reminds me of a Bismuth crystal. Here is some Mathematica code for the first miniature Tetrabrot set:
 (* runtime: 2.6 minutes, increase n for higher resolution *)
 n = 100; norm[x_] := x.x; BicomplexSquare[{x_, y_, z_, w_}] := {x^2 - y^2 - z^2 + w^2, 2 (x y - z w), 2 (x z - y w), 2 (x w + y z)};
 Mandelbrot4D[c_] := Module[{q = {0, 0, 0, 0}, i = 0}, While[i < 20 && norm[q] < 4, q = BicomplexSquare[q] + c; i++]; i];
 image = Table[z = 0.06; While[z >= 0 && Mandelbrot4D[{x, y, z, 0}] < 20, z -=0.12/n]; If[z < 0, -0.06, z], {y, -0.06, 0.06, 0.12/n}, {x, -1.82, -1.7, 0.12/n}];
 ListDensityPlot[image, Mesh -> False, Frame -> False, PlotRange -> {-0.02, 0.06}]
 
 Links
 TetraBrot Explorer - deep zooming program by Etienne Martineau and Dominic Rochon, see the animations
 similar renderings - by Jason McGuiness
 
 | 
 |    The image on the left is the third order variation of the Tetrabrot set. The photo on the right is a Bismuth crystal. 
 | 
 4D Bicomplex Julia Set - C++, 7/29/09
 
 |  {x,y,z,w}2 = {x2-y2-z2+w2, 2(xy-zw), 2(xz-yw), 2(xw+yz)} 
 | 
 Inverse 4D Bicomplex Julia Set - C++, 9/21/09
 
 |        These images were calculated using the inverse Julia set method. Dominic Rochon helped me with the formula for finding the square root of a bicomplex number. This formula has 4 roots, so the total number of points increases by a factor of 4 after each iteration. Here is some Mathematica code: (* runtime: 0.7 second *)
 SqrtW[w1_, w2_] := 0.5{Re[w1] + Re[w2], Im[w1] + Im[w2], Im[w2] - Im[w1], Re[w1] - Re[w2]};
 BicomplexSqrts[{a_, b_, c_, d_}] := Module[{z1 = a + I b, z2 = c + I d, w1, w2}, w1 = Sqrt[z1 - I z2]; w2 = Sqrt[z1 + I z2]; {SqrtW[w1, w2], SqrtW[w1, -w2],SqrtW[-w1, w2], SqrtW[-w1, -w2]}];
 wc = {-0.2, 0.8, 0, 0}; wlist = {{1.0, 1.0, 1.0, 1.0}};
 Do[wlist = Flatten[Map[BicomplexSqrts[# - wc] &, wlist], 1], {6}];
 Show[Graphics3D[{PointSize[0.005], {Hue[#[[4]]], Point[Delete[#, 4]]} & /@ wlist}, Boxed -> False, Background -> RGBColor[0, 0, 0], ViewPoint -> {0, 10, 0}]]
 
 | 
 4D Julibrot Set - C++, 7/6/09
 
 |      The complex Julia set is typically calculated in 2D with a constant complex seed value. However, if you consider the set of all possible Julia sets (all possible seed values), you get the complete 4D Julibrot set. The Mandelbrot set is contained as a 2D slice of the 4D Julibrot set. You can also take 3D slices of the 4D Julibrot set as shown here. The left picture shows what the Mandelbrot set looks like if you vary the real part of the starting value along the third dimension and the right picture shows what it looks like when you vary the imaginary part. Technically, this is not a hypercomplex fractal, but it's still worth mentioning here. Here is some Mathematica code: (* runtime: 18 seconds, increase n for higher resolution *)
 n = 100; Julibrot[z0_, zc_] := Module[{z = z0, i = 0}, While[i <12 && Abs[z] < 2, z = z^2 + zc; i++]; i];
 image = Table[z = 1.5; While[z >= 0 && Julibrot[z, x + I y] < 12, z -= 3/n]; z, {y, -1.5, 1.5, 3/n}, {x, -2, 1, 3/n}];
 ListDensityPlot[image, Mesh -> False, Frame -> False,PlotRange -> {0, 1.5}]
 
 | 
 4D Hopfbrot Set - C++, 12/2/09
 
 3D Exponential Mandelbrot Set - C++, 7/6/09
 
 |    This is what the Mandelbrot set looks like if you vary the exponent along the third dimension. It has a twisted shape because the features of the Mandelbrot set tend to rotate around the origin as the exponent is increased. Technically, this is not a hypercomplex fractal, but it's still worth mentioning here. Here is some Mathematica code: (* runtime: 49 seconds, increase n for higher resolution *)
 n = 100; Mandelbrot[n_, zc_] := Module[{z = 0, i = 0}, While[i <12 && Abs[z] < 2, z = z^n + zc; i++]; i];
 image = Table[y = 1.5; While[y >= 0 && Mandelbrot[z, x + I y] < 12, y -= 3/n]; y, {z, 1, 4, 3/n}, {x, -2, 1, 3/n}];
 ListDensityPlot[image, Mesh -> False, Frame -> False, PlotRange -> {0, 1.5}]
 
 | 
 4D Mandelbulb Set - C++, 7/7/09
 
 |    Here is a 4D version of Daniel White's squarring formula. In 4D, rotation is about a plane and there are 6 possible rotational matrices to choose from: Rxy, Ryz, Rxz, Rxw, Ryw, Rzw. Following this train of thought, I came up with the following 4D analog to Daniel's formula by applying three consecutive 4D rotations: {x,y,z,w}2 = Rxy(2θ)Rxz(2φ)Rxw(2ψ){r2,0,0,0}, r2=sqrt(x2+y2), r3=sqrt(x2+y2+z2), r=sqrt(x2+y2+z2+w2), θ=atan2(y,x), φ=atan2(z,r2), ψ=atan2(w,r3)
 This can be expanded to give:
 {x, y, z, w}2 = r2{cos(2ψ)cos(2φ)cos(2θ), cos(2ψ)cos(2φ)sin(2 θ), -cos(2ψ)sin(2φ), sin(2ψ)}
 This formula reduces to Daniel's squaring formula when w = 0 (which, in turn, reduces to the regular complex squaring formula when z = 0). For faster calculations, this formula can be simplified to:
 {x, y, z, w}2 = {(x2-y2)b, 2xyb, -2r2za, 2r3w}, a=1-w2/r32, b=a(1-z2/r22)
 Using this formula, I rendered 3D slices of this 4D Mandelbrot fractal at w = 0, z = 0, y = 0, x = 0. The animation shows the 3D slice starting from w = 0, and then rotating to z = 0, then to y = 0, then to x = 0, and finally back to w = 0 again.
 
 | 
 3D Nebulabrot - C++, 9/9/09
 
 |    Here are some volumetric renderings. This still has some bugs in it, but I think this will look quite beautiful if I ever get the bugs ironed out. 
 | 
 3D Mandelbrot Set Orbit Traps - C++, 8/18/09
 
 |      Orbit trapping is a popular technique for mapping arbitrary shapes to 2D fractals. These images were created by extending this concept into 3 dimensions. These 3D orbit traps are shaped like spheres. See also my 3D Pickover stalks and 2D orbit trap fractals. Here is some Mathematica code: (* runtime: 2 minutes, increase n for higher resolution *)
 n = 100; norm[x_] := x.x; TricornSquare[{x_, y_, z_}] := {x^2 - y^2 - z^2, 2x y, -2x z};
 Clear[OrbitTrap]; OrbitTrap[p_] := 1.0 - norm[p]/0.3^2;
 Mandelbrot[c_] := Module[{p = {0, 0, 0}, i = 0}, While[i < 10 && norm[p] < 4 && (i < 2 || OrbitTrap[p] < 0.0), p =TricornSquare[p] + c; i++]; OrbitTrap[p]];
 image2 = Table[z = 1.5; While[z > -1.5 && Mandelbrot[{x, y, z}] < 0.0, z -= 3/n]; z, {y, -1.5, 1.5, 3/n}, {x, -2.0, 1.0, 3/n}];
 ListDensityPlot[image2, Mesh -> False, Frame -> False, PlotRange -> All]
 
 Links
 square orbit trap - uses formula for 8th order Mandelbrot set, by Enforcer
 3D Mandelbrot set orbit trap - animation by David Makin
 Quaternion Eggs - I think this is a 3D orbit trap, by Ramiro Perez
 
 | 
 3D Mandelbrot Set, based on Rudy Rucker's formula for squaring a 3D hypercomplex number - C++, 9/6/09
 
 3D Christmas Tree Mandelbrot Set - C++, 12/2/09
 
 |    {x,y,z}n = rn{cos(θ),sin(θ)cos(φ),sin(θ)sin(φ)} r=sqrt(x2+y2+z2), θ=n cos-1(x/r), φ=n atan2(z,y)
 This 3D Mandelbrot set is based on a power formula that travels the same distance around the sphere no matter what the angle. The Cosine method also satisfies this requirement but one advantage of this method is that it reduces to complex arithmetic in the limit when z = 0 and therefore it contains the regular 2D complex Mandelbrot set as a cross-section in the x-y plane. This formula also happens to be identical to the 3D slice at y = 0 of my 4D Hopfbrot formula.
 
 | 
 3D Mandelbrot Set based on my attempt to combine the forward method with a reversed method - C++, 11/20/09
 
 |    {x,y,z}n = rn{cos(θ)cos(φ), sin(θ)cos(φ), cos(θ)sin(φ)} r=sqrt(x2+y2+z2), θ=n atan2(y,x), φ=n sin-1(z/r)
 This was my attempt to combine the forward method with a reversed order method. The 8th order version looks like an explosion.
 
 
 | 
 4D "Roundy" Mandelbrot Set - C++, 7/24/09
 
 | | Mandelbrot Set, imax=24 
  | 
 {x,y,z,w}2 = {x2-y2-z2-w2, 2(xy+zw), 2(xz+yw), 2(xw+yz)}| Minibrot, imax=21 
  | 
 For a brief time, this formula was a popular method for creating 3D Mandelbrot sets. Here is some Mathematica code:
 (* runtime: 1.2 minutes, increase n for higher resolution *)
 n = 100; norm[x_] := x.x; RoundySquare[{x_, y_, z_, w_}] := {x^2 - y^2 - z^2 - w^2, 2 (x y + z w), 2 (x z + y w), 2 (x w + y z)};
 Mandelbrot4D[c_] := Module[{q = {0, 0, 0, 0}, i = 0}, While[i < 24 && norm[q] < 4, q = RoundySquare[q] + c; i++]; i];
 image = Table[z = 1.5; While[z >= 0 && Mandelbrot4D[{x, y, z, 0}] < 24, z -= 3/n]; If[z < 0, -1.5, z], {y, -1.5, 1.5, 3/n}, {x, -2.0, 1.0, 3/n}];
 ListDensityPlot[image, Mesh -> False, Frame -> False, PlotRange -> {-0.1,1.5}]
 
 Links
 minibrot, Mandelbrot, Genuine Mandelbrot 3D, Minibrot 3D, New Mandelbrot 3D, True 3D Julia, Patch with Julia Fractals - volumetric renderings, by Krzysztof Marczak, be sure to see his animation
 Space Station - volumetric mini-minibrot, by Krzysztof Marczak
 3D Mandlebrot - by Thomas Ludwig
 4D Mandy, Minibrot - by David Makin
 
 | 
 4D "Roundy" Julia Set - C++, 7/28/09
 
 Bristorbrot Set - C++, 12/4/09
 
 3D Mandelbrot Set, based on David Makin's formula for squaring a 3D hypercomplex number - C++, 10/14/09
 
 3D Mandelbulb / Juliabulb Sets, based on Roger Bagula's parameters - C++, 12/15/09
 
 3D Mandelbrot Set - C++, 12/9/09
 
 |      The original Mandelbulb triplex power formula applies the azimuthal angle (φ) rotation away from the azimuthal axis (z-axis). However, one day I got confused and I thought that I had to sandwich the y-axis rotation between an un-rotate / re-rotate pair of transformations. Hence was born this formula. 
 Links
 Theme Park Ride - by David Makin
 
 | 
 3D Mandelbulb Orbits - POV-Ray 3.1, 1/5/10
 
 4D "Squarry" Mandelbrot Set - C++, 7/30/09
 
 4D "Squarry" Julia Set - C++, 7/29/09
 
 |  {x,y,z,w}2 = {x2-y2-z2-w2, 2(xy+zw), 2(xz+yw), 2(xw-yz)} c={-0.2,0.8,0,0}, imax=12
 
 | 
 4D "Mandy Cousin" Mandelbrot Set - C++, 7/30/09
 
 4D "Mandy Cousin" Julia Set - C++, 7/29/09
 
 |  {x,y,z,w}2 = {x2-y2-z2+w2, 2(xy+zw), 2(xz+yw), 2(xw+yz)} c={-0.2,0.8,0,0}, imax=12
 
 | 
 4D Mandelbrot Set Variation - C++, 7/30/09
 
 | | Mandelbrot Set, imax=12 
  | 
 {x,y,z,w}2 = {x2-y2-z2-w2, 2(xy+zw), 2(xz-yw), 2(xw+yz)}| Minibrot, imax=20 
  | 
 Here is another formula for squaring a 4D hypercomplex number that David Makin suggested to me.
 
 | 
 4D Julia Set Variation - C++, 7/29/09
 
 |  {x,y,z,w}2 = {x2-y2-z2-w2, 2(xy+zw), 2(xz-yw), 2(xw+yz)} c={-0.2,0.8,0,0}, imax=12
 
 | 
 Quadray Quaternion Julia Set - C++, 10/15/09
 
 BoneYard - C++
 
 |  Here are some other weird looking hypercomplex fractals that I made while experimenting with different formulas. The possibilities are endless, but the challenge is to find formulas that produce intricate 3D details without looking like stretched out taffy. | 
 Other Hypercomplex Fractal Links
 3D Fractals - created by Richard Rosenman using Terry Gintz' QuaSZ software, includes "quaternion, hypercomplex, cubic Mandelbrot, complexified quaternion and octonion"
 
 
 |