Turn Left at the Next Prime

A while ago, John Cook made a post about Gaussian integer walks. See his post for the details, but at the end, he provides some Python code which he uses to generate the walks/plots. However, just for the simulation starting at $127+130i$, it took several minutes for the whole process, including plotting, which just seemed completely unreasonable to me as a Julia user. So of course, I wanted to see if I could reimplement that code in Julia to be much faster and yet not sacrifice ease of readability.

The basic premise is this: start at some point in the complex plane (with integer coefficients), facing east. Move unit-by-unit until the number is a Gaussian prime; if Gaussian prime, turn left.

The first thing we need is a function to determine if a complex number is a Gaussian prime (Julia uses im for the imaginary unit):

function isgaussprime(z::Complex{<:Integer})
    a, b = real(z), imag(z)
    if a * b != 0
        return isprime(a^2 + b^2)
    else
        c = abs(a+b)
        return isprime(c) && c % 4 == 3
    end
end

To do the Gaussian walk, we are going to write the function walk that starts at some point and keeps walking up to some limit (these walks all end up spirally, so we need some limit):

function walk(start, limit = 10)
    points = [start]
    z = start
    delta = 1
    
    while limit > 0
        z = z + delta
        push!(points, z)
        
        isgaussprime(z) && (delta *= im)
        z == start && break
        limit -= 1
    end
    
    points
end

To support the interactivity, I also have some other related functions, see the attached notebook.

(To be fair, the code above doesn't quite do the same thing as the original code in a couple of ways, but it should still be straightforward to understand. There is also a huge speed improvement, to the point that this can be done interactively.)

The cool part is what happens when we plot the walk for various complex numbers. To verify the implementation against the original code, we'll use the same examples. First, for $3+5i$,

The one that took a few minutes to plot:

The behavior honestly is very weird. A larger Gaussian prime doesn't necessarily mean a longer period. Although a lot of these are symmetrical in some way, a lot aren't; for example, $127+144i$ looks very not symmetrical (and is pretty short):

The full Pluto notebook is here for your perusal.

To determine the period of the cycle for a walk starting from some point:

function walklimit(start)
    limit = 1
    z = start
    delta = 1
    
    while true
        z = z + delta
        
        isgaussprime(z) && (delta *= im)
        z == start && break
        limit += 1
    end
    
    limit + 1
end

To determine the bounds of the cycle in order to draw a graph with the best range:

function walkbounds(start)
    min_real, max_real, min_imag, max_imag = real(start), real(start), imag(start), imag(start)
    z = start
    delta = 1
    
    while true
        z = z + delta
        
        min_real = min(min_real, real(z))
        max_real = max(max_real, real(z))
        min_imag = min(min_imag, imag(z))
        max_imag = max(max_imag, imag(z))
        
        isgaussprime(z) && (delta *= im)
        z == start && break
    end
    
    (min_real, max_real), (min_imag, max_imag)
end