On the Integer Spiral

The problem: find a procedure to map the integers 0, 1, 2 … uniquely into points (i, j) ∈ Z² with integer coordinates, starting with 0 ↦ (0, 0), and continuing in a counterclockwise spiral pattern. The first few points are (0, 0), (1, 0), (1, 1), (0, 1), (-1, 1), (-1, 0), (-1, -1), (0, -1), (1, -1), (2, -1), (2, 0), (2, 1), (2, 2), (2, 1), (2, 0) … The idea is to find a bijective mapping between N and Z². Of course there are simpler ones, but the integer spiral is particularly pretty:

An instant's doodling on graphing paper shows that every complete turn of the spiral completely surrounds the origin with a square, each having 1, 8, 16, 24 … new points. The first, or "zero-th" such "square" is the origin itself, having "side" 1. Each new square's boundary points are at a distance 1 from those of the older one (except for the corner), so each new square has side two plus the older square's side, or l = 2⋅n + 1. An appeal to the inclusion-exclusion principle proves that the perimeter of a square is four times the side minus the four corners counted twice; hence, the perimeter is comprised of p = 8⋅n points. The total area of this square is, obviously a = (2⋅n + 1)² points.

The i-th point lands on the n-th square's perimeter if and only if it is not inside it, that is, if it doesn't land on the (n-1)-th square's area. The index starts at zero, so this condition is equivalent to a.(n-1) ≤ i < a.n. Massaging the inequality to extract n:

   a.(n-1) ≤ i < a.n
≡ { Definition a }
   (2⋅(n-1) + 1)² ≤ i < (2⋅n + 1)²
≡
   (2⋅n - 1)² ≤ i < (2⋅n + 1)²
≡ { i is nonnegative }
   2⋅n - 1 ≤ √i < 2⋅n + 1
≡ { Algebra }
   n ≤ (√i + 1)/2 < n + 1
≡ { Definition floor }
   n = ⌊(√i + 1)/2⌋

Let d = i - (2⋅n - 1)² be how far around the n-th square's perimeter the index falls. There are four sides to it, each thus having 2⋅n points. If d is less than 2⋅n it is mapped to the first (i.e., left) side; if it is less than 4⋅n it's mapped to the second (top) one, and so on. Let this side index be k = ⌊d/(2⋅n)⌋.

As the spiral goes, each side has start- and end-points that depend on k like this:

k	Start	End
0	(n, 1-n)	(n, n)
1	(n-1, n)	(-n, n)
2	(-n, n-1)	(-n, -n)
3	(1-n, -n)	(n, n)

So, if sx, sy, dx, dy are the signs and offsets applied to each starting point (sx⋅(n - dx), sy⋅(n - dy)) as a function of k, we have:

k	dx	dy	sx	sy
0	0	1	+1	-1
1	1	0	+1	+1
2	0	1	-1	+1
3	1	0	-1	-1

Note that dy.k = dx.(k - 1) and sy.k = sx.(k - 1), modulo 4. It is not difficult to see that dx.k = k % 2, and sx.k = 1 - k % 4 + k % 2.

Aside: To explain it in full: first, sx.k must have period four, hence the reduction of the argument modulo 4. Second, it takes two consecutive values in pairs; this means that 0 and 1, and 2 and 3 must be treated the same, hence the reduction modulo 2. Finally, the result is shifted so that the range is ±1.

Some algebraic manipulation gives the alternative sx.k = 1 - 2⋅⌊(k % 4)/2⌋, which when implemented with bitwise operators is the very compact 1 - (k & 2).

Similarly, we need to find two characteristic functions, or multipliers, to apply to the corner the side offset, call it o = d % (2⋅n); so that points are traversed first up, then left, then down and finally down:

k	ox	oy
0	0	+1
1	-1	0
2	0	-1
3	+1	0

But this is simply -sx⋅dx (resp. -sy⋅dy)!. Expanding and collecting equal terms, we have that the i-th point has coordinates (sx⋅(n - dx⋅(o + 1)), sy⋅(n - dy⋅(o + 1))).

Aside: Had I chosen instead to subtract the offset from the side's endpoint, I would have arrived to this directly, without finding the need to simplify the expression.

Hence, it all translates directly to OCaml as:

let spiral i =
  let n  = truncate (0.5 *. (sqrt (float i) +. 1.)) in
  let d  = i - (2*n - 1) * (2*n - 1) in
  let k  = d / (2*n)
  and o  = d mod (2*n) in
  let sx = 1 -  k    land 2
  and sy = 1 - (k+3) land 2
  and dx =      k    land 1
  and dy =     (k+3) land 1 in
  (sx * (n - dx * (o + 1)), sy * (n - dy * (o + 1)))