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

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

Tue Feb 15 21:53:23 2005 UTC (15 years ago) by dpvc
File size: 11659 byte(s)
Improved the Real(), Complex(), Point(), Vector(), Matrix() and
String() constructors so that they will process formulas passed to
them as strings rather than requiring perl objects for these.

For example, you can use Real("2/3") rather than Real(2/3) if you
want.  Also, Real("1+x") will return a formula returning a real
(essentially the same as Formula("1+x") in this case).


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