[system] / trunk / pg / lib / Value / Matrix.pm Repository: Repository Listing bbplugincoursesdistsnplrochestersystemwww

# Annotation of /trunk/pg/lib/Value/Matrix.pm

 1 : sh002i 2558 ########################################################################### 2 : # 3 : # Implements the Matrix class. 4 : # 5 : # @@@ Still needs lots of work @@@ 6 : # 7 : package Value::Matrix; 8 : my $pkg = 'Value::Matrix'; 9 : 10 : use strict; 11 : use vars qw(@ISA); 12 : @ISA = qw(Value); 13 : 14 : use overload 15 : '+' => \&add, 16 : '-' => \&sub, 17 : '*' => \&mult, 18 : '/' => \&div, 19 : '**' => \&power, 20 : '.' => \&Value::_dot, 21 : 'x' => \&Value::cross, 22 : '<=>' => \&compare, 23 : 'cmp' => \&compare, 24 : 'neg' => sub {$_[0]->neg}, 25 : 'nomethod' => \&Value::nomethod, 26 : '""' => \&stringify; 27 : 28 : # 29 : # Convert a value to a matrix. The value can be: 30 : # a list of numbers or list of (nested) references to arrays of numbers 31 : # a point, vector or matrix object 32 : # 33 : sub new { 34 : my $self = shift; my$class = ref($self) ||$self; 35 : my $M = shift; 36 : return bless {data =>$M->data}, $class 37 : if (Value::class($M) =~ m/Point|Vector|Matrix/ && scalar(@_) == 0); 38 : $M = [$M,@_] if ((defined($M) && ref($M) ne 'ARRAY') || scalar(@_) > 0); 39 : Value::Error("Matrices must have at least one entry") unless defined($M) && scalar(@{$M}) > 0; 40 : return $self->matrixMatrix(@{$M}) if (ref($M->[0]) eq 'ARRAY'); 41 : return$self->numberMatrix(@{$M}); 42 : } 43 : 44 : # 45 : # (Recusrively) make a matrix from a list of array refs 46 : # and report errors about the entry types 47 : # 48 : sub matrixMatrix { 49 : my$self = shift; my $class = ref($self) || $self; 50 : my ($x,$m); my @M = (); my$d = scalar(@{$_[0]}); my$isFormula = 0; 51 : foreach $x (@_) { 52 : if (Value::isFormula($x)) {push(@M,$x);$isFormula = 1} else { 53 : $m =$pkg->new($x); push(@M,$m); 54 : $isFormula = 1 if Value::isFormula($m); 55 : } 56 : } 57 : my ($type,$len) = ($M[0]->entryType->{name},$M[0]->length); 58 : foreach $x (@M) { 59 : Value::Error("Matrix rows must all be the same type") 60 : unless (defined($x->entryType) && $type eq$x->entryType->{name}); 61 : Value::Error("Matrix rows must all be the same length") unless ($len eq$x->length); 62 : } 63 : return $self->formula([@M]) if$isFormula; 64 : bless {data => [@M]}, $class; 65 : } 66 : 67 : # 68 : # Form a 1 x n matrix from a list of numbers 69 : # (could become a row of an m x n matrix) 70 : # 71 : sub numberMatrix { 72 : my$self = shift; my $class = ref($self) || $self; 73 : my @M = (); my$isFormula = 0; 74 : foreach my $x (@_) { 75 : Value::Error("Matrix row entries must be numbers") unless (Value::isNumber($x)); 76 : dpvc 2579 $x = Value::Real->make($x) if $$Value::context->flag('useFuzzyReals') && !Value::isFormula(x); 77 : sh002i 2558 push(@M,x); isFormula = 1 if Value::isFormula(x); 78 : } 79 : return self->formula([@M]) if isFormula; 80 : bless {data => [@M]}, class; 81 : } 82 : 83 : # 84 : # Recursively get the entries in the matrix and return 85 : # an array of (references to arrays of ... ) numbers 86 : # 87 : sub value { 88 : my self = shift; 89 : my M = self->data; 90 : return @{M} if Value::class(M->[0]) ne 'Matrix'; 91 : my @M = (); 92 : foreach my x (@{M}) {push(@M,[x->value])} 93 : return @M; 94 : } 95 : # 96 : # The number of rows in the matrix (for n x m) 97 : # or the number of entries in a 1 x n matrix 98 : # 99 : sub length {return scalar(@{shift->{data}})} 100 : # 101 : # Recursively get the dimensions of the matrix. 102 : # Returns (n) for a 1 x n, or (n,m) for an n x m, etc. 103 : # 104 : sub dimensions { 105 : my self = shift; 106 : my r = self->length; 107 : my v = self->data; 108 : return (r,) if (Value::class(v->[0]) ne 'Matrix'); 109 : return (r,v->[0]->dimensions); 110 : } 111 : # 112 : # Return the proper type for the matrix 113 : # 114 : sub typeRef { 115 : my self = shift; 116 : return Value::Type(self->class, self->length, Value::Type{number}) 117 : if (Value::class(self->data->[0]) ne 'Matrix'); 118 : return Value::Type(self->class, self->length, self->data->[0]->typeRef); 119 : } 120 : 121 : # 122 : # True if the matrix is a square matrix 123 : # 124 : sub isSquare { 125 : my self = shift; 126 : my @d = self->dimensions; 127 : return 0 if scalar(@d) > 2; 128 : return 1 if scalar(@d) == 1 && d[0] == 1; 129 : return d[0] == d[1]; 130 : } 131 : 132 : # 133 : # True if the matrix is 1-dimensional (i.e., is a matrix row) 134 : # 135 : sub isRow { 136 : my self = shift; 137 : my @d = self->dimensions; 138 : return scalar(@d) == 1; 139 : } 140 : 141 : # 142 : # Make arbitrary data into a matrix, if possible 143 : # 144 : sub promote { 145 : my x = shift; 146 : return pkg->new(x,@_) if scalar(@_) > 0 || ref(x) eq 'ARRAY'; 147 : return x if ref(x) eq pkg; 148 : return pkg->make(@{x->data}) if Value::class(x) =~ m/Point|Vector/; 149 : Value::Error("Can't convert ".Value::showClass(x)." to a Matrix"); 150 : } 151 : 152 : ############################################ 153 : # 154 : # Operations on matrices 155 : # 156 : 157 : sub add { 158 : my (l,r,flag) = @_; 159 : if (l->promotePrecedence(r)) {return r->add(l,!flag)} 160 : (l,r) = (promote(l)->data,promote(r)->data); 161 : Value::Error("Matrix addition with different dimensions") 162 : unless scalar(@{l}) == scalar(@{r}); 163 : my @s = (); 164 : foreach my i (0..scalar(@{l})-1) {push(@s,l->[i] + r->[i])} 165 : return pkg->make(@s); 166 : } 167 : 168 : sub sub { 169 : my (l,r,flag) = @_; 170 : if (l->promotePrecedence(r)) {return r->sub(l,!flag)} 171 : (l,r) = (promote(l)->data,promote(r)->data); 172 : Value::Error("Matrix subtraction with different dimensions") 173 : unless scalar(@{l}) == scalar(@{r}); 174 : if (flag) {my tmp = l; l = r; r = tmp}; 175 : my @s = (); 176 : foreach my i (0..scalar(@{l})-1) {push(@s,l->[i] - r->[i])} 177 : return pkg->make(@s); 178 : } 179 : 180 : sub mult { 181 : my (l,r,flag) = @_; 182 : if (l->promotePrecedence(r)) {return r->mult(l,!flag)} 183 : # 184 : # Constant multiplication 185 : # 186 : if (Value::matchNumber(r) || Value::isComplex(r)) { 187 : my @coords = (); 188 : foreach my x (@{l->data}) {push(@coords,x*r)} 189 : return pkg->make(@coords); 190 : } 191 : # 192 : # Make points and vectors into columns if they are on the right 193 : # 194 : if (!flag && Value::class(r) =~ m/Point|Vector/) 195 : {r = (promote(r))->transpose} else {r = promote(r)} 196 : # 197 : if (flag) {my tmp = l; l = r; r = tmp} 198 : my @dl = l->dimensions; my @dr = r->dimensions; 199 : if (scalar(@dl) == 1) {@dl = (1,@dl); l = pkg->make(l)} 200 : if (scalar(@dr) == 1) {@dr = (1,@dr); r = pkg->make(r)} 201 : Value::Error("Can only multiply 2-dimensional matrices") if scalar(@dl) > 2 || scalar(@dr) > 2; 202 : Value::Error("Matices of dimensions dl[0]xdl[1] and dr[0]xdr[1] can't be multiplied") 203 : unless (dl[1] == dr[0]); 204 : # 205 : # Do matrix multiplication 206 : # 207 : my @l = l->value; my @r = r->value; 208 : my @M = (); 209 : foreach my j (0..dr[1]-1) { 210 : my @row = (); 211 : foreach my i (0..dl[0]-1) { 212 : my s = 0; 213 : foreach my k (0..dl[1]-1) {s += l[i]->[k] * r[k]->[j]} 214 : push(@row,s); 215 : } 216 : push(@M,pkg->make(@row)); 217 : } 218 : return pkg->make(@M); 219 : } 220 : 221 : sub div { 222 : my (l,r,flag) = @_; 223 : if (l->promotePrecedence(r)) {return r->div(l,!flag)} 224 : Value::Error("Can't divide by a Matrix") if flag; 225 : Value::Error("Matrices can only be divided by numbers") 226 : unless (Value::matchNumber(r) || Value::isComplex(r)); 227 : Value::Error("Division by zero") if r == 0; 228 : my @coords = (); 229 : foreach my x (@{l->data}) {push(@coords,x/r)} 230 : return pkg->make(@coords); 231 : } 232 : 233 : sub power { 234 : my (l,r,flag) = @_; 235 : if (l->promotePrecedence(r)) {return r->power(l,!flag)} 236 : Value::Error("Can't use Matrices in exponents") if flag; 237 : Value::Error("Only square matrices can be raised to a power") unless l->isSquare; 238 : return Value::Matrix::I(l->length) if r == 0; 239 : Value::Error("Matrix powers must be positive integers") unless r =~ m/^[1-9]\d*/; 240 : my M = l; foreach my i (2..r) {M = M*l} 241 : return M; 242 : } 243 : 244 : # 245 : # Do lexicographic comparison 246 : # 247 : sub compare { 248 : my (l,r,flag) = @_; 249 : if (l->promotePrecedence(r)) {return r->compare(l,!flag)} 250 : (l,r) = (promote(l)->data,promote(r)->data); 251 : Value::Error("Matrix comparison with different dimensions") 252 : unless scalar(@{l}) == scalar(@{r}); 253 : if (flag) {my tmp = l; l = r; r = tmp}; 254 : my cmp = 0; 255 : foreach my i (0..scalar(@{l})-1) { 256 : cmp = l->[i] <=> r->[i]; 257 : last if cmp; 258 : } 259 : return cmp; 260 : } 261 : 262 : sub neg { 263 : my p = promote(@_)->data; 264 : my @coords = (); 265 : foreach my x (@{p}) {push(@coords,-x)} 266 : return pkg->make(@coords); 267 : } 268 : 269 : # 270 : # Transpose an n x m matrix 271 : # 272 : sub transpose { 273 : my self = shift; 274 : my @d = self->dimensions; 275 : if (scalar(@d) == 1) {@d = (1,@d); self = pkg->make(self)} 276 : Value::Error("Can't transpose ".scalar(@d)."-dimensional matrices") unless scalar(@d) == 2; 277 : my @M = (); my M = self->data; 278 : foreach my j (0..d[1]-1) { 279 : my @row = (); 280 : foreach my i (0..d[0]-1) {push(@row,M->[i]->data->[j])} 281 : push(@M,pkg->make(@row)); 282 : } 283 : return pkg->make(@M); 284 : } 285 : 286 : # 287 : # Get an identity matrix of the requested size 288 : # 289 : sub I { 290 : my d = shift; d = shift if ref(d) eq pkg; 291 : my @M = (); my @Z = split('',0 x d); 292 : foreach my i (0..d-1) { 293 : my @row = @Z; row[i] = 1; 294 : push(@M,pkg->make(@row)); 295 : } 296 : return pkg->make(@M); 297 : } 298 : 299 : # 300 : # Extract a given row from the matrix 301 : # 302 : sub row { 303 : my M = promote(shift); my i = shift; 304 : return if i == 0; i-- if i > 0; 305 : if (M->isRow) {return if i != 0; return M} 306 : return M->data->[i]; 307 : } 308 : 309 : # 310 : # Extract a given element from the matrix 311 : # 312 : sub element { 313 : my M = promote(shift); 314 : return M->extract(@_); 315 : } 316 : 317 : # 318 : # Extract a given column from the matrix 319 : # 320 : sub column { 321 : my M = promote(shift); my j = shift; 322 : return if j == 0; j-- if j > 0; 323 : my @d = M->dimensions; my @col = (); 324 : return if j+1 > d[1]; 325 : return M->data->[j] if scalar(@d) == 1; 326 : foreach my row (@{M->data}) {push(@col,pkg->make(row->data->[j]))} 327 : return pkg->make(@col); 328 : } 329 : 330 : # @@@ removeRow, removeColumn @@@ 331 : # @@@ Det, inverse @@@ 332 : 333 : ############################################ 334 : # 335 : # Generate the various output formats 336 : # 337 : 338 : sub stringify { 339 : my self = shift; 340 : dpvc 2579 my open =$$Value::context->lists->get('Matrix')->{open}; 341 : my $close = $$Value::context->lists->get('Matrix')->{close}; 342 : sh002i 2558 return open.join(',',@{self->data}).close 343 : if (Value::class(self->data->[0]) ne 'Matrix'); 344 : return open.join(",\n ",@{self->data}).close; 345 : } 346 : 347 : sub string { 348 : my self = shift; my equation = shift; 349 : dpvc 2579 my open = shift ||$$Value::context->lists->get('Matrix')->{open}; 350 : my$close = shift || $$Value::context->lists->get('Matrix')->{close}; 351 : sh002i 2558 my @coords = (); 352 : foreach my x (@{self->data}) { 353 : if (Value::isValue(x)) {push(@coords,x->string(equation,open,close))} 354 : else {push(@coords,x)} 355 : } 356 : return open.join(',',@coords).close; 357 : } 358 : 359 : # 360 : # Use \matrix to lay out matrices 361 : # 362 : sub TeX { 363 : my self = shift; my equation = shift; 364 : dpvc 2579 my open = shift ||$$Value::context->lists->get('Matrix')->{open}; 365 : my $close = shift ||$$Value::context->lists->get('Matrix')->{close}; 366 : sh002i 2558$open = '\{' if $open eq '{';$close = '\}' if $close eq '}'; 367 : my$TeX = ''; my @entries = (); 368 : dpvc 2579 if ($self->isRow) { 369 : foreach my$x (@{$self->data}) { 370 : sh002i 2558 push(@entries,(Value::isValue($x))? $x->TeX($equation,$open,$close): $x); 371 : } 372 : dpvc 2579$TeX .= join(' &',@entries) . "\n"; 373 : } else { 374 : foreach my $row (@{$self->data}) { 375 : foreach my $x (@{$row->data}) { 376 : push(@entries,(Value::isValue($x))?$x->TeX($equation,$open,$close):$x); 377 : } 378 : $TeX .= join(' &',@entries) . '\cr'."\n"; @entries = (); 379 : } 380 : sh002i 2558 } 381 : return '\left'.$open.'\matrix{'."\n".$TeX.'}\right'.$close; 382 : } 383 : 384 : ########################################################################### 385 : 386 : 1; 387 :