Let σ(*n*) be the sum of the positive divisors of *n* and let gcd(*a*, *b*) be the greatest common divisor of *a* and *b*.

Form an *n* by *n* matrix *M* whose (*i,* *j*) entry is σ(gcd(*i*, *j*)). Then the determinant of *M* is *n*!.

The following code shows that the theorem is true for a few values of *n* and shows how to do some common number theory calculations in SymPy.

from sympy import gcd, divisors, Matrix, factorial def f(i, j): return sum( divisors( gcd(i, j) ) ) def test(n): r = range(1, n+1) M = Matrix( [ [f(i, j) for j in r] for i in r] ) return M.det() - factorial(n) for n in range(1, 11): print( test(n) )

As expected, the test function returns zeros.

If we replace the function σ above by τ where τ(*n*) is the number of positive divisors of *n*, the corresponding determinant is 1. To test this, replace `sum`

by `len`

in the definition of `f`

and replace `factorial(n)`

by 1.

In case you’re curious, both results are special cases of the following more general theorem. I don’t know whose theorem it is. I found it here.

For any arithmetic function *f*(*m*), let *g*(*m*) be defined for all positive integers *m* by

Let *M* be the square matrix of order *n* with *ij* element *f*(gcd(*i*, *j*)). Then

Here μ is the Möbius function. The two special cases above correspond to *g*(*m*) = *m* and *g*(*m*) = 1.

It’s better to use igcd() to get the gcd of integers. gcd() does polynomial gcd.