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.