The previous post mentioned a Math Overflow question about unexpected mathematical images, and reproduced one that looks like field grass. This post reproduces another set of images from that post.

Start anywhere in the complex plane with integer coordinates and walk west one unit at a time until you run into Gaussian prime [1]. Then turn left (counterclockwise) 90° and keep taking unit steps. Apparently this process will often (always?) return you to your starting point.

Different starting points lead to different patterns. Here’s an example given in the post, starting at 3 + 5i.

Here’s a more complex walk starting at 27 + 30i.

I tried starting at 127 + 131i and got a simple, uninteresting image. I tried again starting at 127 + 130i and got something much more complicated. I didn’t time it, but it took several minutes to plot.

Here’s the code that made the plots. (Note that Python uses j rather than i for imaginary unit.)

from sympy import isprime
import matplotlib.pyplot as plt

def isgaussprime(z: complex):
    a, b = int(z.real), int(z.imag)
    if a*b != 0:
        return isprime(a**2 + b**2)
    else:
        c = abs(a+b)
        return isprime(c) and c % 4 == 3

def connect(z1: complex, z2: complex):
    plt.plot([z1.real, z2.real], [z1.imag, z2.imag], 'b')
    
start = 127 + 130j
#start = 3 + 5j
step = 1
z = start
next = None

while next != start:
    next = z + step
    connect(z, next)
    if isgaussprime(next):
        step *= 1j
    z = next

plt.axes().set_aspect(1)    
plt.show()

Related posts

• Square root of a Gaussian integer
• Random walk mod p
• A knight’s tour magic square

[1] If a and b are integers, then a + bi is called a Gaussian integer. A Gaussian integer is a Gaussian prime if (1) both a and b are non-zero and a² + b² is prime, or (2) one of a or b is zero, and the absolute value of the non-zero part is a prime congruent to 3 mod 4.

Why is this definition so complicated? It’s actually a theorem. There’s a natural generalization of what it means to be prime in a commutative ring, and it works out that an element the Gaussian integers is prime if and only if the above criteria hold.

In general, a non-zero element p of a commutative ring R is prime if whenever p divides a product ab, p must either divide a or divide b.

The function σ_k takes an integer n and returns the sum of the kth powers of divisors of n. For example, the divisors of 14 are 1, 2, 4, 7, and 14. If we set k = 3 we get

σ₃(n) = 1³ + 2³ + 4³ + 7³ + 14³ = 3096.

A couple special cases may use different notation.

• σ₀(n) is the number of divisors n and is usually denoted d(n), as in the previous post.
• σ₁(n) is the sum of the divisors of n and the function is usually written σ(n) with no subscript.

In Python you can compute σ_k(n) using divisor_sigma from SymPy. You can get a list of the divisors of n using the function divisors, so the bit of code below illustrates that divisor_sigma computes what it’s supposed to compute.

    n, k = 365, 4
    a = divisor_sigma(n, k)
    b = sum(d**k for d in divisors(n))
    assert(a == b)

The Wikipedia article on σ_k gives graphs for k = 1, 2, and 3 and these graphs imply that σ_k gets smoother as k increases. Here is a similar graph to those in the article.

The plots definitely get smoother as k increases, but the plots are not on the same vertical scale. In order to make the plots more comparable, let’s look at the kth root of σ_k(n). This amounts to taking the Lebesgue k norm of the divisors of n.

Now that the curves are on a more similar scale, let’s plot them all on a single plot rather than in three subplots.

If we leave out k = 1 and add k = 4, we get a similar plot.

The plot for k = 2 that looked smooth compared to k = 1 now looks rough compared to k = 3 and 4.

A while back I wrote about continued fractions of square roots. That post cited a theorem that if d is not a perfect square, then the continued fraction representation of d is periodic. The period consists of a palindrome followed by 2⌊√d⌋. See that post for details and examples.

One thing the post did not address is the length of the period. The post gave the example that the continued fraction for √5 has period 1, i.e. the palindrome part is empty.

There’s a theorem [1] that says this pattern happens if and only if d = n² + 1. That is, the continued fraction for √d is periodic with period 1 if and only if d is one more than a square. So if we wanted to find the continued fraction expression for √26, we know it would have period 1. And because each period ends in 2⌊√26⌋ = 10, we know all the coefficients after the initial 5 are equal to 10.

[1] Samuel S. Wagstaff, Jr. The Joy of Factoring. Theorem 6.15.

Take any positive integer d that is not a multiple of 2 or 5. Then there is some integer k such that d × k has only 1’s in its decimal representation. For example, take d = 13. We have

13 × 8457 = 111111.

Or if we take d = 27,

27 × 4115226337448559670781893 = 111111111111111111111111111.

Let’s change our perspective and start with the string of 1’s. If d is not a multiple of 2 or 5, then there is some number made up of only 1’s that is divisible by d. And in fact, the number of 1’s is no more than d.

This theorem generalizes to any integer base b > 1. If d is relatively prime to b, then there is a base b number with d or fewer “digits” which is divisible by d [1].

The following Python code allows us to find the number k such that d × k has only 1’s in its base b representation, provided k is relatively prime to b. It returns the number of 1’s we need to string together to find a multiple of k. If k shares a factor with b, the code returns 0 because no string of 1’s will ever be divisible by k.

    from math import gcd
    
    def all_ones(n, b = 10):
        return sum(b**i for i in range(n))
    
    def find_multiple(k, b = 10):
        if gcd(k, b) > 1:
            return 0
        for n in range(1, k+1):
            if all_ones(n, b) % k == 0:
                return n

The two Python functions above default to base 10 if a base isn’t provided.

We could find a multiple of 5 whose hexadecimal representation is all 1’s by calling

    print(find_multiple(5, 16))

and this tells us that 11111_sixteen is a multiple of 5, and in fact

5 × 369d_sixteen = 11111_sixteen.

[1] Elmer K. Hayashi. Factoring Integers Whose Digits Are all Ones. Mathematics Magazine, Vol. 49, No. 1 (Jan., 1976), pp. 19-22

Imagine seeing the following calculation:

The correct result is

and so the first calculation is off by 25 orders of magnitude.

But there’s a variation on the calculation above that is correct! A theorem by Édouard Lucas from 1872 that says for p prime and for any nonnegative integers m and n,

So while the initial calculation was grossly wrong as stated, it is perfectly correct mod 19. If you divide 487343696971437395556698010 by 19 you’ll get a remainder of 10.

A stronger versions of Lucas’ theorem [1] says that if p is at least 5, then you can replace mod p with mod p³. This is a stronger conclusion because it says not only is the difference between the left and right side of the congruence divisible by p, it’s also divisible by p² and p³.

In our example, not only is the remainder 10 when 487343696971437395556698010 is divided by 19, the remainder is also 10 when dividing by 19² = 361 and 19³ = 6859.

More on binomial coefficients

• Binomial coefficient trick
• Maximum gap between binomial coefficients
• Lebesgue’s proof of Weierstrass’ theorem

[1] V. Brun, J. O. Stubban, J. E. Fjeldstad, L. Tambs, K. E. Aubert, W. Ljunggren, E. Jacobsthal. On the divisibility of the difference between two binomial coefficients, Den 11te Skandinaviske Matematikerkongress, Trondheim, 1949, 42–54.

Someone asked an interesting question on MathOverflow: given an odd number, can you always flip a bit in its binary representation to make it prime?

It turns out the answer is no, but apparently it is very often the case an odd number is just a bit flip away from being prime. I find that surprising.

Someone pointed out that 2131099 is not a bit flip away from a prime, and that this may be the smallest example [1]. The counterexample 2131099 is itself prime, so you could ask whether an odd number is either a prime or a bit flip away from a prime. Is this always the case? If not, is it often the case?

The MathOverflow question was stated in terms of Hamming distance, counting the number of bits in which two bit sequences differ. It asked whether odd numbers are always Hamming distance 1 away from a prime. My restatement of the question asks whether the Hamming distance is always at most 1, or how often it is no more than 1.

You could ask more generally about the Hamming distance to the nearest prime. Is it bounded, if not by 1, then by another finite number? If so, what is the smallest such bound? What is the probability that its value is 1? Etc.

This ties into a couple of other things I’ve blogged about. A few weeks ago I wrote about new work on the problem of finding the proportion of odd numbers that can be written as the sum of a power of 2 and a prime. That’s a little different problem since bit flipping is taking the XOR (exclusive or) and not always the same as addition. It also leaves out the possibility of flipping a bit beyond the most significant bit of the number, i.e. adding to a number n a power of 2 greater than n.

Another related post is on the Rowhammer attack on public key cryptography. By flipping a bit in the product of two primes, you can produce a number which is much easier to factor.

These two posts suggest a variation on the original problem where we disallow flipping bits higher than the most significant bit of n. So giving a k-bit number n, how often can we flip one of its k bits and produce a prime?

[1] Note that the bit flipped may be higher than the most significant bit of the number, unless ruled out as in the paragraph above. Several people have asked “What about 85?” It is true that flipping any of the seven lowest bits of 85 will not yield a prime. But flipping a zero bit in a more significant position will give a prime. For example, 1024 + 85 is prime. Bur for 2131099 it is not possible to add any larger power of 2 to the number and produce a prime.

In previous post I showed how to compute the square root of a complex number. I gave as an example that computed the square root of 5 + 12i to be 3 + 2i.

(Of course complex numbers have two square roots, but for convenience I’ll speak of the output of the algorithm as the square root. The other root is just its negative.)

I chose z = x + iy in my example so that x² + y² would be a perfect square because this simplified the exposition.That is, I designed the example so that the first step would yield an integer. But I didn’t expect that the next two steps in the algorithm would also yield integers. Does that always happen or did I get lucky?

It does not always happen. My example was based on the Pythagorean triple (5, 12, 13). But (50, 120, 130) is also a Pythagorean triple, and the square root of z = 50 + 130i does not have integer real and imaginary parts.

At this point it will help to introduce a little terminology. A Gaussian integer is a complex number with integer real and imaginary parts. We could summarize the discussion above by saying the square root of 5 + 12i is a Gaussian integer but the square root of 50 + 120i is not.

While (5, 12, 13) and (50, 120, 130) are both are Pythagorean triples, only the first is a primitive Pythagorean triple, one in which the first two components are relatively prime. So maybe we should look at primitive Pythagorean triples.

There is a theorem going all the way back to Euclid that primitive Pythagorean triples have the form

(m² – n², 2mn, m² + n²)

where m and n are relatively prime integers and one of m and n is even.

Using the algorithm in the previous post, we can show that if x and y are the first components of a triple of the form above, then the square root of the Gaussian integer x + iy is also a Gaussian integer. Specifically we have

ℓ = m² + n²
u = √((m² + n² + m² – n²)/2) = m
v = sign(y) √((m² + n² – m² + n²)/2) = sign(y) n

This means u and v are integers and so the square root of x + iy is the Gaussian integer u + iv.

There is one wrinkle in the discussion above. Our statement of Euclid’s theorem left out the assumption that we’ve ordered the sides of our triangle so that the odd side comes first. Primitive Pythagorean triples could also have the form

(2mn, m² – n², m² + n²)

and in that case our proof would break down.

Surprisingly, there’s an asymmetry between the real and imaginary parts of the number whose square root we want to compute. It matters that x is odd and y is even. For example, √(5 + 12i) is a Gaussian integer but √(12 + 5i) is not!

So a more careful statement of the theorem above is that if x and y are the first two components of a primitive Pythagorean triple, the square root of x + iy is a Gaussian integer.

Suppose you fill two n×n matrices with random integers. What is the probability that the determinants of the two matrices are relatively prime? By “random integers” we mean that the integers are chosen from a finite interval, and we take the limit as the size of the interval grows to encompass all integers.

Let Δ(n) be the probability that two random integer matrices of size n have relatively prime determinants. The function Δ(n) is a strictly decreasing function of n.

The value of Δ(1) is known exactly. It is the probability that two random integers are relatively prime, which is well known to be 6/π². I’ve probably blogged about this before.

The limit of Δ(n) as n goes to infinity is known as the Hafner-Sarnak-McCurley constant [1], which has been computed to be 0.3532363719…

Since Δ(n) is a decreasing function, the limit is also a lower bound for all n.

Python simulation

Here is some Python code to experiment with the math discussed above. We’ll first do a simulation to show that we get close to 6/π² for the proportion of relatively prime pairs of integers. Then we look at random 2×2 determinants.

    from sympy import gcd
    from numpy.random import randint
    from numpy import pi
    
    def coprime(a, b):
        return gcd(a, b) == 1
    
    def random_int(N):
        return randint(-N, N)
    
    def random_det(N):
        a, b, c, d = randint(-N, N, 4)
        return a*d - b*c
    
    count = 0
    N = 10000000 # draw integers from [-N, N)
    num_reps = 1000000
    
    for _ in range(num_reps):
        count += coprime(random_int(N), random_int(N))
    print("Simulation: ", count/num_reps)
    print("Theory: ", 6*pi**-2)

This code printed

    Simulation:  0.607757
    Theory:  0.6079271018540267

when I ran it, so our simulation agreed with theory to three figures, the most you could expect from 10⁶ repetitions.

The analogous code for 2×2 matrices introduces a function random_det.

    def random_det(N):
        a, b, c, d = randint(-N, N, 4, dtype=int64)
        return a*d - b*c

I specified the dtype because the default is to use (32 bit) int as the type, which lead to Python complaining “RuntimeWarning: overflow encountered in long_scalars”.

I replaced random_int with random_det and reran the code above. This produced 0.452042. The exact value isn’t known in closed form, but we can see that it is between the bounds Δ(1) = 0.6079 and Δ(∞) = 0.3532.

Theory

In [1] the authors show that

This expression is only known to have a closed form when n = 1.

Related posts

• Probability of coprime sets
• Perpendicular and relatively prime

[1] Hafner, J. L.; Sarnak, P. & McCurley, K. (1993), “Relatively Prime Values of Polynomials”, in Knopp, M. & Seingorn, M. (eds.), A Tribute to Emil Grosswald: Number Theory and Related Analysis, Providence, RI: Amer. Math. Soc., ISBN 0-8218-5155-1.