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

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

Sun Sep 19 14:27:39 2004 UTC (15 years, 5 months ago) by dpvc
File size: 11496 byte(s)
Added isZero and isOne checks for Parser::Value objects (i.e., for
constants within formulas).  These now correctly handle vector and
matrices, in particular.  The isOne and isZero checks are used in the
reduce() method to simplify formulas.


    1 ###########################################################################
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' => \&Value::cmp, 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->numberMatrix(@{$M}) if Value::isNumber($M->[0]); 41 return$self->matrixMatrix(@{$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$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     $x = Value::Real->make($x) if !ref($x); 77 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 # See if the matrix is an Indenity matrix 143 # 144 sub isOne { 145 my$self = shift;
146   return 0 unless $self->isSquare; 147 my$i = 0;
148   foreach my $row (@{$self->{data}}) {
149     my $j = 0; 150 foreach my$k (@{$row->{data}}) { 151 return 0 unless$k eq (($i ==$j)? "1": "0");
152       $j++; 153 } 154$i++;
155   }
156   return 1;
157 }
158
159 #
160 #  See if the matrix is all zeros
161 #
162 sub isZero {
163   my $self = shift; 164 foreach my$x (@{$self->{data}}) {return 0 unless$x->isZero}
165   return 1;
166 }
167
168 #
169 #  Make arbitrary data into a matrix, if possible
170 #
171 sub promote {
172   my $x = shift; 173 return$pkg->new($x,@_) if scalar(@_) > 0 || ref($x) eq 'ARRAY';
174   return $x if ref($x) eq $pkg; 175 return$pkg->make(@{$x->data}) if Value::class($x) =~ m/Point|Vector/;
176   Value::Error("Can't convert ".Value::showClass($x)." to a Matrix"); 177 } 178 179 ############################################ 180 # 181 # Operations on matrices 182 # 183 184 sub add { 185 my ($l,$r,$flag) = @_;
186   if ($l->promotePrecedence($r)) {return $r->add($l,!$flag)} 187 ($l,$r) = (promote($l)->data,promote($r)->data); 188 Value::Error("Matrix addition with different dimensions") 189 unless scalar(@{$l}) == scalar(@{$r}); 190 my @s = (); 191 foreach my$i (0..scalar(@{$l})-1) {push(@s,$l->[$i] +$r->[$i])} 192 return$pkg->make(@s);
193 }
194
195 sub sub {
196   my ($l,$r,$flag) = @_; 197 if ($l->promotePrecedence($r)) {return$r->sub($l,!$flag)}
198   ($l,$r) = (promote($l)->data,promote($r)->data);
199   Value::Error("Matrix subtraction with different dimensions")
200     unless scalar(@{$l}) == scalar(@{$r});
201   if ($flag) {my$tmp = $l;$l = $r;$r = $tmp}; 202 my @s = (); 203 foreach my$i (0..scalar(@{$l})-1) {push(@s,$l->[$i] -$r->[$i])} 204 return$pkg->make(@s);
205 }
206
207 sub mult {
208   my ($l,$r,$flag) = @_; 209 if ($l->promotePrecedence($r)) {return$r->mult($l,!$flag)}
210   #
211   #  Constant multiplication
212   #
213   if (Value::matchNumber($r) || Value::isComplex($r)) {
214     my @coords = ();
215     foreach my $x (@{$l->data}) {push(@coords,$x*$r)}
216     return $pkg->make(@coords); 217 } 218 # 219 # Make points and vectors into columns if they are on the right 220 # 221 if (!$flag && Value::class($r) =~ m/Point|Vector/) 222 {$r = (promote($r))->transpose} else {$r = promote($r)} 223 # 224 if ($flag) {my $tmp =$l; $l =$r; $r =$tmp}
225   my @dl = $l->dimensions; my @dr =$r->dimensions;
226   if (scalar(@dl) == 1) {@dl = (1,@dl); $l =$pkg->make($l)} 227 if (scalar(@dr) == 1) {@dr = (@dr,1);$r = $pkg->make($r)->transpose}
228   Value::Error("Can only multiply 2-dimensional matrices") if scalar(@dl) > 2 || scalar(@dr) > 2;
229   Value::Error("Matices of dimensions $dl[0]x$dl[1] and $dr[0]x$dr[1] can't be multiplied")
230     unless ($dl[1] ==$dr[0]);
231   #
232   #  Do matrix multiplication
233   #
234   my @l = $l->value; my @r =$r->value;
235   my @M = ();
236   foreach my $j (0..$dr[1]-1) {
237     my @row = ();
238     foreach my $i (0..$dl[0]-1) {
239       my $s = 0; 240 foreach my$k (0..$dl[1]-1) {$s += $l[$i]->[$k] *$r[$k]->[$j]}
241       push(@row,$s); 242 } 243 push(@M,$pkg->make(@row));
244   }
245   return $pkg->make(@M); 246 } 247 248 sub div { 249 my ($l,$r,$flag) = @_;
250   if ($l->promotePrecedence($r)) {return $r->div($l,!$flag)} 251 Value::Error("Can't divide by a Matrix") if$flag;
252   Value::Error("Matrices can only be divided by numbers")
253     unless (Value::matchNumber($r) || Value::isComplex($r));
254   Value::Error("Division by zero") if $r == 0; 255 my @coords = (); 256 foreach my$x (@{$l->data}) {push(@coords,$x/$r)} 257 return$pkg->make(@coords);
258 }
259
260 sub power {
261   my ($l,$r,$flag) = @_; 262 if ($l->promotePrecedence($r)) {return$r->power($l,!$flag)}
263   Value::Error("Can't use Matrices in exponents") if $flag; 264 Value::Error("Only square matrices can be raised to a power") unless$l->isSquare;
265   return Value::Matrix::I($l->length) if$r == 0;
266   Value::Error("Matrix powers must be positive integers") unless $r =~ m/^[1-9]\d*$/;
267   my $M =$l; foreach my $i (2..$r) {$M =$M*$l} 268 return$M;
269 }
270
271 #
272 #  Do lexicographic comparison
273 #
274 sub compare {
275   my ($l,$r,$flag) = @_; 276 if ($l->promotePrecedence($r)) {return$r->compare($l,!$flag)}
277   ($l,$r) = (promote($l)->data,promote($r)->data);
278   Value::Error("Matrix comparison with different dimensions")
279     unless scalar(@{$l}) == scalar(@{$r});
280   if ($flag) {my$tmp = $l;$l = $r;$r = $tmp}; 281 my$cmp = 0;
282   foreach my $i (0..scalar(@{$l})-1) {
283     $cmp =$l->[$i] <=>$r->[$i]; 284 last if$cmp;
285   }
286   return $cmp; 287 } 288 289 sub neg { 290 my$p = promote(@_)->data;
291   my @coords = ();
292   foreach my $x (@{$p}) {push(@coords,-$x)} 293 return$pkg->make(@coords);
294 }
295
296 #
297 #  Transpose an  n x m  matrix
298 #
299 sub transpose {
300   my $self = shift; 301 my @d =$self->dimensions;
302   if (scalar(@d) == 1) {@d = (1,@d); $self =$pkg->make($self)} 303 Value::Error("Can't transpose ".scalar(@d)."-dimensional matrices") unless scalar(@d) == 2; 304 my @M = (); my$M = $self->data; 305 foreach my$j (0..$d[1]-1) { 306 my @row = (); 307 foreach my$i (0..$d[0]-1) {push(@row,$M->[$i]->data->[$j])}
308     push(@M,$pkg->make(@row)); 309 } 310 return$pkg->make(@M);
311 }
312
313 #
314 #  Get an identity matrix of the requested size
315 #
316 sub I {
317   my $d = shift;$d = shift if ref($d) eq$pkg;
318   my @M = (); my @Z = split('',0 x $d); 319 foreach my$i (0..$d-1) { 320 my @row = @Z;$row[$i] = 1; 321 push(@M,$pkg->make(@row));
322   }
323   return $pkg->make(@M); 324 } 325 326 # 327 # Extract a given row from the matrix 328 # 329 sub row { 330 my$M = promote(shift); my $i = shift; 331 return if$i == 0; $i-- if$i > 0;
332   if ($M->isRow) {return if$i != 0; return $M} 333 return$M->data->[$i]; 334 } 335 336 # 337 # Extract a given element from the matrix 338 # 339 sub element { 340 my$M = promote(shift);
341   return $M->extract(@_); 342 } 343 344 # 345 # Extract a given column from the matrix 346 # 347 sub column { 348 my$M = promote(shift); my $j = shift; 349 return if$j == 0; $j-- if$j > 0;
350   my @d = $M->dimensions; my @col = (); 351 return if$j+1 > $d[1]; 352 return$M->data->[$j] if scalar(@d) == 1; 353 foreach my$row (@{$M->data}) {push(@col,$pkg->make($row->data->[$j]))}
354   return $pkg->make(@col); 355 } 356 357 # @@@ removeRow, removeColumn @@@ 358 # @@@ Det, inverse @@@ 359 360 ############################################ 361 # 362 # Generate the various output formats 363 # 364 365 sub stringify { 366 my$self = shift;
367   return $self->TeX(undef,$self->{open},$self->{close}) if $$Value::context->flag('StringifyAsTeX'); 368 return self->string(undef,self->{open},self->{close}) 369 } 370 371 sub string { 372 my self = shift; my equation = shift; 373 my def = (equation->{context} ||$$Value::context)->lists->get('Matrix'); 374 my$open = shift || $def->{open}; my$close = shift || $def->{close}; 375 my @coords = (); 376 foreach my$x (@{$self->data}) { 377 if (Value::isValue($x)) {push(@coords,$x->string($equation,$open,$close))}
378       else {push(@coords,$x)} 379 } 380 return$open.join(',',@coords).$close; 381 } 382 383 # 384 # Use \matrix to lay out matrices 385 # 386 sub TeX { 387 my$self = shift; my $equation = shift; 388 my$def = ($equation->{context} ||$$Value::context)->lists->get('Matrix'); 389 my$open = shift || $def->{open}; my$close = shift || $def->{close}; 390$open = '\{' if $open eq '{';$close = '\}' if $close eq '}'; 391 my$TeX = ''; my @entries = (); my $d; 392 if ($self->isRow) {
393     foreach my $x (@{$self->data}) {
394       push(@entries,(Value::isValue($x))?$x->TeX($equation):$x);
395     }
396     $TeX .= join(' &',@entries) . "\n"; 397$d = scalar(@entries);
398   } else {
399     foreach my $row (@{$self->data}) {
400       foreach my $x (@{$row->data}) {
401         push(@entries,(Value::isValue($x))?$x->TeX($equation):$x);
402       }
403       $TeX .= join(' &',@entries) . '\cr'."\n"; 404$d = scalar(@entries); @entries = ();
405     }
406   }
407   return '\left'.$open.'\begin{array}{'.('c'x$d).'}'."\n".$TeX.'\end{array}\right'.$close;
408 }
409
410 ###########################################################################
411
412 1;
413