Skip to content

Instantly share code, notes, and snippets.

# jcgregorio/gen-example.go

Created October 11, 2021 12:51
Show Gist options
• Save jcgregorio/0c2028a7a27b814641e9f502ed8b8b32 to your computer and use it in GitHub Desktop.
Random Integer Matrices With Inverses That Are Also Integer Matrices, aka Generating Unimodular Matrices
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 // 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)) }
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 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⎦
to join this conversation on GitHub. Already have an account? Sign in to comment