A better 3D generalization of the Mandelbrot set.

February 9, 2021 math and physics play , , , , , , ,

I’ve been exploring 3D generalizations of the Mandelbrot set:

The iterative equation for the Mandelbrot set can be written in vector form ([1]) as:
\begin{equation}
\begin{aligned}
\Bz
&\rightarrow
\Bz \Be_1 \Bz + \Bc \\
&=
\Bz \lr{ \Be_1 \cdot \Bz }
+
\Bz \cdot \lr{ \Be_1 \wedge \Bz }
+ \Bc \\
&=
2 \Bz \lr{ \Be_1 \cdot \Bz }
–
\Bz^2\, \Be_1
+ \Bc
\end{aligned}
\end{equation}
Plotting this in 3D was an interesting challenge, but showed that the Mandelbrot set expressed above has rotational symmetry about the x-axis, which is kind of boring.

If all we require for a 3D fractal is to iterate a vector equation that is (presumably) at least quadratic, then we have lots of options. Here’s the first one that comes to mind:
\begin{equation}
\begin{aligned}
\Bz
&\rightarrow
\gpgradeone{ \Ba \Bz \Bb \Bz \Bc } + \Bd \\
&=
\lr{ \Ba \cdot \Bz } \lr{ \Bz \cross \lr{ \Bc \cross \Bz } }
+
\lr{ \Ba \cross \Bz } \lr{ \Bz \cdot \lr{ \Bc \cross \Bz } }
+ \Bd
.
\end{aligned}
\end{equation}
where we iterate starting, as usual with \( \Bz = 0 \) where \( \Bd \) is the point of interest to test for inclusion in the set. I tried this with
\begin{equation}\label{eqn:mandel3:n}
\begin{aligned}
\Ba &= (1,1,1) \\
\Bb &= (1,0,0) \\
\Bc &= (1,-1,0).
\end{aligned}
\end{equation}
Here are some slice plots at various values of z

and an animation of the slices with respect to the z-axis:

Here are a couple snapshots from a 3D Paraview rendering of a netCDF dataset of all the escape time values

Data collection and image creation used commit b042acf6ab7a5ba09865490b3f1fedaf0bd6e773 from my Mandelbrot generalization experimentation repository.

References

[1] L. Dorst, D. Fontijne, and S. Mann. Geometric Algebra for Computer Science. Morgan Kaufmann, San Francisco, 2007.

Slicing of the 3D Mandelbrot set, and analysis.

February 8, 2021 math and physics play , , , , , , ,

Some slices.

As followup to:

here is a bit more experimentation with Paraview slice filtering. This time I saved all the point data with the escape time counts, and rendered it with a few different contours

The default slice filter places the plane in the x-y orientation:

but we can also tilt it in the Paraview render UI

and suppress the contour view to see just the slice

As a very GUI challenged user, I don’t find the interface particularly intuitive, but have at least figured out this one particular slicing task, which is kind of cool.  It’s impressive that the UI can drive interesting computational tasks without having to regenerate or reload any of the raw data itself.  This time I was using the MacOSX Paraview client, which is nicer looking than the Windows version, but has some weird glitches in the file dialogues.

Analysis.

The graphing play above shows some apparent rotational symmetry our vector equivalent to the Mandelbrot equation
\begin{equation}
\Bx \rightarrow \Bx \Be_1 \Bx + \Bc.
\end{equation}
It was not clear to me if this symmetry existed, as there were artifacts in the plots that made it appear that there was irregularity. However, some thought shows that this irregularity is strictly due to sampling error, and perhaps also due to limitations in the plotting software, as such an uneven surface is probably tricky to deal with.

To see this, here are the first few iterations of the Mandlebrot sequence for an arbitary starting vector \( \Bc \).
\begin{equation}
\begin{aligned}
\Bx_0 &= \Bc \\
\Bx_1 &= \Bc \Be_1 \Bc + \Bc \\
\Bx_2 &= \lr{ \Bc \Be_1 \Bc + \Bc } \Be_1 \lr{ \Bc \Be_1 \Bc + \Bc } + \Bc \Be_1 \Bc + \Bc.
\end{aligned}
\end{equation}

Now, what happens when we rotate the starting vector \( \Bc \) in the \( y-z \) plane. The rotor for such a rotation is
\begin{equation}
R = \exp\lr{ e_{23} \theta/2 },
\end{equation}
where
\begin{equation}
\Bc \rightarrow R \Bc \tilde{R}.
\end{equation}
Observe that if \( \Bc \) is parallel to the x-axis, then this rotation leaves the starting point invariant, as \( \Be_1 \) commutes with \( R \). That is
\begin{equation}
R \Be_1 \tilde{R} =
\Be_1 R \tilde{R} = \Be_1.
\end{equation}
Let \( \Bc’ = R \Bc \tilde{R} \), so that
\begin{equation}
\Bx_0′ = R \Bc \tilde{R} = R \Bx_0 \tilde{R} .
\end{equation}
\begin{equation}
\begin{aligned}
\Bx_1′
&= R \Bc \tilde{R} \Be_1 R \Bc \tilde{R} + R \Bc \tilde{R} \\
&= R \Bc \Be_1 \Bc \tilde{R} + R \Bc \tilde{R} \\
&= R \lr{ \Bc \Be_1 \Bc R + \Bc } \tilde{R} \\
&= R \Bx_1 \tilde{R}.
\end{aligned}
\end{equation}
\begin{equation}\label{eqn:m2:n}
\begin{aligned}
\Bx_2′
&= \Bx_1′ \Be_1 \Bx_1′ + \Bc’ \\
&= R \Bx_1 \tilde{R} \Be_1 R \Bx_1 \tilde{R} + R \Bc \tilde{R} \\
&= R \Bx_1 \Be_1 \Bx_1 \tilde{R} + R \Bc \tilde{R} \\
&= R \lr{ \Bx_1 \Be_1 \Bx_1 + \Bc } \tilde{R} \\
&= R \Bx_2 \tilde{R}.
\end{aligned}
\end{equation}

The pattern is clear. If we rotate the starting point in the y-z plane, iterating the Mandelbrot sequence results in precisely the same rotation of the x-y plane Mandelbrot sequence. So the apparent rotational symmetry in the 3D iteration of the Mandelbrot vector equation is exactly that. This is an unfortunately boring 3D fractal. All of the interesting fractal nature occurs in the 2D plane, and the rest is just a consequence of rotating that image around the x-axis. We get some interesting fractal artifacts if we slice the rotated Mandelbrot image.

Some 3D renderings of the Mandelbrot set.

February 7, 2021 math and physics play , , , , ,

As mentioned previously, using geometric algebra we can convert the iterative equation for the Mandelbrot set from complex number form
\begin{equation}
z \rightarrow z^2 + c,
\end{equation}
to an equivalent vector form
\begin{equation}
\mathbf{x} \rightarrow \mathbf{x} \mathbf{e} \mathbf{x} + \mathbf{c},
\end{equation}
where \( \mathbf{e} \) represents the x-axis (say). Geometrically, each iteration takes \( \mathbf{e} \) and reflects it about the direction of \( \mathbf{x} \), then scales that by \( \mathbf{x}^2 \) and adds \( \mathbf{c} \).

To get the usual 2D Mandelbrot set, one iterates with vectors that lie only in the x-y plane, but we can treat the Mandelbrot set as a 3D solid if we remove the x-y plane restriction.

Last time I animated slices of the 3D set, but after a whole lot of messing around I managed to save the data for all the interior points of the 3D set in netCDF format, and render the solid using Paraview. Paraview has tons of filters available, and experimenting with them is pretty time consuming, but here are some initial screenshots of the 3D Mandelbrot set:

It’s interesting that much of the characteristic detail of the Mandelbrot set is not visible in the 3D volume, but if we slice that volume, you can then you can see it again.  Here’s a slice taken close to the z=0 plane (but far enough that the “CN tower” portion of the set is not visible)

You can also see some of that detail if the opacity of the rendering is turned way down:

If you look carefully at the images above, you’ll see that the axis labels are wrong.  I think that I’ve screwed up one of the stride related parameters to my putVar call, and I end up with x+z transposed in the axes labels when visualized.

Visualizing the 3D Mandelbrot set.

February 1, 2021 C/C++ development and debugging.

In “Geometric Algebra for Computer Science” is a fractal problem based on a vectorization of the Mandelbrot equation, which allows for generalization to \( N \) dimensions.

I finally got around to trying the 3D variation of this problem.  Recall that the Mandlebrot set is a visualization of iteration of the following complex number equation:
\begin{equation}
z \rightarrow z^2 + c,
\end{equation}
where the idea is that \( z \) starts as the constant \( c \), and if this sequence converges to zero, then the point \( c \) is in the set.

The idea in the problem is that this equation can be cast as a vector equation, instead of a complex number equation. All we have to do is set \( z = \Be_1 \Bx \), where \( \Be_1 \) is the x-axis unit vector, and \( \Bx \) is an \(\mathbb{R}^2\) vector. Expanding in coordinates, with \( \Bx = \Be_1 x + \Be_2 y \), we have
\begin{equation}
z
= \Be_1 \lr{ \Be_1 x + \Be_2 y }
= x + \Be_1 \Be_2 y,
\end{equation}
but since the bivector \( \Be_1 \Be_2 \) squares to \( -1 \), we can represent complex numbers as even grade multivectors. Making the same substitution in the Mandlebrot equation, we have
\begin{equation}
\Be_1 \Bx \rightarrow \Be_1 \Bx \Be_1 \Bx + \Be_1 \Bc,
\end{equation}
or
\begin{equation}
\Bx \rightarrow \Bx \Be_1 \Bx + \Bc.
\end{equation}
Viola! This is a vector version of the Mandlebrot equation, and we can use it in 2 or 3 or N dimensions, as desired.  Observe that despite all the vector products above, the result is still a vector since \( \Bx \Be_1 \Bx \) is the geometric algebra form of a reflection of \( \Bx \) about the x-axis.

The problem with generalizing this from 2D is really one of visualization. How can we visualize a 3D Mandelbrot set? One idea I had was to use ray tracing, so that only the points on the surface from the desired viewpoint need be evaluated. I don’t think I’ve ever written a ray tracer, but I thought that there has got to be a quick and dirty way to do this.  Also, figuring out how to make a ray tracer interact with an irregular surface like this is probably non trivial!

What I did instead, was a brute force evaluation of all the points in the upper half plane in around the origin, one slice of the set at a time. Here’s the result

Code for the visualization can be found in github. I’ve used Pauli matrices to do the evaluation, which is actually pretty quick (but slower than plain std::complex< double> evaluation), and the C++ ImageMagick API to save individual png files for the slices. There are better coloring schemes for the Mandelbrot set, and if I scrounge up some more time, I may try one of those instead.

Example of PL/I macro

January 27, 2021 Mainframe , , ,

Up until last week PL/I macros were a bit of a mystery.  Most of the ones that I’d seen in customer code were impressively inscrutable, and if I had to look at any of them, my reaction was to throw my hands in the air and plead with the compiler backend guys for help.  Implementing one such macro has been very helpful to understanding how these work.

Here is a C program that roughly models some PL/I code of interest

The documentation for the ‘foo’ function says of the final return code parameter that it is 12 bytes long, and that the ‘rcvalues.h’ header file has a set of RCNNN constants and a RCCHECK macro that can be used to test for any one of those constants.  A possible C implementation of that header might look something like:

/* rcvalues.h */
#define RC000 0x0000000000000000LL
#define RC001 0x0000000123456789LL
/* ... */

#define RCCHECK( urc, crc ) ( memcmp( &(urc), &(crc), 8 ) == 0 )

PL/I APIs do not typically use modern constructs like typedefs.  The closest that I have seen is for an API header file (copybook in the mainframe lingo) is to declare a variable (which becomes a local variable with a specific name in the including module), which the programmer can refer to using the LIKE keyword, as in the following example:

I believe there is also a DEFINE keyword available in newer PL/I compilers, which provides a typedef like mechanism, but most existing code probably doesn’t use such new-fangled nonsense, when cut and paste has far superior maintenance characteristics.  For that reason, the API would be unlikely to have a typedef equivalent for the return code structure.  Instead, the PL/I equivalent of the C code above, would probably look like:

(i.e. the C code is really modeled on the PL/I code of this form, and if this was a C API, the API would have a struct declaration or a typedef for the return code structure)

The RCNNN constants would actually be found as named variables (not immutable constants) in the copy book, perhaps declared something like:

I struggled a bit to figure out what the PL/I equivalent of my C RCCHECK macro would be.  The following inner function correctly did the required type casting and comparisons:

The implementation is very long, since the entire declaration of the input parameter type has to be duplicated.

If I was to put this RCCHECK implementation above into my RCVALUES.inc header file, it would only work if all the customer declaration of their return code structure objects were field by field compatible.  What I really want is for my RCCHECK function to take the address of the parameter, and pass that instead of the underlying type.  That was not at all obvious to figure out how to do, but with some help, I was eventually able to construct a PL/I macro (with helper inner-function) of the following form:

It’s clearly no longer a one liner.  Some notes on this PL/I macro:

  • The PL/I macro body looks like a regular PL/I function, but the begin-PROCEDURE and END statements start with % (% is not part of the PROC name.)
  • Macro parameters and return values are explicit strings, regardless of the types of the parameters that were actually passed.
  • In PL/I the || symbol is used for string concatenation, so this constructs output that inserts an ADDR() call around ARG1 token and then passes the ARG2 token as is.
  • I don’t know if there’s a way to implement this macro in a way that doesn’t require a helper function, and still have the output work in the context of an IF statement.
  • You have to explicitly enable the macro, using %ACTIVATE.  In my case, without %ACTIVATE, the RCCHECK symbol ends up as an undeclared external entry, and was no call to the @RCCHECK_HELPER function \({}^{[1]}\).
  • Observe that the PL/I macro provides a mechanism to jam whatever you want into the code, as the compiler’s macro preprocessor replaces the macro call tokens with the string that you have provided, leaving that string for the final compiler pass to interpret instead.

If I compile the code using this macro version of RCCHECK, the preprocessor output looks like:

I’m still pretty horrified at some of the macros that I’ve seen in customer code — they almost seem like the source equivalent of self modifying code.  You can’t figure out what is going on without also looking at all the output of the precompiler passes.  This is especially evil, since you can write PL/I preprocessor macros that generate preprocessor macros and require multiple preprocessor passes to produce the final desired output!

Footnotes

[1] note that @ is a valid PL/I character to use in a symbol name, as is # and $ — so if you want your functions to look like swear words, this is a language where that is possible.  Something like the following is probably valid PL/I :

V = #@$1A#@(1);

For added entertainment, your file names (i.e. PDS member names) can also be like ‘#@$1A#@’. Storing files with names like that on a Unix filesystem results in hours of fun, as you are then left with the task of figuring out how to properly quote file names with embedded $’s and #’s in scripts and makefiles.