The Integer Circle Algorithm

Consider the following pseudocode:

x <- 32;
y <- 0;
while(true){
  x <- x - floor((1/8) * y);
  y <- y + floor((1/8) * x);
  plotPointAt(x, y); 
}

x and y are always integers. This draws a sort of oval in a loop:

It’s a bit lumpy, though. One might try changing the 1/8ths to 1/16ths to draw a rounder circle — but they’ll get an octagon instead:

But by using a different initial value for x and y, say 64 and 0 respectively, we can get a rounder circle:

There’s an important point easily missed here: in the pseudocode above, we update x, and then the next line uses the new value of x to update y. This is key to the stability of the algorithm.

We can say that this integer-based circle-drawing algorithm is parameterized by four values, x0, y0, d (for “delta”), and e (for “epsilon”), as follows:

x <- x0;
y <- y0;
while(true){
  x <- x - floor(d * y);
  y <- y + floor(e * x);
  plotPointAt(x, y); 
}

Steve Witham discovered, and Corey Ziegler Hunts, Julian Ziegler Hunts, Bill Gosper, and Jack Holloway rediscovered that this “circle algorithm” can draw things that are distinctly noncircular. For instance, if x0=233 (this value of x0, which is also a Fibonacci number, is important in this case), y0=0, d=1, and e=1.38197...=4 sin(pi/5)^2, then this draws the following:

In this case, it takes 17070 steps for the (x, y) point to return to its initial position.

Why Floor?

This algorithm – in an earlier context where it was designed for drawing circles, d=e, and e was small – appears to have been discovered first on PDP-1 systems. (It looks like it was independently discovered by two separate groups and then rediscovered several times over the years, but the important part here is the PDP-1.)

The original versions of this algorithm used the PDP-1’s DIV instruction, which rounded towards negative infinity (i.e. it acted like the floor in the pseudocode above). Most widely used programming languages today (such as C++) instead implement integer division using rounding towards zero (although as usual, C didn’t specify how to round integer divisions until C99).

The PDP-1 was an unusual computer in different ways as well: its words of memory were 18 bits long (this is one example of why a char in C or C++ actually has an implementation-defined number of bits that is greater than or equal to 8). It also was the size of a room, cost $120,000, and was released in 1960, but that won’t be relevant to the rest of this article.

Here’s another image of a “circle” generated by a similar algorithm that removes this skew.

Reversibility

This integer circle algorithm is reversible. That means a few things:

  • We can run (x, y) points both forwards and backwards in steps.
  • Point paths never converge: there’s no way to have two points (x1, y1) and (x2, y2) such that running one step of the integer circle algorithm maps them to the same point (x3, y3).
  • In particular, if a point’s path ever loops, it must return to its starting point.

To see this, we can write the inner loop of the integer circle algorithm as follows:

x[n+1] = x[n] - floor(d y[n])
y[n+1] = y[n] + floor(e x[n+1])

We then have

y[n] = y[n+1] - floor(e x[n+1])
x[n] = x[n+1] + floor(d y[n])

which gives an iteration for obtaining x[n] and y[n] from x[n+1] and y[n+1].

Un-Floored Circle Algorithm Behavior

Before we get to the main portion of this article – about the Mandelbrot set-like images that can be generated using this algorithm – here’s some analysis of the theoretical behavior of this algorithm without floors, and using real numbers instead of integers.

In this case, we can solve the resulting system of recurrence relations exactly (where x[0], y[0], d, and e are specified):

x[n+1] = x[n] - d y[n]
y[n+1] = y[n] + e x[n+1]

Substituting the first equation into the second gives the system

x[n+1] = x[n] - d y[n]
y[n+1] = y[n] + e x[n] - d e y[n]

This can be expressed as a matrix recurrence:

( x[n+1] )   ( 1, -d   ) ( x[n] )
( y[n+1] ) = ( e, 1-de ) ( y[n] )

This gives the nth coordinate in terms of x[0] and y[0] using a matrix power. (Also note that this matrix has a determinant of 1.)

( x[n] )   ( 1,    d )^n ( x[0] )
( y[n] ) = ( e, 1-de )   ( y[0] )

To find a closed form for this matrix power, we diagonalize it:

( 1,    d )^n
( e, 1-de )   = P D^n P^(-1)

where

    ( d(1/2 - a), d(1/2 + a) )
P = (          1,          1 ),

    ( 1 - de(1/2 + a), 0               )
D = (               0, 1 - de(1/2 - a) ),

and

a = sqrt(1/4 - 1/(de)).

where I’ve made some assumptions in simplifying that d>0 and e>0.

If d*e <= 4, then a is pure imaginary, so the diagonal elements of D have norm (1 - de/2)^2 - (de)^2(1/4 - 1/(de)) = 1, and so this algorithm draws an ellipse (with some amount of skew). If d*e > 4, though, point trajectories will follow a different conic section and be unbounded.

Many of these observations also apply to the integer circle algorithm – but the additional floors often throw a wrench into these conclusions, as we’ve seen from the d=1, e=4 sqrt(pi/5)^2 case above. There can also be cases where d*e < 4 but the trajectory explodes, and cases where d*e > 4 but the trajectory is bounded for the integer circle algorithm. For more information, please see Ziegler Hunts, Ziegler Hunts, Gosper, and Holloway’s book.

Ziegler Hunts, Ziegler Hunts, Gosper, and Holloway also provide an equivalent expression for this recurrence without complex numbers:

x[n] = x[0] cos(n w) + (x[0]/2 - y[0]/e) sin(n w)/b
y[n] = y[0] cos(n w) + (x[0]/d - y[0]/2) sin(n w)/b

where

w = 2 arcsin(sqrt(d e)/2)

and

b = sqrt(1/(d e) - 1/4).

From this, the four authors define the period, the (real) value of n > 0 at which the trajectory first returns to its starting point in this system, p = pi/arcsin(sqrt(d e)/2). This happens to correspond to approximate rotational symmetry in the plots below. (The inverse of this equation is d * e = 4 sin(pi/p)^2 – this is where the value of e in the above plot came from.

In 1990, Neal and Pitteway described (among other algorithms) a three-step circle algorithm that removes the skew in the two-step circle algorithm above in their article “Yet More Circle Generators” (The Computer Journal, Volume 33, Issue 5, 1990, Pages 408-411, https://doi.org/10.1093/comjnl/33.5.408):

while(true){
  x <- x - (p/2) y;
  y <- y + p x;
  x <- x - (p/2) y;
  plotPointAt(x, y);
}

This draws an ellipse with major and minor axes parallel to the coordinate axes, and is equivalent to an un-skewed version of the two-step integer circle algorithm. (Thanks to https://blog.hrvoje.org/2020/05/drawing-circles/ for the tip.)

Ziegler Hunts, Ziegler Hunts, Gosper, and Holloway introduce equivalents of the d and e parameters into this algorithm and use it for many of their illustrations:

while(true){
  x <- x - (a/2) y;
  y <- y + b x;
  x <- x - (a/2) y;
  plotPointAt(x, y);
}

Rugplots: Visualizing Choices of Parameters in the Integer Circle Algorithm

As described above, the (two-step) integer circle algorithm is parameterized by four values: x0, y0, d, and e. How can we visualize how the behavior of the algorithm – which at time can draw unusual shapes – changes as we modify these parameters?

Much like how many visualizations of the Mandelbrot set count and plot the number of iterations it takes for each complex starting point to escape, Ziegler Hunts, Ziegler Hunts, Gosper, and Holloway choose to associate each (x0, y0, d, e) tuplet with the number of iterations it takes for the initial point (x0, y0) to return to its starting point.

They then plot two-dimensional slices of this four-dimensional space – such as plotting the number of iterations varies with x and y, holding d and e constant.

The results are quite striking: here’s the first XY plot they plotted, using d=1 and e = 4 sin(pi/9)^2 (i.e. d=1 and p=9):

Period 9

Here’s d = 1, e = 0.91939:

Readers may notice something at this point...

And here’s a version of the d=1, p=5 system we originally presented:

Since all (x, y) points on a trajectory have the same period, individual trajectories can be easily seen on these rugplots. This also provides a way to optimize rugplot rendering – once one trajectory has been calculated, all other pixels along that trajectory can be filled with the same value. However, modern graphics processors are also fast enough that they can visualize rugplots in real time – see this Shadertoy for an example.

Mapping Periods to Colors

Such a visualization algorithm needs a way to turn a period – the number of iterations the algorithm takes before it loops – to a color. These methods have a large impact on structure visibility. A previous version of this article presented the mapping (with some modification, presented as pseudocode):

r = 65536 * log(period)/log(maxPeriod);
return { (r % 1), (r % 256)/256, r }; 

where period is the period of a pixel, maxPeriod is the maximum period in the image, % is the modulus operator, and this returns sRGB colors.

For d=100, e=1/99, this generates an image like this:

Made with a slightly modified version

Ziegler Hunts, Ziegler Hunts, Gosper, and Holloway use a different method:

L = log(max(maxPeriod, 2));
d = max(6, ceiling(log(8, maxPeriod)));
r = (8^d - 1)*(1 - log(min(period, maxPeriod)))/L;
Express r as a 3*d-bit number ... b_2 b_1 b_0 in base 2, where b_0 is the least significant bit;
Deinterleave these bits to get three d-bit numbers:
R = ... b_8 b_5 b_2;
G = ... b_7 b_4 b_1;
B = ... b_6 b_3 b_0;
return { R/2^d, G/2^d, B/2^d}

where period is the period of a pixel, maxPeriod is the maximum expected period, and the output is once again in sRGB. This gives the following image:

D-E Plots

Instead of holding d and e constant and varying x0 and y0, we can instead hold x0 and y0 constant and vary d and e, visualizing a different plane in this four-dimensional space.

Here’s what it looks like if we vary d and e, using the initial conditions x0=1, y0=0:

Rather large (4001×4001), so viewing in hi-res is recommended. 1px lines are often in fact infinitesimally small and due to scan granularity.

Gray areas denote areas where the period was greater than 5000 (or the system never looped), and white areas have period 1.

The hyperbola d*e = 4 is visible – but we can see areas (such as on the upper-right) where (d, e) pairs beyond this limit nevertheless looped. As mentioned above, there are also areas inside this hyperbola that never loop!

The rectangular structure visible here may be surprising. A given (x, y) trajectory is bounded, and only takes integer values; this means that the trajectory (and thus period) doesn’t depend on small changes in d and e in these cases, creating rectangular blocks – and this knowledge can be used to optimize rendering of these types of slices.

Here are some close-ups of this structure:

Generated by Julian Ziegler Hunts.

Nevertheless, we can see accumulation points in the above diagram where trajectories get arbitrarily close to diverging, and periods approach infinity. One might also ask what happens exactly on the boundaries between two of these rectangular regions, where the period changes from one block to another. We can find the boundaries of these regions exactly (they’re rational numbers) – and it turns out they often hold more structure and limit points. Here’s a plot by Bill Gosper of a region with several points with large periods denoted.

Double-click to see in high resolution
Double-click to see in high resolution

There’s tons more to study of course, but in reality all of the pictures we see of these plots, whether x-y, d-e, or some combination of the two, are all projections of a four-dimensional, jagged, infinite object. Here are a few of its properties, and conjectures about this space:

  • Most of it is contained within the region d * e <= 4 – but not all of it!
  • There appear to be no closed trajectories inside the region d * e < 0.
  • We can talk about congruence regions, which are regions where the trajectories are the same (ignoring translation in the x and y axes). These roughly correspond to the rectangles in the d-e plots above. When sliced through with a 3-dimensional plane, they look like polyhedra – but some of their sides are curved!
  • According to Ziegler Hunts, Ziegler Hunts, Gosper, and Holloway, these congruence regions are Cartesian products of two polygons with hyperbolically curved sides. They conjecture that these polygons have at most 6 sides.
  • There can be infinitely fine detail around points for which d * e is irrational – but there’s also a conjecture that there exists no (x, y, d, e) tuplet such that d * e is irrational, 0 < d * e <= 4, and the trajectory never repeats.

I’ve hosted a version of the (now somewhat old) “Problems posed and solved” list at neilbickford.com/ica/problems.txt.

It’s interesting to look back and see how much mathematics can be produced from an algorithm that uses only a few lines. We still don’t know everything about this space, and there are still many problems waiting to be uncovered and solved. In addition, I hope readers will find trying to visualize this space fun and interesting.

Why This Article Was Renamed (Content Warning)

Please note that the content of this section is relatively intense, and is unrelated to the rest of the article above; it’s fine with me to republish this text without this section.

This article was originally published in April 2011 under the name “The Minsky Circle Algorithm”. So far as I can tell, the earliest version of this algorithm (described at the top of this post) was discovered independently by David Mapes (see this interview here; thank you to this page for the link) and Marvin Minsky. Minsky published this in HAKMEM, which became a popular document and wound up being the first place I saw this algorithm. (Ziegler Hunts, Ziegler Hunts, Gosper, and Holloway’s book Minsky’s and Trinsky’s was published in December 2010, and almost all of this article starting from the second code sample was originally sourced from there and did not involve Minsky to the best of my knowledge.)

But revelations about Minsky’s past wound up being rather horrific. As described with sources here, Minsky had affiliations with Jeffrey Epstein, the infamous convicted sex trafficking ring operator; he held two conferences on Epstein’s private island, one (in 2011) after Epstein was a registered sex offender (he was first convicted in 2008), and received the first, $100,000 donation from Epstein to MIT. (The MIT Media Lab was largely responsible for courting $850,000 of donations from Epstein and a total of $7.5 million from donors acting on Epstein’s guidance, according to Business Insider. Minsky founded the MIT AI Lab, which merged in 2003 to become MIT’s CSAIL). In addition, in a 2015 deposition, Virginia Giuffre testified that Epstein associate Ghislaine Maxwell directed her to have sex with him, which is horrific. Minsky’s widow disputed this in an interview to the New York Post (but it should be noted that a deposition and a news article are rather different). This horrified me when I heard about it, and it didn’t sit well with me that I had an article that helped to prop up Minsky’s name.

I think this algorithm is interesting on its own, but it’s difficult to say the least to separate a person from the algorithm when the algorithm, related algorithms, and the related mathematical objects are named after the person.

I switched the article to private mode, and reworked the article to make it possible to consider the algorithm (and the results, which form most of this article and are largely not due to Minsky) without considering Minsky by removing reference to him (except in this note).

I also took the opportunity to rewrite much of the article, redoing the introduction with new static images, adding a more thorough derivation of properties of the un-floored circle algorithm, including the un-skewed algorithm (the previous version called this the “Trinsky algorithm”, but I now know it was published by Neal and Pitteway in 1993), adding a better description of the four authors’ coloring method, fixing a number of errors in the last section, and fixing the link to the “Problems posed and solved” article (which broke after I moved hosting servers; thank you to the Internet Archive WayBack Machine). My style and formatting has also changed a lot in the last 9 years, and so I changed some things I thought were cool 9 years ago but have since reconsidered. I think the article is much better as a result, but this took longer than I expected. I probably should have changed it earlier, really.

6 thoughts on “The Integer Circle Algorithm

  1. As a note on the coloring methods above, Dr. Fibdork uses the Processing code
    double LogMax=log(maximum);
    for(int i=0;i<xSize;i++){
    for(int j=0;j<ySize;j++){
    result=(int)(16777216*(log(periods[i][j])/LogMax)); //Log transformation
    fill(result%256,((result/256)%256),result/65536); //Choose the color
    rect(i,j,1,1); //Draw a pixel with that color
    }
    }

    where maximum is the maximum value encountered, periods is a matrix containing the periods of the points, and xSize and ySize are the sizes of the window.
    This coloring is pretty fast, but leaves out many small features. However, Corey and Julian use the Mathematica code
    (Returns a function for converting periods into colors. maxperiod is the largest period in the list of periods the function is to be applied to.)
    ColorAlgorithm[maxperiod_]/;And[IntegerQ[maxperiod],maxperiod>0]:=
    Block[{d,L},
    FromDigits[Partition[Rest[IntegerDigits[8^d+Floor[(8^d-1)*(1-Log[Min[maxperiod,#]]/L)],2]],3],2]/2.^d&/.{L->Log[N[Max[maxperiod,2]]],d->Max[6,Ceiling[Log[8,maxperiod]]]}
    ]

    They then apply this function to the array of periods, and then plot the colors they get. I don’t know if either is the best coloring method possible, so I would encourage others to come up with their own.
    –Neil Bickford

  2. That is a fascinating post. Yesterday at U of Arkansas, I gave a talk encouraging people to use computers to look at obscure problems. The computers of today are a million times as powerful as computers from a few decades ago. Thus, if a problem vulnerable to computer study hasn’t been looked at recently, it basically hasn’t been looked at.
    Stage 2 in a modern computer analysis is to generate some big numbers or number sequences that come out of the computer searches, and to post results online. Okay, that is definitely done.
    These are some beautiful results here.

  3. Wonderful write up! (Just one minor request: please put a pause button on the munching squares animation!) Thanks!

  4. Several decades ago I spent a good bit of time performing computer experiments with 2D cellular automata. I went on a long journey of testing various grid geometries, systems of phased updates and rules modeled on chemical reactions.

    Also, coincidentally, I also spent a good bit of time working on very close relatives of the Minsky algorithm you’ve discussed above, though in the context of physics simulations rather than a study of the detailed patterns one can observe.

    There is no end to the experimental fun of that kind of activity. But does it all lead anywhere? The combination of empirical observations interspersed with some mathematical deductions now makes me a little bit jittery, because it seems to spread out, both geometrically and intellectually, in all directions without any horizon or sight of land.

    In the not so distant past, number theorists spent disturbingly long weeks and months creating lists of primes, catalogs of polynomials and so on. Presumably they didn’t embark on such time-consuming studies without some idea about where they hoped to arrive. Now that we have these marvelous computers, the possibilities are (countably) infinite, but the goals are off somewhere in the mist.

    This is decidedly not meant as a critique. It is just my idea of occasionally putting up a periscope to look around and see where we might be heading. Any comments would be appreciated.

    • The inexhaustibility of mathematics going off in all directions can be exhilarating but also disheartening. I was reading about ways to add/multiply 1s to build up any integer. So many variations can be done, using other digits, operations, +/-, ^, mod, binary, bases, etc…, in how many ways, min number of symbols, of 1s, of parentheses, …
      then I read
      http://www.maa.org/editorial/mathgames/mathgames_04_12_04.html
      and I just felt overwhelmed!

  5. Pingback: TIL #262: Minsky Circle Algorithm – Jason Neal's Blog

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

This site uses Akismet to reduce spam. Learn how your comment data is processed.