Random Integer Matrices With Inverses That Are Also Integer Matrices, aka Generating Unimodular Matrices

gen-example.go

// Demo application for generating random Unimodular Matrices: https://en.wikipedia.org/wiki/Unimodular_matrix
package main

import (
	"fmt"
	"math/rand"
	"time"

	"gonum.org/v1/gonum/mat"
)

// eye returns a new identity matrix of size n×n.
func eye(n int) *mat.Dense {
	ret := mat.NewDense(n, n, nil)
	for i := 0; i < n; i++ {
		ret.Set(i, i, 1)
	}
	return ret
}

// randUnimodular returns a random nxn matrix A this has all integer
// coefficients >=0, and it's inverse, which is comprised of integers.
//
// permutationLoops is the number of times the function runs through a
// permutation of the rows and applies row operations for each.
func randUnimodular(n int, permutationLoops int) (*mat.Dense, *mat.Dense) {
	A := eye(n)
	AInv := eye(n)
	// You can keep applying as many row operations as you like, but n/2
	// permutations seems to produce good matricies with numbers that aren't
	// huge, and very few zeroes.
	for loops := 0; loops < permutationLoops; loops++ {
		// We want a random matrix, but making sure that we apply row operations
		// to every row gives nicer results, so we generate a permutation of all
		// the rows.
		for x, y := range rand.Perm(n) {
			// Don't use a permutation if x == y, since that row operation will
			// actually change the determinant of the matrix.
			if x == y {
				continue
			}
			for index := 0; index < n; index++ {
				// Apply the row operation to A.
				A.Set(y, index, A.At(y, index)+A.At(x, index))

				// Apply the inverse column operation to A^{-1}.
				AInv.Set(index, x, AInv.At(index, x)-AInv.At(index, y))
			}
		}
	}
	return A, AInv
}

func main() {
	rand.Seed(time.Now().UnixNano())

	A, AInv := randUnimodular(8, 5)
	fmt.Printf("A\n%v\n\n", mat.Formatted(A))
	fmt.Printf("A^{-1}\n%v\n\n", mat.Formatted(AInv))
	var i mat.Dense
	(&i).Mul(A, AInv)
	fmt.Printf("A A^{-1} = I\n%v\n", mat.Formatted(&i))
}

output.txt

A
⎡ 9  21  23   8  18   7  12   7⎤
⎢ 8  13  15   5  10   5   9   6⎥
⎢ 0  13  11   2  13   1   3   0⎥
⎢13  11  16  13   6  11  13  10⎥
⎢11  21  23   7  17   7  13   8⎥
⎢10   8  12   9   4   8  10   8⎥
⎢11  16  19   7  12   7  12   8⎥
⎣ 6  12  14   6  10   5   8   5⎦

A^{-1}
⎡ 0   2  -1   2   1  -3  -2   0⎤
⎢-1   1   2  -1  -2   2   2  -1⎥
⎢ 2   1   1  -1  -5   1   4  -2⎥
⎢-2   4   0   2   0  -3  -2   2⎥
⎢ 0  -1  -3   2   6  -3  -6   2⎥
⎢ 4  -7   0  -3   0   5   3  -4⎥
⎢-3  -4   1  -1   1   1   3   3⎥
⎣ 0   3  -2   1   4  -1  -7   1⎦

A A^{-1} = I
⎡1  0  0  0  0  0  0  0⎤
⎢0  1  0  0  0  0  0  0⎥
⎢0  0  1  0  0  0  0  0⎥
⎢0  0  0  1  0  0  0  0⎥
⎢0  0  0  0  1  0  0  0⎥
⎢0  0  0  0  0  1  0  0⎥
⎢0  0  0  0  0  0  1  0⎥
⎣0  0  0  0  0  0  0  1⎦