# At the next prime, turn left

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 . 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

 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.

# Sum of divisor powers

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

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

A couple special cases may use different notation.

• σ1(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.

# An almost integer pattern in many bases

A few days ago I mentioned a passing comment in video by Richard Boucherds. This post picks up on another off-hand remark from that post.

Boucherds was discussing why exp(π √67) and exp(π √163) are nearly an integer.

exp(π √67) = 147197952744 – ε1
exp(π √163) = 262537412640768744 – ε2

where ε1 and ε2 and on the order of 10-6 and 10-12 respectively.

He called attention to the 744 at the end and commented that this isn’t just an artifact of base 10. In many other bases, these numbers end in that bases’ representation of 744. This is what I wanted to expand on.

To illustrate Boucherds’ remark in hexadecimal, note that

    147197952744 -> 0x2245ae82e8
262537412640768744 -> 0x3a4b862c4b402e8
744 -> 0x2e8


Boucherds’ comment is equivalent to saying

147197952744 = 744 mod m

and

262537412640768744 = 744 mod m

for many values of m. Equivalently 147197952000 and 262537412640768000 have a lot of factors; every factor of these numbers is a base where the congruence holds.

So for how many bases m are these numbers congruent to 744?

The number of factors of a number n is written d(n). This is a multiplicative function, meaning that for relatively prime numbers a and b,

d(ab) = d(a) d(b).

Note that

147197952000 = 215 33 53 113

and

262537412640768000 = 218 33 53 233 293

It follows that

d(147197952000) =
d(215 33 53 113) =
d(215) d(33) d(53) d(113).

Now for any prime power pk

d(pk) = k + 1,

and so

d(147197952000) = 16 × 4 × 4 × 4.

Similarly

d(262537412640768000) = 19 × 4 × 4 × 4 × 4.

For more about almost integers, watch Boucherds’ video and see this post.

# Continued fractions with period 1

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  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. Samuel S. Wagstaff, Jr. The Joy of Factoring. Theorem 6.15.

# A bevy of ones

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 .

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 11111sixteen is a multiple of 5, and in fact

5 × 369dsixteen = 11111sixteen.

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

# Binomial coefficients mod primes

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  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

 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.

# Bit flipping to primes

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 . 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?

 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.

# Square roots of Gaussian integers

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 zxiy 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

(2mnm² – 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 xiy is a Gaussian integer.

# How probable is a probable prime?

A probable prime is a number that passes a test that all primes pass and that most composite numbers fail. Specifically, a Fermat probable prime is a number that passes Fermat’s primality test. Fermat’s test is the most commonly used, so that’s nearly always what anyone means by probable prime unless they’re more specific.

A number n is a Fermat probable prime for base b if

bn-1 = 1 (mod n).

This test isn’t conclusive, but it can be implemented efficiently and it weeds out most composite numbers. To read more on probable primes, see this post.

If a number is a probable prime, how probable is it that it really is prime? This post will briefly summarize some results from a paper  that makes this precise. From that paper:

… let P(x) denote the probability that an integer n is composite given that

1. n is chosen at random with 1 < nx, n odd,
2. b is chosen at random with 1 < b < n − 1, and
3. n is a probable prime to the base b.

The authors give upper bounds on P(x) for x equal to various powers of 2. In particular they report

P(2250) ≤ 5.876 × 10−6

and

P(2260) ≤ 3.412 × 10−6

and so the chances that a 256-bit probable prime is composite are in the neighborhood of 4 in a million. In practice, one would test with multiple b‘s. The tests for various b‘s aren’t entirely independent, but running the tests for multiple bases does mean that fewer composite numbers slip through. There are a few pesky numbers, the Carmichael numbers, that are Fermat probable primes for nearly all bases (see footnote  here for more details), but these are rare.

I looked through the paper for results for larger powers of 2 to get results that would be applicable to, for example, RSA keys. The largest result explicit in the paper is

P(2330) ≤ 8.713 × 10−8

though I would like to know P(21024) and P(21536) since RSA keys are the products of (probable) primes this large. Presumably the results in  could be used to compute P(x) for larger values of x, but I haven’t read the paper closely enough to how much personal effort or computational effort that would require. I imagine it would be difficult or else the authors would have reported results for probable primes of the size frequently used in applications.

## Related posts

 Jared Ducker Lichtman and Carl Pomerance. Improved error bounds for the Fermat primality test on random inputs. Mathematics of Computation. Volume 87, Number 314, November 2018, pages 2871–2890.

# Relatively prime determinants

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 , 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 106 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  the authors show that This expression is only known to have a closed form when n = 1.

## Related posts

 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.