# matrix.tcl - matrix math routines # Jonas Beskow 2000 # the internal representation of a matrix is an array structured as below: # m(rows) - number of rows # m(cols) - number of columns # m(0,0) m(0,1) ... m($i,$j) - elements of the matrix # # matrices are passed as lists (i.e. array get/set is used to # serialize/restore the matrix array) namespace eval matrix { } # mult - multiply two matrices proc matrix::mult {m1l m2l} { array set m1 $m1l array set m2 $m2l if {$m1(cols) != $m2(rows)} { error "matrix inner dimensions must agree" } for {set i 0} {$i<$m1(rows)} {incr i} { for {set j 0} {$j<$m2(cols)} {incr j} { set iprod 0 for {set k 0} {$k<$m1(cols)} {incr k} { set iprod [expr {$iprod + $m1($i,$k)*$m2($k,$j)}] } set m($i,$j) $iprod } } set m(rows) $m1(rows) set m(cols) $m2(cols) array get m } # transpose - transpose a matrix proc matrix::transpose {ml} { array set m $ml for {set i 0} {$i<$m(rows)} {incr i} { for {set j 0} {$j<$m(cols)} {incr j} { set n($j,$i) $m($i,$j) } } set n(rows) $m(cols) set n(cols) $m(rows) array get n } # eye - generate a unit matrix of a given dimension proc matrix::eye {dim} { set m(rows) $dim set m(cols) $dim for {set i 0} {$i<$dim} {incr i} { for {set j 0} {$j<$dim} {incr j} { if {$i==$j} {set m($i,$j) 1.0} {set m($i,$j) 0.0} } } array get m } # det - calculate the determinant of a matrix proc matrix::det {ml} { array set m $ml if {$m(rows)!=$m(cols)} {error "matrix must be square"} set n $m(rows) set det 0 for {set i 0} {$i<$n} {incr i} {lappend colperm0 $i} foreach {colperm alpha} [_permute2 $colperm0] { set row 0 set prod 1 foreach col $colperm { set prod [expr {$prod*$m($row,$col)}] incr row } if {$alpha%2==0} {set sign 1} {set sign -1} set det [expr {$det+$sign*$prod}] } return $det } # det2 - calculate the determinant of a matrix (recursive algorithm) # algorithm: # Choose a fixed row value i. # The determinant can be calculated emanating from the ith row. # |A| = Ai,1 * ai,1 + Ai,2 * ai,2 + Ai,3 * ai,3 + ... Ai,n * ai,n # The value of Ai,j = (-1)^(i+j) * (the determinant of the sub-matrix of A, # obtained from A by crossing out the ith row and the jth column) proc matrix::det2 {ml} { array set m $ml set self [lindex [info level 0] 0] if {$m(rows)!=$m(cols)} {error "matrix must be square"} set n $m(rows) if {$n==1} {return $m(0,0)} set i 0 set det 0 for {set j 0} {$j<$n} {incr j} { if {($i+$j)%2==0} {set sign 1} {set sign -1} set subm [_wipe $ml $i $j] set subdet [$self $subm] set det [expr {$det + $sign * $subdet * $m($i,$j)}] } return $det } # lset - convert from list to internal representation proc matrix::lset {list} { set i 0 foreach row $list { puts row=$row set j 0 foreach elem $row { set m($i,$j) $elem incr j } incr i } set m(rows) $i set m(cols) $j array get m } # lget - convert from internal to list representation proc matrix::lget {ml} { array set m $ml set list "" for {set i 0} {$i<$m(rows)} {incr i} { set row {} for {set j 0} {$j<$m(cols)} {incr j} { lappend row $m($i,$j) } puts row=$row lappend list $row } set list } # cset - convert from string to internal representation proc matrix::cset {s} { set list "" regsub -all {([[:blank:]])*(\n)+([[:blank:]])*} [string trim $s] \n s puts s=$s regsub -all {[[:blank:]]+} $s " " s puts s=$s foreach row [split $s \n] { lappend list [split $row " "] } lset $list } # lget - convert from internal to string representation proc matrix::cget {ml} { array set m $ml set str "" for {set i 0} {$i<$m(rows)} {incr i} { set row {} for {set j 0} {$j<$m(cols)} {incr j} { append str $m($i,$j) if {$j!=[expr $m(cols)-1]} {append str \t} } if {$i!=[expr $m(rows)-1]} {append str \n} } string trim $str } # # utility routines to generate common 4x4 matrices for 3D transformation # proc matrix::translate {x y z} { lset [list \ "1 0 0 $x" \ "0 1 0 $y" \ "0 0 1 $z" \ "0 0 0 1"] } proc matrix::scale {x y z} { lset [list \ "$x 0 0 0" \ "0 $y 0 0" \ "0 0 $z 0" \ "0 0 0 1"] } proc matrix::rotX {a} { set ca [expr {cos($a)}] set sa [expr {sin($a)}] set nsa [expr {-$sa}] lset [list \ "1 0 0 0" \ "0 $ca $nsa 0" \ "0 $sa $ca 0" \ "0 0 0 1"] } proc matrix::rotY {a} { set ca [expr {cos($a)}] set sa [expr {sin($a)}] set nsa [expr {-$sa}] lset [list \ "$ca 0 $sa 0" \ "0 1 0 0" \ "$nsa 0 $ca 0" \ "0 0 0 1"] } proc matrix::rotZ {a} { set ca [expr {cos($a)}] set sa [expr {sin($a)}] set nsa [expr {-$sa}] lset [list \ "$ca $nsa 0 0" \ "$sa $ca 0 0" \ "0 0 1 0" \ "0 0 0 1"] } # _permute2 # - recursive help function to calculate all possible permutations of a list # returns a list of the form {perm1 n1 perm2 n2 ... permN nN} where # permK is a sub-list containing a permutation and nK is the number of # exchanges necessary to bring the original sequence to permK proc matrix::_permute2 {list} { set self [lindex [info level 0] 0] if {[llength $list]==1} {return [list $list 0]} set result "" set first [lindex $list 0] set tail [lrange $list 1 end] foreach {p n} [$self $tail] { for {set i 0} {$i<=[llength $p]} {incr i} { lappend result [linsert $p $i $first] [expr $i+$n] } } return $result } # _wipe - help func to det2 proc matrix::_wipe {ml rr cc} { array set m $ml set m2(rows) [expr {$m(rows)-1}] set m2(cols) [expr {$m(cols)-1}] for {set r 0} {$r<$m(rows)} {incr r} { if {$r==$rr} continue if {$r>$rr} {set r2 [expr {$r-1}]} else {set r2 $r} for {set c 0} {$c<$m(cols)} {incr c} { if {$c==$cc} continue if {$c>$cc} {set c2 [expr {$c-1}]} else {set c2 $c} set m2($r2,$c2) $m($r,$c) } } array get m2 } return # Testing below # ----------------------------------------------------------------------------- # # magic squares set m10 [m::cset " 92 99 1 8 15 67 74 51 58 40 98 80 7 14 16 73 55 57 64 41 4 81 88 20 22 54 56 63 70 47 85 87 19 21 3 60 62 69 71 28 86 93 25 2 9 61 68 75 52 34 17 24 76 83 90 42 49 26 33 65 23 5 82 89 91 48 30 32 39 66 79 6 13 95 97 29 31 38 45 72 10 12 94 96 78 35 37 44 46 53 11 18 100 77 84 36 43 50 27 59"] set m7 [m::cset " 30 39 48 1 10 19 28 38 47 7 9 18 27 29 46 6 8 17 26 35 37 5 14 16 25 34 36 45 13 15 24 33 42 44 4 21 23 32 41 43 3 12 22 31 40 49 2 11 20"] set m5 [m::cset " 17 24 1 8 15 23 5 7 14 16 4 6 13 20 22 10 12 19 21 3 11 18 25 2 9"] set m3 [m::cset " 8 1 6 3 5 7 4 9 2"] set a [m::lset [list {1 2} {3 4}]] set a [m::lset [list {1 2 3 4} {5 6 7 8} {9 10 11 12} {13 14 15 16}]] set b [m::lset [list {9 2 1 6} {6 4 7 23} {23 1 3 6} {98 5 3 1}]] puts [m::lget [m::mult $a $b]] set v [m::lset {1 0 0 1}] set T [m::translate 0 .5 0] set R [m::rotZ [expr 2*atan(1)]] puts [m::lget [m::mult $R $v]]