Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Random Integer Matrices With Inverses That Are Also Integer Matrices, aka Generating Unimodular Matrices
// 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))
}
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⎦
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment