Decomposing the Ulam spiral

The Ulam spiral is a well-known mystery to mathematicians: draw a 1 in a grid, 2 in the cell to the right, 3 above the 2, and so on in a spiral:

Colour the primes:

For reasons unknown, long diagonal chains of primes appear. In this very small Ulam spiral, we don’t really see them (19, 7, 23), but we will.

How can we draw our own spiral? The essential elements of the spiral are two-fold: a primality test, and winding the line of positive integers round and round the origin of the complex plane. Hm, composing functions sounds promising!

Our primality test is a function from the positive integers to (True, False). Our traversal/map is a function from the positive integers to C, the complex plane. That is, given an integer n it tells us where on the plane that point will be after traversing the spiral n steps.

Squeak images have Integer >> #isPrime built in, so we need just write the map function.

Integer >> ulam

    "Ulam maps integers to (the integral values of) the complex plane.

     1 maps to 0@0, 2 maps to 1@0, 3 maps to 1@1, and so on in an

     anticlockwise square spiral."

    | n lastCorner squareSize sqModulus bottomLeft topLeft topRight |

    n := self leastSatisfier: [:k | ((2 * k) + 1) squared].

    squareSize := (2*n) + 1.

    lastCorner := squareSize squared.

    bottomLeft := lastCorner - (1 * (squareSize - 1)).

    topLeft := lastCorner - (2 * (squareSize - 1)).

    topRight := lastCorner - (3 * (squareSize - 1)).

    sqModulus := squareSize // 2.

    self = lastCorner ifTrue: [^ sqModulus @ sqModulus negated].

    bottomLeft <= self ifTrue:

        ["bottom edge of square"

        ^ (sqModulus negated + (self - bottomLeft)) @ sqModulus negated].

    topLeft < self ifTrue:

        ["left edge of square"

        ^ sqModulus negated @ (sqModulus - (self - topLeft))].

    topRight < self ifTrue:

        ["top edge of square"        

        ^ (sqModulus - (self - topRight)) @ sqModulus].

    ^ sqModulus @ (sqModulus - (topRight - self)).



Integer >> leastSatisfier: aUnaryBlock

    "Can we find an integer value for which (aUnaryBlock value: n) = anInteger?"

    | k n max |

    "While this prevents nontermination, it means not being able to generate

     large spirals. Hold your nose and move along."

    max := 10000.

    n := 0.

    [k := aUnaryBlock value: n.

    k >= self ifTrue: [^ n].

    n = max ifTrue: [Error signal: 'Couldn''t satisfy equation ',

        aUnaryBlock decompile printString, ' with value ', self printString].

    n := n + 1] repeat.

With this in hand, the solution is simple: create a canvas, select the primes, and paint those points representing primes.

| form primes width |

width := 400.



primes := (1 to: width squared) select: #isPrime.

form := Form extent: width@width depth: 1.

primes do: [:p | | u |

    u := p ulam.

    "ulam returns origin-centred values, so we translate the points to the

     middle of the form. Too, we flip the points in the x-axis to match

     Form's layout, which as 1@1 at the top left."

    u := (u flipBy: #vertical centerAt: 0@0)

        translateBy: (width // 2)@(width // 2).



    "Filling always fills Rectangles. The smallest Rectangle is

     thus 1px wide. Let u = (x,y). Then u + 1 = (x + 1, y + 1)."

    form fill: (u corner: (u + 1)) fillColor: Color black].

PNGReadWriter putForm: form onFileNamed: 'ulam.png'

And result!

As a result of our insight – decomposing the problem into separate traversal, find and do-something functions, we may now construct different spirals: rectangular spirals, clockwise spirals, and so on.