# Ulam Spiral Variation

Its always fun to play with prime numbers.

Why not. In math they mean a lot, in computers they are important, and in recreational mathematics they are… ahem… interesting.

## Ulam Spiral

In 1963 mathematician Stanislaw Ulam published this spiral in Martin Gardner’s Mathematical games in scientific american:

From that small Scientific american cover we can see some interesting trends DIAGONALS. Some interesting diagonals appear. What is this? there is some bizarre order in prime numbers? You can read more in wikipedia.

Now lets (for the fun of it) try to see what happens with other forms of turning the number line into s surface.

## Diagonalization

Lets try to simply ‘diagonalize’ the number-line. Here ‘1’ starts in the top left corner then ‘2’ is placed in its left and subsequent numbers follow a diagonal:

Now lets mark the prime numbers:

Lets juice it up to visualize the prime numbers up to one million (999983). What? there are diagonals? however here the diagonals follow the numbers. Meaning these diagonals are only some prime dense regions. However there are some subtle vertical and horizontal lines. As I said prime numbers are always interesting.

## Hilbert

Hilbert curve is a space filling fractal (there are higher dimensions variants). I have played with them for bioinformatics purposes.

Now lets try to put some primes in here.

What? Diagonals? again.

Well, I am not sure what this means, but I will never stop playing with this stuff.

## Code Time

Here is the code needed for the Ulam spiral. The important part is the last few lines

``````# -*- coding: utf-8 -*-
"""
Created on Sun Dec 30 16:28:25 2018

@author: lahat
"""
import matplotlib.pyplot as plt
import sympy
import numpy as np

"""
This is a module to convert between one dimensional distance along a
`Hilbert curve`_, :math:`h`, and N-dimensional coordinates,
:math:`(x_0, x_1, ... x_N)`.  The two important parameters are :math:`N`
(the number of dimensions, must be > 0) and :math:`p` (the number of
iterations used in constructing the Hilbert curve, must be > 0).

We consider an N-dimensional `hypercube`_ of side length :math:`2^p`.
This hypercube contains :math:`2^{N p}` unit hypercubes (:math:`2^p` along
each dimension).  The number of unit hypercubes determine the possible
discrete distances along the Hilbert curve (indexed from :math:`0` to
:math:`2^{N p} - 1`).  The image below illustrates the situation for
:math:`N=2` and :math:`p=3`.

.. figure:: nD=2_p=3.png

This is the third iteration (:math:`p=3`) of the Hilbert curve in two
(:math:`N=2`) dimensions.  Distances, :math:`h`, along the curve are
labeled from 0 to 63 (i.e. from 0 to :math:`2^{N p}-1`).  The provided
functions translate between N-dimensional coordinates and the one
dimensional distance.  For example, between (:math:`x_0=4, x_1=6`) and
:math:`h=36`.

Reference
=========

This module is based on the C code provided in the 2004 article
"Programming the Hilbert Curve" by John Skilling,

I was also helped by the discussion in the following stackoverflow post,

- `mapping-n-dimensional-value-to-a-point-on-hilbert-curve`_

which points out a typo in the source code of the paper.  The Skilling code
provides two functions ``TransposetoAxes`` and ``AxestoTranspose``.  In this
case, Transpose refers to a specific packing of the integer that represents
distance along the Hilbert curve (see below for details) and
Axes refer to the N-dimensional coordinates.  Below is an excerpt of the docs
from that code that appears in the paper by Skilling, ::

//+++++++++++++++++++++++++++ PUBLIC-DOMAIN SOFTWARE ++++++++++++++++++++++++++
// Functions: TransposetoAxes  AxestoTranspose
// Purpose:   Transform in-place between Hilbert transpose and geometrical axes
// Example:   b=5 bits for each of n=3 coordinates.
//            15-bit Hilbert integer = A B C D E F G H I J K L M N O is stored
//            as its Transpose
//                   X = A D G J M                X|
//                   X = B E H K N    <------->       | /X
//                   X = C F I L O               axes |/
//                          high  low                    0------ X
//            Axes are stored conveniently as b-bit integers.
// Author:    John Skilling  20 Apr 2001 to 11 Oct 2003

.. _Hilbert curve: https://en.wikipedia.org/wiki/Hilbert_curve
.. _hypercube: https://en.wikipedia.org/wiki/Hypercube
.. _mapping-n-dimensional-value-to-a-point-on-hilbert-curve: http://stackoverflow.com/questions/499166/mapping-n-dimensional-value-to-a-point-on-hilbert-curve/10384110#10384110
"""

def _binary_repr(num, width):
"""Return a binary string representation of `num` zero padded to `width`
bits."""
return format(num, 'b').zfill(width)

def _hilbert_integer_to_transpose(h, p, N):
"""Store a hilbert integer (`h`) as its transpose (`x`).

:param h: integer distance along hilbert curve
:type h: ``int``
:param p: number of iterations in Hilbert curve
:type p: ``int``
:param N: number of dimensions
:type N: ``int``
"""
h_bit_str = _binary_repr(h, p*N)
x = [int(h_bit_str[i::N], 2) for i in range(N)]
return x

def _transpose_to_hilbert_integer(x, p, N):
"""Restore a hilbert integer (`h`) from its transpose (`x`).

:param x: the transpose of a hilbert integer (N components of length p)
:type x: ``list`` of ``int``
:param p: number of iterations in hilbert curve
:type p: ``int``
:param N: number of dimensions
:type N: ``int``
"""
x_bit_str = [_binary_repr(x[i], p) for i in range(N)]
h = int(''.join([y[i] for i in range(p) for y in x_bit_str]), 2)
return h

def coordinates_from_distance(h, p, N):
"""Return the coordinates for a given hilbert distance.

:param h: integer distance along the curve
:type h: ``int``
:param p: side length of hypercube is 2^p
:type p: ``int``
:param N: number of dimensions
:type N: ``int``
"""
x = _hilbert_integer_to_transpose(h, p, N)
Z = 2 << (p-1)

# Gray decode by H ^ (H/2)
t = x[N-1] >> 1
for i in range(N-1, 0, -1):
x[i] ^= x[i-1]
x ^= t

# Undo excess work
Q = 2
while Q != Z:
P = Q - 1
for i in range(N-1, -1, -1):
if x[i] & Q:
# invert
x ^= P
else:
# excchange
t = (x ^ x[i]) & P
x ^= t
x[i] ^= t
Q <<= 1

# done
return x

def distance_from_coordinates(x, p, N):
"""Return the hilbert distance for a given set of coordinates.

:param x: coordinates len(x) = N
:type x: ``list`` of ``int``
:param p: side length of hypercube is 2^p
:type p: ``int``
:param N: number of dimensions
:type N: ``int``
"""
M = 1 << (p - 1)

# Inverse undo excess work
Q = M
while Q > 1:
P = Q - 1
for i in range(N):
if x[i] & Q:
x ^= P
else:
t = (x ^ x[i]) & P
x ^= t
x[i] ^= t
Q >>= 1

# Gray encode
for i in range(1, N):
x[i] ^= x[i-1]
t = 0
Q = M
while Q > 1:
if x[N-1] & Q:
t ^= Q - 1
Q >>= 1
for i in range(N):
x[i] ^= t

h = _transpose_to_hilbert_integer(x, p, N)
return h

def hilbert(l):
w = int((round(len(l)**0.5+0.5))*1.4)

to_matrix={}
for n,i in enumerate(l):
x,y = coordinates_from_distance(n,w,2)
to_matrix[x,y]=i

matrix = np.zeros([1+max([i for i  in to_matrix.keys()]),
1+max([i for i  in to_matrix.keys()])])
for a,b in to_matrix.items():
matrix[a]=b
return matrix

def HilberColor(*l):
if len(l)==1:
a=l
b=a
c=a
elif len(l)==2:
a=l
b=l
c=b
elif len(l)==3:
a=l
b=l
c=l
elif len(l)>3: return None
a,b,c = [list(l) for l in [a,b,c]]
M = (max([len(i) for i in [a,b,c]]))
a += *(M-len(a))
b += *(M-len(b))
c += *(M-len(c))
a,b,c = [hilbert(i) for i in [a,b,c]]
a,b,c = [i-np.nanmin(i) for i in [a,b,c]]
a,b,c = [i/(np.nanmax(i)+0.0001) for i in [a,b,c]]

return np.dstack([a,b,c])

def all_primes(start, end):
return set(sympy.sieve.primerange(start, end))

#Important stuff

N=7 # how many iterations of Ulaminess? we want
primes = all_primes(0,4**N) # makes list of all primes needed
plt.figure(figsize = (15,15)) # figure size
plt.imshow(HilberColor([i not in primes for i in range(1,4**N)]))

# this line needs some unpaking.
# what matters is:
# [i not in primes for i in range(1,4**N)]
# This checks if 'i' is in the 'primes' list and adds the result
# (true or false) to a list that then gets turned into a Hilbert curve
# ('HilbertColor') then the Hilbert curve is plotted ('plt.imshow')

plt.axis('off') # removes the axis numbers

``````

Now the code for the ‘diagonalization’ of primes.

``````
import matplotlib.pyplot as plt
import sympy
import numpy as np
from random import random
import math
#%%
def all_primes(start, end):
return set(sympy.sieve.primerange(start, end))
#%%
x,y=0,0
n=1
N=1000000
D={}
primes=all_primes(1,N)

for i in range(1,N):
if i%1000000==0:print('checking coordinates: {}'.format(i/N))
D[x,y]=i
if x==0:
x=n
n+=1
y=0
else:
x-=1
y+=1

i=0
matrix=np.zeros([n,n])
for a,b in D.items():
if i%1000000==0:print('checking primes: {}'.format(i/N))
i+=1
if b in primes:
matrix[a]=1
#%%
plt.figure(figsize = (25,25))
plt.imshow(matrix,
cmap='Greys')
plt.axis('off')
plt.xlim([0,n/2])
plt.ylim([0,n/2])``````