Here's some Python code that plots this prime spiral path for primes less than num, with the grid cell colour determined by the number of times the path has visited it. Unvisited cells are black, visited cells are coloured by cycling through {red, yellow, green, cyan, blue, magenta}. Thus cells that have been visited once are red, cells that have been visited twice are yellow, thrice-visited cells are green, etc. When we get to the end of the colour list we go back to red, so a cell that's been visited 7, 13, ... times is red.
The origin cell (i.e., the cell containing 1) is re-coloured white at the end of the plotting process, but it's not easy to find when num is large.
This code uses the 3rd-party modules Numpy and Pillow to do the plotting. It was written for Python 3, but it will still run correctly on Python 2.
#!/usr/bin/env python3
''' New prime spiral by Karl
Saves path as a PNG
See http://math.stackexchange.com/q/2072308/207316
Written by PM 2Ring 2017.1.1
'''
from collections import defaultdict
from itertools import cycle, islice
from operator import itemgetter
import numpy as np
from PIL import Image
def primes_bool(n):
''' Return a boolean list of all primes 0 <= p < n '''
s = [True] * n
s[:2] = [False] * 2
for i in range(2, int(n**0.5) + 1):
if s[i]:
s[i*i : n : i] = [False] * (1 + (n - 1)//i - i)
return s
def vadd(u, v):
return u[0] + v[0], u[1] + v[1]
step = cycle([(1, 0), (0, 1), (-1, 0), (0, -1)])
st = next(step)
current = (0, 0)
points = defaultdict(int)
points[current] = 1
# perform steps
num = 1000000
for v in islice(primes_bool(num), 2, None):
current = vadd(current, st)
points[current] += 1
if v:
st = next(step)
# find bounds, adding a border of 1
keys = points.keys()
ig0, ig1 = itemgetter(0), itemgetter(1)
lox = min(keys, key=ig0)[0] - 1
hix = max(keys, key=ig0)[0] + 1
loy = min(keys, key=ig1)[1] - 1
hiy = max(keys, key=ig1)[1] + 1
width = 1 + hix - lox
height = 1 + hiy - loy
print('x:', lox, hix, width)
print('y:', loy, hiy, height)
print('max crossings:', max(points.values()))
# convert to image grid
#KRYGCBMW
palette = np.array([
[0, 0, 0],
[255, 0, 0],
[255, 255, 0],
[0, 255, 0],
[0, 255, 255],
[0, 0, 255],
[255, 0, 255],
[255, 255, 255],
], dtype='uint8')
grid = np.zeros((height, width, 3), dtype='uint8')
for (x, y), v in points.items():
grid[y - loy, x - lox] = palette[(v - 1) % 6 + 1]
del points
# reset location of 1, i.e.,(0, 0) to white
grid[-loy, -lox] = palette[7]
img = Image.fromarray(grid)
#img.show()
img.save('primespiral.png')
Uncomment the #img.show() line to display the image at full size with an external viewer. This normally uses display from the ImageMagick package on *nix systems; on Windows, it saves the image to a temporary BMP file, and allegedly uses the standard BMP display utility to show it.
Here's the output for num = 1000000:
And here's the output for num = 10000000, scaled by a factor of 1/16:
If you'd like to see when the current record number of visits is reached, replace the main path generating loop with this code:
# perform steps
num = 10000000
maxvisits = 1
for i, v in enumerate(islice(primes_bool(num), 2, None), 2):
current = vadd(current, st)
points[current] += 1
if points[current] > maxvisits:
maxvisits = points[current]
print('{} visits at {}'.format(maxvisits, i))
if v:
st = next(step)
Here's what it prints for num = 10000000
2 visits at 15
3 visits at 35
4 visits at 47
5 visits at 81
6 visits at 7087
7 visits at 7399
8 visits at 19865
9 visits at 19913
10 visits at 24087
11 visits at 24279
12 visits at 408257
13 visits at 409303
14 visits at 2042205
15 visits at 5262017
16 visits at 5262089
17 visits at 6189393
18 visits at 6435851
19 visits at 6435983
Here's a version that uses turtle graphics, so you can watch the path being drawn. The turtle module is included in most modern Python installations.
#!/usr/bin/env python3
''' New prime spiral by Karl
See http://math.stackexchange.com/q/2072308/207316
https://oeis.org/A268463
Written by PM 2Ring 2017.1.1
turtle version 2017.1.3
'''
import sys
import turtle
from collections import defaultdict
from itertools import islice
from operator import itemgetter
def primes_bool(n):
''' Return a boolean list of all primes 0 <= p < n '''
s = [True] * n
s[:2] = [False] * 2
for i in range(2, int(n**0.5) + 1):
if s[i]:
s[i*i : n : i] = [False] * (1 + (n - 1)//i - i)
return s
def vadd(u, v):
return u[0] + v[0], u[1] + v[1]
steps = [(1, 0), (0, 1), (-1, 0), (0, -1)]
step_index = 0
st = steps[step_index]
current = (0, 0)
points = defaultdict(int)
points[current] = 1
# turtle stuff
palette = [
'black',
'red',
'yellow',
'green',
'cyan',
'blue',
'magenta',
'white',
]
win_width, win_height = 3000, 3000
pensize = 3
turtle.setup()
turtle.screensize(win_width, win_height, palette[0])
t = turtle.Turtle()
t.hideturtle()
t.speed(0)
t.pensize(pensize)
t.color(palette[-1])
t.forward(pensize)
# perform steps
num = 100000
maxvisits = 1
for i, v in enumerate(islice(primes_bool(num), 2, None), 2):
print('\r {} '.format(i), end='', file=sys.stderr)
current = vadd(current, st)
points[current] += 1
pc = points[current]
if pc > maxvisits:
maxvisits = pc
print('{} visits at {}'.format(maxvisits, i))
if v:
step_index = (step_index + 1) % 4
st = steps[step_index]
t.right(90)
t.color(palette[(pc - 1) % 6 + 1])
t.forward(pensize)
turtle.done()
print()
To change the scale of the drawing just change the pensize variable to another positive integer. Note that the pen is rounded, so the path will have noticeable rounded corners for pensize > 3; this may or may not be desirable. :)
You can easily modify this script to draw some alternate prime paths. All primes > 3 are of the form $6n\pm 1$. So we can turn right if prime $i \mod 6 = 1$, or turn left for the other primes.
To do that change the code in the if v: block to
if v:
dx = 1 if i % 6 == 1 else -1
step_index = (step_index + dx) % 4
st = steps[step_index]
t.right(90 * dx)
All primes > 2 are of the form $4n\pm 1$. This is perhaps more important than the $6n\pm 1$ distinction since primes of the form $4n-1$ are also prime in the Gaussian integers; 2 and primes of the form $4n+1$ can be written as the sum of two squares and can thus be factorised in the Gaussian integers.
So we can make the path turn right for $4n+1$ primes and left for the other primes by changing the modulus in the dx assignment:
dx = 1 if i % 4 == 1 else -1
Here's a live version that runs on the SageMathCell server. If you make num too high the server will abort the job, but you should be able to go up to 10000000 or so.
