I just ran across this lovely page about guillochés, and something in it bothered me slightly. I reflected on that feeling, and I realized what it was: I felt like the complexity of using transcendental functions was somehow too large for the spare and parsimonious loveliness that is a guilloché. That is, no reflection on Max Irwin, but we ought to be able to generate these designs on a computer without invoking floating-point sums of Tschebysheff polynomials or table lookups, because of course the natural world does not work in that way.

After some thought, this led to the JS implementation of guillochés as linear dynamical systems that you see in the gray canvas here.

However, this system is much more general than guillochés. You can explore its 12-dimensional space by randomizing the parameters repeatedly and mousing over the canvas; if it stops being interesting, you can reset it to return to the original guilloché.

My initial uneasiness led me to think about Minsky’s circle algorithm in HAKMEM item 149, (cf. AIM-239 item 149 on “page 73”, the 74th page of 107) which approximates sines and cosines discretely using a simple linear dynamical system:

NEW X = OLD X - ε * OLD YWhich is to say,

NEW Y = OLD Y + ε * NEW(!) X

1 | -ε | ||

ε | 1-ε² |

It seems to me
that we ought to be able to generalize this system
to produce guillochés,
simply by changing the transition matrix
and using a state vector
with more dimensions,
of which the first two will be used for plotting,
and the others are “hidden variables”
that merely affect the future states of the system.
For example,
the matrix

gives us two independent Minsky oscillators
the sum of whose values
produces the next point to plot,
giving us
what Irwin calls
“the basic rosette guilloché”.

0 | 0 | 1 | 0 | 1 | 0 | ||

0 | 0 | 0 | 1 | 0 | 1 | ||

0 | 0 | 1 | -ε₀ | 0 | 0 | ||

0 | 0 | ε₀ | 1-ε₀² | 0 | 0 | ||

0 | 0 | 0 | 0 | 1 | -ε₁ | ||

0 | 0 | 0 | 0 | ε₁ | 1-ε₁² |

There are a couple of interesting things about this matrix. First, it’s nearly diagonal: a third of its diagonal elements are 1, and another third are nearly 1, while of its 30 off-diagonal elements, 22 are 0 and four more are nearly 0. (This results in its spatial behavior being fairly coherent, since the state-variable changes from one step to the next are small, and its causality graph being fairly simple.) Second, as a result, it is a fairly sparse matrix. Third, it is a function of only two variables, ε₀ and ε₁, (and nearly a linear one); these two represent the angular velocities of its two oscillators. The space of basic guilloché rosettes up to uniform scaling and translation is two-dimensional, but these two variables do not represent its two dimensions; rather, their ratio represents one of the two dimensions, a common factor between them is a step size parameter controlling the granularity of the approximation (which may result in chaotic behavior!) while the relative radii or amplitudes of the oscillators is determined by their initial state. And it is singular, because its rank is only 4, which is to say that it collapses a six-dimensional input space into four dimensions (because its first two columns are entirely zero). But within those four dimensions (that is, if we amputate the first two rows and columns) it is unitary.

These curious properties imply that we are not likely to stumble across a matrix like this if we simply randomly explore the 36-dimensional space of 6×6 matrices, because almost all matrices are invertible and almost no invertible matrices are unitary. (I hope the reader will forgive me if this is too basic; my grasp of linear algebra is fairly weak. Conversely, I hope the reader will forgive me if this is too complicated, and perhaps take a detour through the basics of linear algebra, which is super fun!) But, if we want to explore similar systems, we could perhaps decompose it into a nearly-diagonal N×N hidden-state-transition submatrix and an arbitrary 2×N projection into pixel space. In this way, the original Minsky circle algorithm becomes capable of drawing arbitrarily scaled and rotated ellipses while its state space remains two-dimensional.

A slightly different generalization of the circle algorithm
which assigns a lower Kolmogorov complexity
to the “basic rosette guilloché” variant above
is a sequence of instructions of the form
*xᵢ* ← *xᵢ* ± *xⱼ* `>>` *k*,
where *i* ≠ *j*,
followed by a selection of disjoint subsets of the *xᵢ*
for the pixel coordinates.
Among other advantages,
this constrains the system to be invertible
(as long as (*a*+*b*)-*b* = *a* is valid,
which is true for integer math but not floating point),
since this basic instruction preserves *xⱼ*,
and converts a uniform distribution on *k*
into a roughly exponential distribution on ε.
It has the disadvantage that
To represent this as a series of numeric (*i*, *j*, *k*) tuples,
we can arbitrarily choose to represent the ± as the sign of *k*,
so the original Minsky circle system is simply the sequence
[(0, 1, -lg ε), (1, 0, lg ε)].

It happens that this particular sequence of instructions
has an exhaustive set of (*i*, *j*) pairs
lexicographically sorted,
as does the “basic guilloché” matrix above,
and if we restrict our attention to sequences that share this property,
we can elide the (*i*, *j*) pairs entirely,
leaving us with only the *k*s,
which additionally need only 4 to 7 bits each.
(If *k* is 5 bits, its minimum value is -32,
and we can take -*xⱼ* `>>` 32
to be 0
if our state variables are 32-bit integers;
but 16-bit integers may be adequately interesting.)

This reduces the representation of the underlying dynamical system of Minsky’s circles to some 8 to 14 bits, of which some 2⁴ to 2⁷ instances are instances of Minsky’s circles, which is small enough that you could imagine finding the Minsky behavior either by exhaustive enumeration or simply by search. The “basic guilloché” 4×4 matrix, then, becomes 4×3×[4 to 7] = [48 to 84] bits.

We should expect the constraint to invertibility to result in a more interesting state space because it excludes systems with point attractors — the system can have fixed points, but no other states can feed into them. This should become less important as the dimensionality increases, but I expect many of the more interesting systems to appear at low dimensionalities because our minds and especially visual cortices are not equipped to infer large numbers of hidden variables from projections into two dimensions.

If we were to allow the *k*s to take arbitrary real values
instead of simple integers,
then we could explore continuous variation over this state space,
but because the orbits are generally going to be chaotic,
this will not result in continuous variation of the image.

Systems in this family are eminently suited to simple digital hardware implementation, and audio output from them should also be interesting, since it should be periodic and produce a variety of somewhat distorted sinewaves. Continuous-time analog implementations may also be possible, for example, with each state variable represented by a capacitor and each non-diagonal matrix element implemented by the gain on an op-amp, or in current-mode with each state variable represented by an inductor. The correspondence between such analog implementations and the discrete-time algorithm is necessarily imperfect for a variety of reasons, but both may be interesting systems.

Now, therefore, be it resolved that here I will write code to systematically explore and visualize this family of low-Kolmogorov-complexity dynamical systems.

I posted this comment on Mark VandeWettering’s post about the algorithm:

With fixed-point (or generally non-arbitrary-precision) arithmetic, this algorithm is a deterministic finite-state automaton with no input and, say, 16 bits of state in X and 16 bits of state in Y. Any DFA with no input has eventually-periodic behavior, perhaps with a degenerate period-1 single-state loop, but it doesn't necessarily follow that its initial state is part of whatever periodic attractor it ends up in; it could be a sequence like A, B, C, B, (repeat). But for this to be true, some state must have more than one predecessor — in this example, B is preceded by both A and C. This means that the DFA's transition function does not have a deterministic inverse: if you construct a time-reversed version of it, it's nondeterministic. If the transition function is reversible (information-preserving), then you can run the DFA backward to the initial state (and beyond!) from any eventual successor state.But the state transition function implemented by the HAKMEM-149 algorithm is reversible, as noted in HAKMEM-150, because each of the two steps in the loop, “x -= ε·y” and “y += ε·x”, are individually invertible — in fixed-point arithmetic, anyway. (In floating point, maybe you lose low-precision bits if the change results in floating the point one binary place to the left, so this doesn't apply. And of course in “modern C” signed overflow is undefined, because to compiler vendors, the correctness of your code is less important than fractional percentage point performance increases.)

Therefore, with fixed-point arithmetic with the traditional overflow semantics, every possible initial state is a point on its own attractor — we eventually return to the original point, QED.

(However, as so vividly illustrated in “Minskys and Trinskys”, the periodic behavior from roundoff error and overflow need not look anything like a circle.)

In floating-point, the algorithm is not necessarily reversible, but the above argument demonstrates that it is still necessarily periodic — although conceivably with periods of up to 2⁶⁴ or 2¹²⁸. In practice, it converges to periodicity very rapidly.

This reversibility argument applies equally well to the
second and third
generalizations of Minsky’s algorithm
described above;
as long as the instruction
that computes the amount by which to change *xᵢ*
does not depend on the previous value of *xᵢ*,
it is guaranteed to be reversible and periodic
and therefore eventually return to its starting point.
However,
you don’t need to increase the dimensionality much
before the potential periodicity time
becomes large compared to the lifetime of the universe.

Babbage designed his Difference Engine
to tabulate polynomials
using the Method of Finite Differences,
a sort of generalization of Pascal’s Triangle
sometimes called the Method of Vanishing Triangles
or the Method of Difference Tables.
This is a very simple way to extrapolate a polynomial sequence.
First, for a degree-N sequence, you write down N+1 terms of the sequence.
Suppose we want to tabulate *y* = *x*² - 2x + 1;
first we write down its values for *x* = 0, 1, 2,
and just for fun, 3 and 4:

1 0 1 4 9Now we calculate the first forward differences Δ

These are the first differences Δ

1 0 1 4 9 1 -1 1 3 5

And Δ³

1 0 1 4 9 1 -1 1 3 5 1 -2 2 2 2

At this point our finite differences have petered out in an infinite sequence of zeroes; Δ³

1 0 1 4 9 1 -1 1 3 5 1 -2 2 2 2 1 -3 4 0 0

This difference table is also known as a “vanishing triangle” because (if you ignore the gray elements on the left of each row that depend on things before the beginning of the sequence, which are normally not considered part of it) it has a triangular shape, and if the input sequence is polynomial, sooner or later you get to a row where the triangle is all zeroes. Then the leftmost diagonal column of this triangle down to where it vanishes — here, (1, -1, 2) — provides enough information to reconstruct the whole original sequence.

And the relationship between each such diagonal
and the diagonal to its right is linear,
so it can be written as a matrix multiplication,
and in fact a very simple one
that can be computed without scalar multiplication:

=

, and then

=

,
and by multiplying by that same matrix again and again
we can get further state vectors
whose first element is the value of that polynomial
at some integer.

1 | 1 | 0 | ||

0 | 1 | 1 | ||

0 | 0 | 1 |

1 | ||

-1 | ||

2 |

0 | ||

1 | ||

2 |

1 | 1 | 0 | ||

0 | 1 | 1 | ||

0 | 0 | 1 |

0 | ||

1 | ||

2 |

1 | ||

3 | ||

2 |

And that same matrix generates not just *this*
quadratic polynomial sequence,
but *every* quadratic polynomial sequence,
with the coefficients and starting *x*
being encoded in the initial 3-vector.
And of course the relationships between those first three terms,
between the coefficients and the terms,
the coefficients and the differences,
and the differences and those first three grayed-out terms
that would be part of Δ³*y*
but are normally not considered part of it,
are all invertible linear transformations as well,
which means that you can put any combination of them
into an equivalent transition matrix.

What this means
is that our first and second generalizations of HAKMEM-149
it is capable of calculating arbitrary
degree-*n* polynomial trajectories
with *n*+1 state variables,
although clearly they is not limited to that — they lack that constraint,
and indeed *one single matrix*
out of all the possible matrices they can implement
generates all those polynomial trajectories
under different initial conditions!
This might seem to conflict with the earlier claims of periodicity
made for the second and third generalizations,
since we don’t normally think of polynomials as periodic,
but remember that those claims
depend on the finite size of the state variables,
enforced by overflow.