P3.2.2The Gaussian Prime spiral
A Gaussian integer is a complex number whose real and imaginary parts are both integers. A Gaussian prime is a Gaussian integer \(x+iy\) such that either:
- one of \(x\) and \(y\) is zero and the other is a prime number of the form \(4n+3\) or \(-(4n+3)\) for some integer \(n \ge 0\); or
- both \(x\) and \(y\) are non-zero and \(x^2+y^2\) is prime.
Consider the sequence of Gaussian integers traced out by an imaginary particle, initially at \(c_0\), moving in the complex plane according to the following rule: it takes integer steps in its current direction (\(\pm 1\) in either the real or imaginary direction), but turns left if it encounters a Gaussian prime. Its initial direction is in the positive real direction (\(\Delta c = 1 + 0i\Rightarrow \Delta x = 1, \Delta y = 0\)). The path traced out by the particle is called a Gaussian prime spiral.
Write a program to plot the Gaussian prime spiral starting at \(c_0 = 5 + 23i\).
Show solution
The following code plots the Gaussian prime spiral starting at \(c_0 = 5 + 23i\).
import numpy as np
import matplotlib.pyplot as plt
def is_prime(a):
"""Return True if a is a prime number."""
if a == 2:
return True
if not a % 2:
return False
for d in range(3, int(np.sqrt(a)) + 1, 2):
if not a % d:
return False
return True
def is_gaussian_prime(c):
"""Return True if c = x + iy is a Gaussian prime (otherwise False)."""
x, y = c
if x == 0 or y == 0:
a = abs(x) or abs(y)
return is_prime(a) and a % 4 == 3
a = x**2 + y**2
return is_prime(a)
def turn_left(c):
"""Permute (1,0) -> (0,1) -> (-1,0) -> (0,-1) -> (1,0)."""
dx, dy = c
return -dy, dx
# Starting point and direction
x, y = x0, y0 = 5, 23
dx, dy = 1, 0
# Keep track of the points on our spiral path in these lists
pathx, pathy = [x], [y]
# It isn't known whether every iteration produces a closed loop so
# set a maximum number of steps to take
max_steps = 10000
step = 0
while step < max_steps:
step += 1
x, y = x + dx, y + dy
pathx.append(x)
pathy.append(y)
if (x, y) == (x0, y0):
print(f"Returned to ({x}, {y}) in {step} steps.")
break
if is_gaussian_prime((x, y)):
dx, dy = turn_left((dx, dy))
else:
print(f"max_steps = {max_steps} reached.")
plt.plot(pathx, pathy, lw=2)
plt.savefig("p3/p3_2_2_1__gaussianspiral.png") # noweb
plt.show()