# Generating an Ulam spiral

An Ulam Spiral is quite an interesting construction, revealing unexpected features in the distribution of primes.

Here is a related topic with one answer by Pinguin Dirk, who has provided one approach.

I've created a faster implementation, but I believe it can surely be speeded up:

ClearAll[spiral]
spiral[n_?OddQ] := Nest[
With[{d = Length@#, l = #[[-1, -1]]},
Composition[
Insert[#, l + 3 d + 2 + Range[d + 2], -1] &,
Insert[#\[Transpose], l + 2 d + 1 + Range[d + 1], 1]\[Transpose] &,
Insert[#, l + d + Range[d + 1, 1, -1], 1] &,
Insert[#\[Transpose], l + Range[d, 1, -1], -1]\[Transpose] &
][#]] &,
{{1}},
(n - 1)/2];

spiral[5] // MatrixForm
spiral[201] // PrimeQ // Boole // ArrayPlot


This is good enough for my purposes, and I don't have time to focus on this now, However, it might be interesting for future visitors, so the question is how to make it faster by improving/rewriting.

Note: There may be confusion concerning whether I need a matrix populated with primes or one populated only with 0 and 1, 1 indicating the prime positions. Let's assume it doesn't matter. Maybe someone has a neat idea for the first representation that does not look good in the second one, so I'm leaving it open.

• Question: is there any concrete use for that spiral? Or is it purely recreational?
– Öskå
Commented Jun 9, 2014 at 14:23
• @Öskå Well, it will be a part of my final project for mma course :) So now I have to focus on the rest while I'm curious what folks can show here :)
– Kuba
Commented Jun 9, 2014 at 14:24
• Well, university is purely recreational anyway ;o)
– Öskå
Commented Jun 9, 2014 at 14:27
• @eldo The questions is not how but how faster :) but thanks for the link to demos I've not seen.
– Kuba
Commented Jun 9, 2014 at 15:47
• Related: (8209264) Commented Jul 12, 2015 at 2:53

The function findPrimePosInBoundarys below finds out which primes are on which square and where they are on such squares. The coordinate is just an integer and for layer = 3, the coordinates indicate the following positions.

It does so for each layer<=layers

findPrimePosInBoundarys[layers_] :=
Block[
{intRoot, primePiClean, primes, ran, squares, oneMostSquares,
primePis, primeSpans, primesPerLayer}
,
intRoot = 2*layers - 1;
primePiClean = PrimePi[intRoot^2];
primes = Prime /@ Range[primePiClean];
ran = Range[3, intRoot, 2];
squares = ran^2;
oneMostSquares = Prepend[Most@squares, 1];
primePis = PrimePi /@ squares;
primeSpans =
MapThread[Span, {Most[Prepend[primePis, 0] + 1], primePis}];
primesPerLayer = Part[primes, #] & /@ primeSpans;
primesPerLayer - oneMostSquares
]


Below is a compiled function. Given the positions on the square, and which layer we are working on, it transforms this integer position into a {x,y} coordinate.

cfu =
Compile[
{{primePos, _Integer, 1}, {layerMO, _Integer, 0}}
,
Block[
{qRs, varCo, fixedCo, case, pos, len}
,
qRs = {Quotient[#, 2 layerMO, 1], Mod[#, 2 layerMO, 1]} & /@
primePos;
len = Length@primePos;
Table[
{case, pos} = qRs[[i]];

varCo = layerMO - pos;
If[
case == 0 || case == 3
,
varCo = -varCo;
];

If[
case == 0 || case == 1
,
fixedCo = layerMO;
,
fixedCo = -layerMO
];

If[
case == 0 || case == 2
,
{fixedCo, varCo},
{varCo, fixedCo}
]
,
{i, len}
]
]
,
CompilationTarget -> "C"
];


The function below maps the compiled function over the layers

findCoords[pPIB_, layers_] :=

findCoords[layers_] :=
cfu, {findPrimePosInBoundarys[layers], Range[layers - 1]}];


Here is a function to display the result

showSpiral[coords_] :=
Graphics[{Point[{0, 0}], Point /@ coords}]


Now we can do

layers = 101;
(coords = findCoords[layers]) // Timing // First
showSpiral[coords]


Giving

 0.051310


For a larger value of layers, we have

layers = 1001;
(coords = findCoords[layers]) // Timing // First

1.015547

• @Öskå thanks :) Commented Jun 9, 2014 at 18:46
• That's 10 times longer than my unanswered code but 10 times faster.. :)
– Öskå
Commented Jun 9, 2014 at 18:46
• @Öskå so, what is your favorite length to efficiency ratio? ;) Commented Jun 9, 2014 at 18:47
• @JacobAkkerboom Closer to 0 than to infinity. :)
– Kuba
Commented Jun 9, 2014 at 18:49

I cannot take credit for this, but here is a rather interesting and compact implementation that I found somewhere:

F[n_] := {Re[#], Im[#]} & /@ Fold[Join[#1, Last[#1] + I^#2 Range[#2/2]] &, {0}, Range[n]]
G[n_] := Table[#[[Prime[k]]], {k, 1, PrimePi[n^2/4 + 1]}] &[F[n]]
ListPlot[G[500], AspectRatio -> Automatic, Axes -> None, Frame -> True]


Whether that is faster than the other implementations given here, I do not know; but I do think it could be tweaked to run faster.

• Wow... It's really a perfect solution! Did you find it from here ? Commented Jun 25, 2014 at 12:07

Here's my version using Join instead of Insert to build the spiral of numbers. By my measurements it's the fastest method for the smaller spiral of 200 layers, taking just 0.025 seconds to complete compared to 0.046 with Kuba's code and 0.055 with Jacob's code. But Jacob's solution wins hands down on the larger spiral with 1001 layers. That test takes just 1.14 seconds for Jacob's code while it takes 56 seconds for mine and 115 seconds for Kuba's.

addSpiralLayer[] := Module[{
matrix = {{5, 4, 3}, {6, 1, 2}, {7, 8, 9}}, layer = 2, start = 9,
perimeter, shortSide, longSide
},
Function[
perimeter = 8 layer;
shortSide = (perimeter - 4)/4;
longSide = shortSide + 2;

layer++;

matrix = Join[
{Reverse[start + shortSide + Range[longSide]]},
Transpose@Join[
{start + shortSide + longSide + Range[shortSide]},
Transpose@matrix,
{Reverse[start + Range[shortSide]]}
],
{start + 2 shortSide + longSide + Range[longSide]}
];
start = start + 2 shortSide + 2 longSide;
matrix
]
]


Here's a test with as many layers as in Kuba's example (the visualization is with one layer more).

spiral = addSpiralLayer[];
Do[spiral[];, {101}] // Timing
spiral[] // PrimeQ // Boole // Image // ColorNegate


(As an aside for future visitors: //Image//ColorNegate is faster than ArrayPlot, if you want to optimize every aspect of the problem.)

One point of confusion for some readers in regards to the code above could be that it is implemented as a closure which, as they say in the Mathematica Coobook, is "a bit outside garden-variety Mathematica". I'm referencing that book because, for those interested, its author has implemented a closure in an answer of his on this site. That example is simpler so it could be a good place to start learning how closures work.

One improvement of my code that I know but did not implement is that Range can list numbers in reversed order by having a negative step size. This could be used instead of Reverse.

Here is a method that treats the Ulam spiral as a permutation of the numbers Range[n^2]. Start with number one at the centre of the list, then go 1 right, 1 up, 2 left, 2 down, 3 right, 3 up, or

Accumulate[{(n^2 + 1)/2, +1, -n, -1, -1, +n, +n, +1, +1, +1, -n, -n, -n, ...}]


until you have a list of length n^2. This can be used in Permute to shuffle a Range[n^2] around.

ulamSpiral[n_] :=
Permute[Range[n^2],
Accumulate@
Take[Join[{n^2 + 1}/2,
Flatten@Table[(-1)^j i, {j, n}, {i, {-1, n}}, {j}]], n^2]]~
Partition~n


The final Partition[#, n] turns the flat list into the usual matrix form, and can of course be left out.

In[]:= ulamSpiral[5] // MatrixForm
Out[]//MatrixForm= (
17  16  15  14  13
18  5   4   3   12
19  6   1   2   11
20  7   8   9   10
21  22  23  24  25
)


As for the performance, it is not optimised but about twice as fast as the implementation of top on my machine.

(This is from a code golf challenge.)

• I often wondered whether there is a better way to generate lists with repeating patterns, such as {a, b, a, a, b, b, a, a, a, b, b, b, ...}, which would make this approach much clearer to the reader.
– user34757
Commented Oct 12, 2015 at 6:53
• Good question. Perhaps, you should post it as one. Commented Oct 12, 2015 at 13:06
• Here's an alternative: ulamSpiral[n_] := Partition[InversePermutation[FoldList[Plus, (n^2 + 1)/2, Take[Flatten[Table[(-1)^j i, {j, n}, {i, {-1, n}}, {j}]], n^2 - 1]]], n] Commented Mar 23, 2018 at 2:16

Yet an other approach using Graph:

Basically I'm taking advantage of PathGraph:

PathGraph@Range@25


sp[n_?OddQ] := Module[{rmat, newmat},
rmat = Reverse /@ (({1, 1} + #) & /@
Rationalize@
AbsoluteOptions[PathGraph[Range[n^2]], VertexCoordinates][[2]]);
newmat = SparseArray[Rule[#, #2] & @@@ Thread[{rmat, Range@(n^2)}]];
Permute[newmat, Cycles[Thread@{Range[1, n - Ceiling[n/2]], Range[n, Ceiling[n/2] + 1, -1]}]]]
sp@n; // AbsoluteTiming


It takes more time if n < 501 compared to kuba's code and it's faster if n > 501 (~10 sec vs ~30 sec for n = 1001). But it's obviously much slower than Jakob's solution :)

project euler 28

f[x_] := 1/4 (-1)^x (1 + 3 (-1)^x + 2 (-1)^x x)

Which[type == 1, Append[x\[Transpose], pad]\[Transpose], type == 2,