Investigating The Properties Of A Square Spiral

At work the other day, I needed an algorithm that could visit an arbitrary number of neighbouring grid squares in a closest-neighbours-first pattern. The easiest way to do this was to enumerate a square spiral starting at the origin (0,0) and spiraling outwards.

I eventually decided on this Python procedure which exploits a quirk of complex arithmetic whereby a counter-clockwise turn of 𝛑/2 radians can be simply implemented by multiplication by the complex number 0+i. The other insight the code uses is that the length of each side of a spiral increases by 1 for every half-turn.

The outer loop of the function iterates across the turns of the spiral until a maximum number of points has been visited. The inner loop enumerates all the points on a single quarter-turn of the spiral.

def spiral(max):
 “””
     Yields an iteration through a counter-clockwise, 
     square-spiral starting at 0,0 and extending in 
     increments of 1 heading east initially and turning 
     counter-clockwise at each corner of the spiral.
 “”” 
 position = complex(0,0) # origin
 motion = complex(1,0)   # initially, head east
 turn = complex(0,1)     # always turn counter-clockwise
 extent = 0              # initially, zero extent
 
 count = max
 while count > 0:
    if motion.imag == 0:   # if about to head east or west...
      extent = extent + 1  # ...extend by one
    k = min(count,extent)
    count = count — k
    while k > 0:
      yield (position.real, position.imag)
      position = position + motion
      k = k — 1
    motion = motion * turn

Here’s the result of running plotting the output of spiral(625).

For some reason, I was curious as to how the average position of the points in the spiral would change as the number of points in the spiral increased. My initial intuition was that the average would converge to the origin, however, my intuition would prove to be wrong.

The following graph shows how the average position of all points previously plotted evolves:

There are some interesting points to note about this graph:

• it becomes smoother as the number of points in the spiral increases
• it seems to regularly pass through the fixed points (0,0), (0,1/2), (1/2, 0), (1/2, 1/2)
• the curve does seem to converge to a smooth curve, but what is that curve — it looks vaguely elliptical, but is it? (spoiler: it isn’t)

The fixed points in the 1/2 unit square centered on (1/4, 1/4) are interesting, and perhaps not too surprisingly correspond to points where the number of points in the curve has some simple relationship to a perfect square. Let N be the number of points in the spiral, then:

• if N=n² and n is odd, then the average is always (0,0)
• if N= n² and n is even, then the average is always (1/2,1/2)
• if N=(n+1)*n and n is odd, then the average is always (1/2, 0)
• if N=(n+1)*n and n is even, then the average is always (0, 1/2)

In the end, it is no surprise that the average doesn’t converge — perfect squares of odd numbers are pulling the average in one way, perfect squares of even numbers are pulling the average in a different way and there is a transition between the two extremes on each turn of the spiral. The curve of convergence is becoming smoother as the number of points in each quarter-turn of the spiral increases, and the general symmetry of the spiral is preventing arbitrary divergence from the observed curve.

Interestingly, the average of the averages does converge to (1/4, 1/4) as seen with this plot. Once we discover that parametric equation for the curve of convergence, it will become obvious why this convergence happens.

So What Is The Curve Of Convergence?

The curve of averages appears to be converging towards an ellipse-like curve. When I originally plotted the curve, it appeared very ellipse-like, but that was an artifact of scaling decisions made by the graphing software I was using.

How to calculate the curve of convergence?

We need an expression for each point of the curve as N, the number of points in the curve approaches infinity.

Let’s start with a spiral of N points, where N=n² and n is an odd number,

We know that the average for this spiral is 0. This is a consequence of all the points in the spiral being evenly distributed around the origin, so the contribution from each point is canceled out by a point exactly opposite to that point through the origin of the spiral.

Extending the spiral involves adding n points to the spiral starting at ((n+1)/2,-(n-1)/2) and extending to ((n+1)/2,(n-1)/2). These points can be described by the set: ((n+1)/2, -(n-1)/2+k-1) : for k in [1,n] which can be simplified as: ((n+1)/2, -(n+1)/2+k) : for k in [1,n].

The sum of all points of the first N+k points in the spiral is thus:

This simplifies to:

which further simplifies to:

The average for each in point the spiral is obtained by dividing these points by n²+k:

A quick sanity check, if we substitute k=0, the average is (0,0) which is the average at the starting position. If we substitute k=n, the average is (1/2,0), so we have preserved the expected results at the limits.

Now, we define a ratio r=k/n and rewrite the equations in terms of r and n.

Simplifying this we get:

Taking limits as n->∞:

Now, substituting x=r/2, we get:

Or, in the form of a quadratic equation:

Overlaying this over the plot of averages for finite spirals we get the orange curve segment, indicating that the equation is a good fit for the observed curve.

With symmetry arguments, you can derive the 3 equations for the other parts of the curve of convergence.

Plotting these curves yields:

Overlaid over the average of a spiral of finite length. we get:

The derivatives of the curve segments at (0,0) and (1/2,1/2) are -1 and at (1/2,0) and (0, 1/2) are +1 indicating that the complete curve is continuously differentiable and suggesting that it may be possible to derive a single equation for the complete curve, perhaps by rendering it into polar coordinates.

The parabolic equations can be expressed in this form:

which leads to the following 4 equivalent equations:

The focal points of these parabolas are: (1/4, 0), (1/4, 1/2), (0,1/4), (1/2, 1/4)

It is interesting to consider how close these jointed parabolas are to a pure circle centred on (1/4,1/4) that passes through each of the fixed points of the curve. That is a circle with this equation:

It turns out, it is quite close:

The curve of convergence is bounded on the outside by:

Or, shown together:

Approximating The Equation of the Curve

With a little bit of fiddling, I was able to find what appears to be a single equation of the curve that apparently fits 4 parabolic equations almost exactly.

The basic technique was to calculate the ratio of the interval between the parabola and the point (1/4,1/4) and √2/4 for 𝜃 in the range π/4 and 3π/4 and then plot this ratio against 𝜃. This resulted in this plot:

The next step was to fit the equation to the above curve:

The 4𝜃 is related to the 4-fold symmetry of the parabolic approximation and the fact that the ratio between the circle and one parabola cycles from 1 back to 1 again over the 90 degrees between 45 degrees and 135 degrees and does it 3 more times for each of the other parabolas.

The terms (3/(2√2)–1)/2 and (3/(2√2)+1)/2 are related to key dimensions derived from the parabola and the circle it approximates. The first is half the maximum difference of the ratio of the distances between the parabola and circle and the center of the circle along a radius of the circle. The second is the midpoint between these two ratios.

Folding the equation into the equation of a circle of radius 2√2/4 centered on (1+i)/4 and simplifying that results in the final equation.

Note that I have not yet derived this equation directly from the parabolic equations themselves. To assess possible errors between my guesstimate and the parabolic curves, I calculated the distance between the parabolic curves and the point (1/4,1/4) and also between the equation above and (1/4+1/4i) and then calculated the difference between the two values. These are plotted below. It is possible that these errors are due to a higher-order term that the equation above omits, but it may also be due to numerical errors however, I haven’t yet been able to rule either explanation out.

So, still to be done is to derive the final equation of the curve from the parabolic equations and so rule either in or out the presence of higher-order terms.

Another interesting thing to explore is constructing parabolic approximations with a larger number of parabolas. The key will likely be to ensure that the curve is continuously differentiable at the junction between each parabola.