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

$\Delta(n) = \prod_{p \text{ prime}}\left[ 1- \left(1 - \prod_{k=1}^n (1-p^{-k})\right )^2 \right ]$

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

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

Python simulation

Theory

Related posts