Category Archives: Mathematics

Quintics Redux

A couple of years ago I wrote a blog post about finding the closest point on a cubic bezier to a given point in which I showed how to reduce the problem to solving a quintic equation and provided an implementation in C# that finds the roots of the quintic using Laguerre’s Method. I also mentioned at the end of that post that there is a fancier algorithm for solving specifically quintics described in Doyle and McMullen’s “Solving the Quintic by Iteration” [1] and speculate that it might be the fastest way to find the closest point to a bezier.

Recently I spent the time to actually implement the Doyle and McMullen algorithm in C++ as a numeric algorithm. I don’t know of any other implementations. My code is here.

It is an interesting algorithm. Its existence implies that although there is no general formula, akin to the quadratic formula, that you can plug into to get the roots of a 5th degree equation, and in fact although the roots of most 5th degree equations cannot even be represented as expressions composed of ordinary arithmetic operations plus radicals — you can represent numbers with radical expressions that are arbitrarily close to these roots. Each of these the roots can be defined exactly in terms of the limit a recursive function converges on as the depth of recursion approaches infinity; each step along the way approximates with increasing accuracy the root of an unsolvable quintic expressed using arithmetic and radicals.

Above I specify my implementation is “a numeric algorithm”, I mean as opposed to an implementation using a symbolic math package; I use ordinary floating point numbers. This is an important point: [1] is from the wilds of number theory, not from a CS text book, and it thus defines an algorithm over real numbers which of course have arbitrary precision. It was an open question to me whether the algorithm has value as a numeric algorithm. The issue I saw is that it is a solution to quintics in the single parameter Brioschi form. The Brioschi form, by magic, collapses every general quintic equation defined by six complex coefficients into a single complex number; it is impossible for such a reduction to not be limited by finite precision. That is, to some extent the finite precision of floating point numbers must determine to which general quintics a numeric implementation of [1] will return sensible results.

I will leave an analysis of this question to someone who is an actual mathematician or computer scientist but empirically the Doyle and McMullen algorithm does seem to me to have merit as a numeric algorithm at double precision. On randomly generated general quintics with coefficients uniformly distributed between -1000 and 1000, my implementation of solving the quintic by iteration returns results that look to me to be about as good or better than the implementation of Laguerre’s Method from Data Structures and Algorithms in C++ but is about 500 times faster. Both algorithms perform better in terms of speed and correctness when the coefficients are lower. When they are between -10 and 10 my implementation is more accurate and about 100 times faster than the Data Structures and Algorithms code.

My implementation works as follows:

  1. Given a general quintic, convert to principal form. I do this conversion as explained here. I found the resultant and solved for the c₁ and c₂ being canceled out using SageMath.
  2. Convert the principal quintic to Brioschi form. I followed the explanation here.
  3. Find two solutions to the Brioschi form quintic via iteration as described in [1] pages 32 to 33. The only difficulty here was that I needed to find the derivative to the function g(Z,w). I did this symbolically again via SageMath; however, a speed optimization would be to replace use of an explicit C++ function for g’ and instead evaluate g(Z,w) and g'(Z,w) simultaneously as described in the chapter in Numerical Recipes: The Art of Scientific Computing on polynomials. Also just as a note to anyone else who may want to write an implementation of this algorithm, the wikipedia article here has a mistake in the definition of h(Z,w). Use the original paper. (Or better yet cut-and-paste from Peter Doyle’s macsyma output and overload C++ such that the caret operator is exponentiation, which is what I did)
  4. Convert the two roots back to the general quintic.
  5. Test both roots. If the error is less than threshold k pass along both roots. If one is less than k pass along the one good root. If both roots yield errors greater than k perform n iterations of Halley’s Method and retest for one or two good roots. If neither root has error less than k pass both along anyway.
  6. If you have two good roots v₁ and v₂, perform synthetic division by (zv₁ )(zv₂ ) yielding a cubic and solve the cubic via radicals. If you have only one good root divide the quintic by (z-v) and solve the resulting quartic. I’m using the cubic solving procedure as described in Numerical Recipes and the quartic formula as described here.

My implementation is a header-only C++17 library (“quintic.hpp” in the github repo i linked to above) parametrized on the specific floating point type you want to use. Single precision is not good enough for this algorithm. Double precision works. I didn’t test on long doubles because Visual Studio does not support them.

“Triangular Life”

Recently I looked through a bunch of triangular cellular automata in which each (1) uses two states, (2) uses the simple alive cell count type of rule, and (3) uses the neighborhood around a cell c that is all the triangles that share a vertex with c; that is, the 12 shaded triangles below are the neighborhood around the yellow triangle:
12-cell triangular neighborhood
These cellular automata have state tables that can be thought of as 13 rows of 2 columns: there are 12 possible non-zero alive cell counts plus thee zero count and each of these counts can map to either alive or dead in the next generation depending on whether the cell in the current generation is alive or dead (column 1 or column 2). I looked at each of the 4096 cellular automata you get by filling the third through eighth rows of these state tables with each possible allocations of 0s and 1s and letting all other rows contain zeros.

A handful of these 4096 feature the spontaneous generation of gliders but one rule is clearly the triangular analog of Conway’s Life. I have no idea if this rule has been described before in the literature but it is the following:

On a triangular grid

  • If a cell is dead and it has exactly four or six vertex-adjacent alive neighbors then it alive in the next generation.
  • If a cell is alive and it has four to six vertex-adjacent alive neighbors, inclusive, then it remains alive in the next generation.
  • Otherwise it is dead in the next generation.

The above has a glider shown below that is often randomly generated and exhibits bounded growth.
Here it running in Jack Kutilek’s web-based CA player:

Tri Life gliders are slightly rarer than in Conway life because they are bigger in terms of number of alive cells in each glider “frame”. If you don’t see a glider in the above, stir it up by dragging in the window.

frames of the glider in cannonical triangular life


Over the past couple of weeks I wrote some code in C# to generate dissections of the rhombic triacontahedron into golden rhombohedrons. George Hart discusses these types of dissections here  and also talks about the problem of enumerating them in an appendix here — briefly, all this material by Hart and others is about how the fact that the rhombic triacontahedron and the rhombic enneacontahedron are zonohedra lead to both having interesting combinatoric properties which can be explored by coloring their dissections.

I was, however, more interested in how such dissections could be turned into an interlocking puzzle, akin to a traditional burr puzzle. amd as such needed code to generate 3D models of the dissections. My generation code is a dumb, constructive, brute force approach in which I just traverse the search space adding rhombohedrons to a candidate dissection in progress and backtracking when reaching a state in which it is impossible to add a rhombohedron without intersecting the one that was already added or the containing triacontahedron, keeping track of configurations that have already been explored.

Dissections of the rhombic triacontahedron into golden rhombohedrons (hereafter “blocks”) turns out to always need 10 and 10 of the two types of blocks that Hart refers to in the above as the “pointy” and “flat” varieties (and which I refer to as yellow and blue). Further it turns out that in all of these dissections there are four blocks that are completely internal, i.e. sharing no face with the triacontahedron; I also believe that the four internal blocks are always three blue and one yellow, but I’m not sure about that.

My strategy for finding an interlocking puzzle was the following:

  • Generate a bunch of raw dissections into blocks
  • For each dissection, search the adjacency graph for four pieces, the union of sets of five blocks, such that
    • Each piece forms a simple path in the dissection; that is, each block in the piece
      • is either an end block that is face adjacent to a next or previous block in the piece or is a non-end block that is face adjacent to a next block and a previous block.
      • and does not share any edges with other blocks in the piece except for the edges of the face adjacencies.
    • Each piece contains at least one fully internal block.
    • Each piece is “single axis disentangle-able” from each other piece, where we mean by that that there exists some edge e in the complete construction such that if given piece p1 and piece p2, if you offset p1 in the direction of  e by a small amount p1 does not intersect p2.
    • Each piece is not single axis disentangle-able from the union of the other three pieces.

I never managed to succeed in doing a complete enumeration, generating all of the dissections for reasons that I don’t feel like going into. (As I said above, I did not do anything fancy and it would be easier to just be smarter about how I do the generation than to make what I have more efficient; i.e. could have done the George Hart algorithm if I had known abouyt that or there are ways of transforming one dissection into another that I don’t do — I do an exhaustive search, period — but I never did the smarter stuff because I found what I was looking for, see below)

But from about 10 dissections I found one set of pieces that uniquely satisfies all of the above:

Here’s some video. (The individual blocks were 3D printed and super glued together)

I’m calling the above “rhombo”. Those pieces are rough because I only 3D printed the individual rhombohedrons and then superglued them together into the pieces, which is imprecise. I had to sand them heavily to get them to behave  nicely. I’ll eventually put full piece models up on Shapeways.

In the course of doing this work, it became apparent that there is no good computational geometry library for C# to use for something like this. There is one called Math.Net Numerics along with Math.Net Spatial that will get you vectors and matrices but not with all the convenience routines you’d expect to treat vectors like 3D points and so forth. What I ended up doing was extracting the vectors and matrices out of monogame and search-and-replacing “float” to “double” to get double precision. Here is that code on github. I also included in there 3D line segment/line segment intersection code and 3D triangle/triangle intersection code which I transliterated to C#. The line segment intersection code came from Paul Bourke’s web site. And the triangle intersection code came from running Tomas Moller’s C code through just a C preprocessor to resolve all the macros and then transliterating the result to C#.

Is there hexagonal analog of Conway’s Game of Life?

The short answer is that in the hexagonal case the best analog of Conway’s Game of Life — in my opinion as someone who has been a CA hobbyist for 30 years or so — is an original creation which I will describe for the first time in detail in this blog post. (Note: The links to cellular automata in this post go to a custom web-based player which may not run correctly on mobile devices.)

Regarding a hexagonal Game of Life, a key thing to understand is that Conway Life isn’t just the rule; it is the rule (Life 2333) running on a square grid with an 8-cell neighborhood i.e. the neighborhood of four squares directly adjacent to a square plus those diagonally adjacent. You can, of course, apply the same rule on a hexagonal grid using the natural six hexagon neighborhood but what you will get won’t look anything like Life. It will look like this: Life 2333 on a hexagon grid.

So if the above does not qualify as a hexagonal Game of Life then what would? Well, at the very least we need a glider. Carter Bays, a professor at the University of South Carolina, presented a Conway-like rule that admits a glider on the hexagonal lattice in 2005, Life 3,5/2 in his notation. However, by Bays’ own admission “this rule is not as rich as Conway’s Life” and indeed when we run Bays’ Hex Life on random input we do not see gliders: Life 3,5/2 on random input. The problem is that its glider is too big to occur randomly. Its glider is in a sense artificial. Part of the beauty of Conway Life is that gliders are frequently spontaneously generated. The other thing you’ll notice about Bays’ Life 3,5/2 is that it descends into still-lifes and oscillators too quickly. Conway Life is dynamic. It sprawls and grows, descends into bounded chaos, before finally decaying into still-lifes and oscillators.

To summarize, we want a hexagonal cellular automaton that

  1. Has a glider that is frequently generated by random initial input.
  2. Frequently exhibits bounded growth from random initial input.

It is my contention that there is no simple totalistic rule over two states and the hexagonal grid using the natural 6-cell neighborhood that exhibits both 1. and 2.

In order to achieve what we want, we need to drop one of the constraints. Using my cellular automata breeding application Lifelike, I explored dropping the constraint that the rule must be a simple totalistic rule over two states. What I have come up with is a cellular automaton that uses a simple totalistic rule over three states,  states 0 to 2. One way to view this move is to view the live cells in Conway Life as counters — beans, pennies, whatever — and imagine dropping the constraint that a cell can only contain one counter. Anyway, my rule is as follows:

  • Take the sum S of the 6-cell neighborhood where death=0, yellow=1, and red=2.
    • If the cell is currently dead it comes to life as yellow if S is exactly 4.
    • If the cell is yellow it goes to red if S is 1 to 4 inclusive or is exactly 6.
    • If the cell is red it stays red if S is 1 or 2 and goes to yellow if S is 4.
    • Otherwise it is dead in the next generation.

The above, which I call “Joe Life”, is as rich as Conway life: click this link to see it run. It exhibits bounded growth with about the same burn rate as Conway Life and features two gliders, the little fish and the big fish, which are frequently generated spontaneously.

The little fish

The big fish

I concede Conway’s rule is more elegant than Joe Life’s rule but if one thinks about it, Conway Life’s neighborhood is less natural than the six hexagon neighborhood in that it is kind of weird to include the the diagonally adjacent squares. So in my opinion, elegance that Joe Life lacks in its state table it makes up for in its neighborhood.

A Woven Icosahedron

You can weave an icosahedron from 10 strips of paper — or anyway you can weave a construction that has icosahedral geometry; it is actually more of a snub icosahedron.

From about 11″ strips of printer paper folded lengthwise to have two layers with a little overlap to lock into rings, as pictured below, locks well and is rigid:

Construction follows the pattern implied by the following: