Reenigne blog

I have been poking again at my fractal plotting program recently. The latest version is available here (source). I've added multi-threading (so all available cores will be used for calculation) and Direct3D rendering (which makes zooming much smoother and allows us to rotate in real time as well with the Z and X keys - something I've not seen any other fractal program do). There is still some jerkiness when zooming in or out (especially past a power-of-two boundary) but it's greatly reduced (especially when there's lots of detail) and with a little more work (breaking up the global matrix operations) it should be possible to eliminate it altogether.

The linked grid data structure in the last version I posted worked reasonably well, but it used too much memory. At double-precision, every block needed 72 bytes even after it had been completed. This means that an image that covered a 1600x1200 screen with detail at 16 texels per pixel for anti-aliasing would use more than 2Gb. I could use a different block type for completed points but there would still need to be at least 4 pointers (one for each direction) at each point, which could take 500Mb by themselves.

So I decided to switch to a data structure based on a quadtree, but with some extra optimizations - at each node, one can have not just 2x2 children but 2ⁿ by 2ⁿ for any positive integer n. When a node contains 4 similar subtrees (each having the same size and subnode type) it is consolidated into a 2ⁿ⁺¹ by 2ⁿ⁺¹ grid. If we need to change the type of one of the subnodes of this grid, we do the reverse operation (a split). In this way, if we have a large area of completed detail it will use only 4 bytes per point (storing the final number of iterations). This works much better - memory usage can still be quite high during processing (but seems to top out at ~700Mb rather than ~2Gb for the linked grid data structure) and goes right back down to ~200Mb when the image is completed (only about 35Mb of which is the data matrix, the rest is Direct3D textures).

One disadvantage of the grid tree method is that it doesn't naturally enforce the rule preventing 4 or more blocks being along one side of one other block. That rule is important for ensuring that we don't lose details that protrude through a block edge. However, this turns out to be pretty easy to explicitly add back in. This gives a more efficient solid-guessing method than I've seen in any other fractal program (by which I mean it calculates the fewest possible points in areas of solid colour).

This data structure was surprisingly difficult to get right. The reason for this is that with all the splitting and consolidating of grids, one must avoid keeping grid pointers around (since they can get invalidated when a split or consolidation happens in a nearby part of the tree). That means that many algorithms that would be natural to write recursively need to be written iteratively instead (since the recursive algorithm would keep grid pointers on the call stack).

Another tricky thing to get right was the memory allocation. My first attempt just used malloc() and free() to allocate and deallocate memory for the grids. However, the way grids are used is that we allocate a whole bunch of them and then deallocate most of them. However, they are quite small and the ones that don't get deallocated are spread throughout memory. This means that very little memory can be reclaimed by the operating system, since nearly every memory page has at least one remaining block (or part of one) on it.

The problem here is that the malloc()/free() interface assumes that a block, once allocated, never moves. However, grids turn out to be quite easy to move compared to general data structures. This is because the grid tree is only modified by one thread at once (it must be, since modifying it can invalidate any part of it) and also because we know where all the grid pointers are (we have to for split and consolidate operations). Only a dozen or so different grid sizes are currently used in practice, so it's quite practical to have a separate heap for each possible grid size. If all the objects in a heap are the same size, there's an extremely simple compaction algorithm - whenever we delete an item, just move the highest item in the heap down into the newly free space. So a heap consists of a singly-linked list of items. These are allocated in chunks of about 1Mb. Each chunk is contiguous region of memory, but chunks within a heap don't have to be contiguous. This allows the heaps to interleave and make the best use of scarce (on 32-bit machines) address space.

Sounds simple enough but again there turned out to be complications. Now it's not just split and consolidate that can invalidate grid pointers - now deleting a grid can also invalidate grid pointers. This means more recursive algorithms that need to be made iterative (such as deleting an entire subtree) and certain operations need to be done in a very specific order to avoid keeping more grid pointers around than are necessary. I now understand why such techniques are so rarely used - they're really difficult to get right. I do think that a general garbage-collected heap would be slower and use more memory than this special purpose allocator, but I'm not sure if difference would be large enough to make all this effort worth it. Perhaps someone will port my code to a garbage-collected language so the performance can be compared.

I'd like to be able to add some kind of sound output to my fractal program. What I'd like to do is generate a click every time a point is completed. So, when it's working on a difficult area it would emit a clicking sound but when it's finishing a lot of points it would generate white noise. The principle is very similar to that of a geiger counter's audio output.

This is trivial to do under DOS, but in Windows it is surprisingly difficult. Ideally we want to generate a waveform at 44.1KHz. When we finish a point, we flip the current sample in the audio buffer between its maximum and minimum positions. The problem is determining what the current sample should be (or, alternatively, moving a pointer to the next sample 44,100 times per second). What we need is a high-resolution timer, but none of the high-resolution timers seem to be particularly reliable. timeGetTime() is too low-resolution. QueryPerformanceCounter() has bugs on some systems which cause it to skip backwards and forwards a few seconds on occasion. The RDTSC CPU instruction isn't guaranteed to be synchronized between CPUs, and its count rate may vary due to power saving features on some CPUs. Another possibility is reading the playback cursor position from DirectSound, but that is unreliable on OSes before Vista.

What I may end up doing is faking it - measure the average number of blocks completed over some small period of time and use that to generate the clicks via a Poisson distribution.

Lots of videos have been made of zooms into the Mandelbrot set. These videos almost all follow the same pattern. The initial image is the familiar one of the entire set. One zooms in on some tendril or other, gathering details, passing various minibrots and possibly some embedded Julia sets on the way until one reaches a final minibrot and the zoom stops.

The reason for this is that these sequences are relatively quick to calculate, they show a variety of the different types of behaviors that the Mandelbrot set has to offer, and they have a natural endpoint which causes the video not to seem like it has been cut off part-way through.

If we relax these conditions, what interesting zoom sequences open up to us?

To answer this question requires classification of the complex numbers according to the roles they play in the Mandelbrot set.

The first taxon separates exterior points (whose z sequences are unbounded, example: +1), interior points (limit of z sequence is a finite attracting cycle, example: 0) and boundary points (all other z sequence behaviours). Zooming into either of the first two will leave us with an image of uniform colour after a finite amount of time. Most zoom movies actually target the second of these.

Zooming in on the boundary will leave us with an infinite sequence of "interesting" images, so let's explore this group further. The second taxon separates shore points (those which are on the edge of an atom), Feignbaum points (which are at the limit of a sequence of atoms of increasing period) and filament points (all other boundary points).

Feignbaum points make quite interesting infinite zoom sequences. The canonical example is ~-1.401155, the westernmost point of the main molecule. Each time you zoom in by a factor of ~4.669201, the lakes look the same but the the filaments get denser. Filaments connect to molecules at these points. These are the only places filaments attach to the continent but each mu-molecule has one additional filament emerging from its cardioid cusp.

Filament points can be divided into two types based on orbit dynamics - those whose orbit is chaotic (irrational external angle) and those whose orbit converges to a finite repelling cycle. The latter includes two subtypes - branch points (where there are more than one but fewer than infinity paths to other points in the set) and terminal points (where there is only one). Well-known terminal points include -2 and +i. Zooming into terminal points and branch points gives a video which is infinite but which devolves into very same-y spirals after a short time as the mu-molecules shrink to less than a pixel.

I don't know an exact formula for any chaotic filament points, but they're easy to find using a fractal plotter - just zoom in, keeping away from any mu-molecules, branch points and terminal points. As with terminal and branch points, the minibrots tend to shrink to less than a pixel quite quickly, so the images obtained will become same-y after a while - the structures formed by the filaments tend to look like L-system fractals.

Shore points can be divided into two categories, parabolic points and Siegel points. Parabolic points come in two flavors, cardioid cusps (the point of maximum curvature of a cardioid, example: +0.25) and bond points (where two atoms meet, example: -0.75). If you zoom in a cusp you'll technically always have detail on the screen but the picture gets very boring very quickly as the cusps shrink to sub-pixel widths.

This leaves the Siegel points. You can find these by zooming into the edge of an atom but avoiding the cusps. For example:

1. Start off with the large circle to the west (call it A) and the slightly smaller circle to the north (call it B).
2. Zoom in until you can pick out the largest circle C between A and B on the boundary of the set.
3. Rename B as A and C as B.
4. Go to step 2.

This particular algorithm takes you to the point whose Julia set is known as the Siegel disk fractal, which is about -.390541+.586788i.

All points on the boundary of the continent cardioid can be found using the formula . For rational , this gives a parabolic point and for irrational this gives a Siegel point (a point whose Julia set has Siegel discs). The "most irrational" number (in the sense that it is least well approximated by rational numbers) is the Golden ratio . The Siegel point corresponding to this number is the point reached by the algorithm given above. In some sense then, this area is the most difficult area of the Mandelbrot set to plot, and indeed if you attempt to zoom into it you'll quickly discover that you have to increase the maximum number of iterations exponentially to avoid losing detail. I don't think it's really practical to explore this area very far without a Mandelbrot set program that can do SOI, such as Fractint, Fractal Witchcraft (which is a pain to get working now as it's a DOS program), almondbread (which is a pain to get working now as one of its prerequisites seems to no longer be available) or the program I'm writing.

If you could zoom a long way into it, what would the Siegel area of the Mandelbrot set look like? Well, we know the part inside the cardioid is inside the set. We also know it's not a bond point so at least part of the remainder is outside the set. The pattern of circles on the cardioid is self-similar (it's related to Farey addition) so the continent part would look similar to a shallow zoom on the same point. The rest is a bit hazy, though. My experiments suggest that the density and relative size of mu-molecules gets greater as one zooms in on this area, so an interesting question to ask is what the limit of this process is? Would the area outside the continent be covered with (possibly highly deformed) tesselating mu-molecules, leading to a black screen? Or would only part of the image be covered by mu-molecules, the rest being extremely dense filament structures? Both possibilities seem to be compatible with the local Hausdorff dimension being 2. Either way I suspect there are some seriously beautiful images hiding here if we can zoom in far enough.

Summary of the taxonomy:

• Exterior points (example: infinity)
• Member points
  • Interior points
    • Superstable points (examples: 0, -1)
    • Non-superstable points
  • Boundary points
    • Filaments
      • Repelling points
        • Terminal points (examples: -2, +/- i)
        • Branch points
      • Chaotic, irrational external angle filaments
    • Feignbaum points (example: ~-1.401155)
    • Shore points
      • Parabolic points
        • Cardioid cusps (examples, 0.25, -1.75)
        • Bond points (examples, -0.75, -1 +/- 0.25i, -1.25, 0.25 +/- 0.5i)
      • Siegel points (example: ~= -.39054087021840008+.58678790734696872i)

I'm writing a super-fast, super-accurate Mandelbrot set plotter. This one is unusual in that no explicit iteration limit is ever set - it will keep iterating any point until it either escapes or is determined to be in the set. This works pretty well and produces much sharper cusps than any other program I've seen.

However, the downside is that it can take an arbitrarily long to determine if any given point is in the set or not. Points on the edge of the set are (as always) the particularly tricky ones - some of them exhibit chaotic behavior for a very long time before becoming periodic or escaping.

Currently I use two methods for determining if a point is in the set. I have explicit formulae for the main cardioid and the period-2 bulb, and I look for periodic orbits using Floyd's algorithm.

I'm wondering if a third method might be practical - one that determines if a point is in the set without waiting for it to become periodic. Here's the method I'm considering.

Let the Mandelbrot iteration function f(z) = z² + c. Suppose we have a point c, and we suspect that it might be periodic with period N. We can determine if it is or not by finding a point z_f that is unchanged under N applications of f, i.e. a solution to f^N(z_f) = z_f. There will generally be 2^N such solutions, of which only N will be relevant to the case at hand. However, if we have a z that is likely to be near one of these fixed points we should be able to converge on one of them quickly by using Newton's method or something similar.

Once we have our fixed point, we take the derivative of f^N with respect to c and evaluate it at z_f. This is easy to do numerically - just evaluate f^N at z_f and a nearby point, take the difference and divide by the original distance. If the modulus of the result is 1 or less, then the attractor is stable and c is in the set. If it isn't, the point could still be in the set as we might have picked the wrong z_f or N - that is we are likely to get false negatives but no false positives using this method (which is exactly what we want - if we get a false negative we can
iterate some more, try a new N and/or a new z_f and eventually get the right answer).

Working through these steps symbolically gives us closed form solutions for the period-1 and period-2 attractors but if we try it for period 3 we get a polynomial of order 6 in z (2³ = 8 minus the two from factoring out the period-1 attractor) - I'm not sure if we can get a closed form solution for period-3 attractors from this.

The tricky part is choosing an appropriate N and an approximation to z_f for any given point. I think the way to do it is this: whenever we find a periodic point with Floyd's algorithm, iterate a further N times to figure out the period, and store the period and fixed point. When we start a new set of iterations on a point, check to see if any neighbouring points are periodic. If they are, try this method using the N from that point and the z_f (appropriately adjusted for the new c). When we find a periodic point this way, we store the found N and z_f just as for Floyd's algorithm.