Thick End of the Wedg
Thick End of the Wedg
I came across Ulam spirals and Sacks spirals on this numberphile video. The Ulam spiral (or prime spiral or ulam cloth) is, according to Wikipedia,
a method of visualizing the prime numbers that reveals the apparent tendency of certain quadratic polynomials to generate unusually large numbers of primes.
The story goes that
it was discovered by the mathematician Stanislaw Ulam in 1963, while he was doodling during the presentation of a "long and very boring paper" at a scientific meeting.
A reveal of the Ulam spiral is some pattern to the prime number distribution relating to some "rich" (in primes) quadratics. Maybe it shouldn't be too surprising that despite the prime number distribution being random (without thinking too much about what that means) there are some patterns. After all some pattern examples of the prime distribution are quite straighforward e.g. primes after 2 must be odd, or primes must be of the form 6n+1 or 6n-1. But the Ulam spiral was a surprise when discovered and according to the numberphile video is not very well understood.
This page shows a way to create visualisations of the Ulam spirals and the Sacks spirals (an Archimedean spiral variant of a Ulam spiral) using Julia.
We'll use the Primes Package to test whether a number is prime, the Images package and Colors package for plotting the Ulam spirals, and (further down) the Gadfly package for plotting the Sacks spirals.
using Primes, Images, Colors
The spiral(n) function returns a spiral matrix of size $n \times n$ (spiralling out from the center). There are options to choose the starting number (defaults to 1) and increments (defaults to 1).
function spiral(n::Int; start::Int = 1, incr::Int = 1)
# Create vector and empty matrix
finish = start + incr * (n^2 - 1)
nums = collect(start:incr:finish)
S = zeros(Int, n, n)
# Starting location
row = ceil(Int, ((n + 1) / 2))
col = floor(Int, ((n + 1) / 2))
# Wrap numbers around matrix
S[row, col] = splice!(nums, 1)
for i = 1:n-col
k = i-1
length = 2i
S[(row - i):(row + k), (col + i)] = reverse(splice!(nums, 1:length))
if size(nums,1) == n-1
S[(row - i), (col - k):(col + k)] = reverse(nums)
break
end
S[(row - i), (col - i):(col + k)] = reverse(splice!(nums, 1:length))
S[(row - k):(row + i), (col - i)] = splice!(nums, 1:length)
S[(row + i), (col - k):(col + i)] = splice!(nums, 1:length)
end
# Return spiral matrix
S
end;
A 10 x 10 spiral matrix looks like this:
spiral(10)
The ulam(n) function first creates a spiral matrix and then plots the Ulam spiral.
function ulam(n::Int; col_prime = colorant"navy", col_background = colorant"white", start = 1, incr = 1)
S = spiral(n, start = start, incr = incr)
canvas = [isprime(S[i,j]) ? col_prime : col_background for i = 1:size(S, 1), j = 1:size(S, 2)]
display("image/png", imresize(canvas, (500, 500)))
end;
ulam(200)
You can start at any number and still find that primes appear along the diagonals. The image below is for a Ulam sprial starting at 41. Notice the very long low-left to top-right diagonal. These points all lie on $x^2-x+41$ which produces 41 primes in a row for $x=0,1,2,\dots,40$.
ulam(200, start = 41, col_prime = colorant"yellow", col_background = colorant"black")
Playing around a bit, let's see what happens when starting at 3 and incrementing by 2 i.e. only odd numbers. I can't see any obviously interesting pattern emerging.
ulam(200, start=3, incr=2, col_prime = colorant"maroon", col_background = colorant"whitesmoke")
Another thing we can do is highlight some of the quadratic polynomials that are quite rich in primes e.g. Eulers prime generating polynomial $x^2 - x + 41$ or the highly rich $4x^2 - 2x + 41$.
function ulam(n::Int, f::Function;
col_prime = colorant"red2", col_func = colorant"navy", col_background = colorant"white", start = 1, incr = 1)
S = spiral(n)
p = [f(x) for x in -n:n]
canvas = [isprime(S[i,j]) ? any(p .== S[i,j]) ? col_func : col_prime : col_background for i = 1:size(S,1), j = 1:size(S,2)]
display("image/png", imresize(canvas, (500, 500)))
end;
ulam(200, x -> 4x^2 - 2x + 41)
The Sacks spiral is a variant of the Ulam spiral where, instead of plotting numbers on a rectangular spiral, they are plotted on an Archimedean spiral.
An Archimedean spiral is the locus of points leading away from a central point with fixed speed and fixed angular velocity i.e. as you turn the angle $(\theta)$ round and round, the radius $(r)$ (or distance from the origin) grows at a constant rate. The relationship between radius and angle is
$$r(\theta) = a + b \, \theta$$with $a$ & $b$ $\in \mathbb{R}$. Parameter $a$ turns the spiral and parameter $b$ controls the distance between the spiral arms.
To create a Sacks spiral we need to first place all the whole numbers $(i = 0, \, 1, \, 2, \, 3, \, \dots)$ along an Archimedean spiral such that all square numbers $(i = 0, \, 1, \, 4, \, 9, \, \dots)$ are on the positive x axis. Then, with all the numbers in place, we plot only the prime number positions.
For plotting it will be easier to use a plotting package instead of an image package as now we are plotting point locations as opposed to colouring rectangular canvas pixels. The Gadfly package is easy to use and produces beautiful quality plots.
using Gadfly
Before plotting we need to work out the Sacks version of the Archimedean spiral and tackle the polar coordinates system. Hopefully this is a quick digression.
If the square numbers are all placed on the positive x axis then
$$\theta(i^2)=2\pi i \quad \text{or} \quad \theta(i)=2\pi \sqrt{i}$$
i.e. $i^2$ is placed (on the x axis) after $i$ rotations. Substituting into the equation for $r$ gives
.
To calculate the value of $a$ note that $$r(0) = 0 \implies a + b \, 2\pi\sqrt{0} = 0 \implies a=0$$.
Since $b$ just controls the distance between spiral arms (effectively the "zoom") we can set it to any value and it's convenient to let $b=\frac{1}{2\pi}$. We then have:
$$r(i) = \sqrt{i} \quad \text{and} \quad \theta(i) = 2\pi \, r(i)$$For plotting, we can calculate the Cartesian coordinates from the polar coordinates using:
$$x = r \cos(\theta) \quad \text{and} \quad y = r \sin(\theta)$$The chart below illustrates this method of placing the whole numbers on the Archemedean spiral. The function archimedean_pts($n$) creates the location of the whole numbers - will plot the points. And archimedean_line($n$) creates the location of a range of real numbers - will plot the path. If you count the dots (starting with $0$) from the origin you'll see that the square numbers $(0, \, 1, \, 4, \, 9, \, \dots)$ are all on the positive x axis.
function archimedean_pts(n)
r(i) = sqrt(i)
θ(i) = 2π * r(i)
range = 0:n
x = [r(i) * cos(θ(i)) for i in range]
y = [r(i) * sin(θ(i)) for i in range]
return x, y
end;
function archimedean_line(n)
r(i) = sqrt(i)
θ(i) = 2π * r(i)
range = linspace(0, n, 10000)
x = [r(i) * cos(θ(i)) for i in range]
y = [r(i) * sin(θ(i)) for i in range]
return x, y
end;
x1, y1 = archimedean_pts(100)
x2, y2 = archimedean_line(100)
set_default_plot_size(12cm, 12cm)
plot(layer(x = x1, y = y1, Geom.point, Theme(default_color = colorant"blue4")),
layer(x = x2, y = y2, Geom.path, Theme(default_color = colorant"grey")))
Now we're ready to create the Sacks spiral.
The function sacks(n) is similar to archimedean(n) but only produces the prime number points. And the function plot_sacks(n; ...) plots the spiral.
function sacks(n)
x = Array{Float64}(0)
y = Array{Float64}(0)
for i in 1:n
if isprime(i)
r = sqrt(i)
θ = 2π * r
push!(x, r * cos(θ))
push!(y, r * sin(θ))
end
end
return x, y
end;
function plot_sacks(n; point_size = 0.05cm, highlight_width = 0.005cm, default_color = colorant"blue4")
xs, ys = sacks(n)
theme = Theme(point_size = point_size, highlight_width = highlight_width, default_color = default_color)
guide = Guide.xticks(ticks = nothing), Guide.yticks(ticks = nothing), Guide.xlabel(nothing), Guide.ylabel(nothing)
return plot(x = xs, y = ys, guide..., theme)
end;
Numbers up to 10,000:
set_default_plot_size(15cm, 15cm)
plot_sacks(10000)
Numbers up to 40,000:
set_default_plot_size(25cm, 25cm)
plot_sacks(40000, point_size = 0.04cm, highlight_width = 0.005cm, default_color = colorant"darkred")