[system] / trunk / pg / lib / Value / Complex.pm Repository:
ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2678 - (download) (as text) (annotate)
Mon Aug 23 23:55:37 2004 UTC (15 years, 6 months ago) by dpvc
File size: 9416 byte(s)
Modified the parser so that the classes for the various object
constructors are stored in the context table rather than hard-coded
into the parser.  That way, you can override the default classes with
your own.  This gives you even more complete control to modify the
parser.  (You had been able to replace the definitions of operators,
functions and list-like objects, but could not override the behaviour
of numbers, strings, variables, and so on.  Now you can.)

This effects most of the files, but only by changing the name of the
calls that create the various objects.

There are also a couple of other minor fixes.

    1 ##########################################################################
    2 
    3 package Value::Complex;
    4 my $pkg = 'Value::Complex';
    5 
    6 use strict;
    7 use vars qw(@ISA $i $pi);
    8 @ISA = qw(Value);
    9 
   10 use overload
   11        '+'   => \&add,
   12        '-'   => \&sub,
   13        '*'   => \&mult,
   14        '/'   => \&div,
   15        '**'  => \&power,
   16        '.'   => \&Value::_dot,
   17        'x'   => \&Value::cross,
   18        '<=>' => \&compare,
   19        'cmp' => \&Value::cmp,
   20        '~'   => sub {$_[0]->conj},
   21        'neg' => sub {$_[0]->neg},
   22        'abs' => sub {$_[0]->norm},
   23        'sqrt'=> sub {$_[0]->sqrt},
   24        'exp' => sub {$_[0]->exp},
   25        'log' => sub {$_[0]->log},
   26        'sin' => sub {$_[0]->sin},
   27        'cos' => sub {$_[0]->cos},
   28      'atan2' => \&atan2,
   29   'nomethod' => \&Value::nomethod,
   30         '""' => \&Value::stringify;
   31 
   32 #
   33 #  Check that the inputs are:
   34 #    one or two real numbers, or
   35 #    an array ref of one or two reals, or
   36 #    a Value::Complex object
   37 #  Make a formula if either part is a formula
   38 #
   39 sub new {
   40   my $self = shift; my $class = ref($self) || $self;
   41   my $x = shift; $x = [$x,@_] if scalar(@_) > 0;
   42   $x = $x->data if ref($x) eq $pkg || Value::isReal($x);
   43   $x = [$x] unless ref($x) eq 'ARRAY'; $x->[1] = 0 unless defined($x->[1]);
   44   Value::Error("Can't convert ARRAY of length ".scalar(@{$x})." to a Complex Number")
   45     unless (scalar(@{$x}) == 2);
   46   $x->[0] = Value::makeValue($x->[0]); $x->[1] = Value::makeValue($x->[1]);
   47   Value::Error("Real part can't be ".Value::showClass($x->[0]))
   48      unless (Value::isRealNumber($x->[0]));
   49   Value::Error("Imaginary part can't be ".Value::showClass($x->[1]))
   50      unless (Value::isRealNumber($x->[1]));
   51   return $self->formula($x) if Value::isFormula($x->[0]) || Value::isFormula($x->[1]);
   52   bless {data => $x}, $class;
   53 }
   54 
   55 #
   56 #  Create a new a+b*i formula from the two parts
   57 #
   58 sub formula {
   59   my $self = shift; my $value = shift;
   60   my $formula = Value::Formula->blank;
   61   my ($l,$r) = Value::toFormula($formula,@{$value});
   62   my $parser = $formula->{context}{parser};
   63   my $I = $parser->{Value}->new($formula,$i);
   64   $r = $parser->{BOP}->new($formula,'*',$r,$I);
   65   $formula->{tree} = $parser->{BOP}->new($formula,'+',$l,$r);
   66 #   return $formula->eval if scalar(%{$formula->{variables}}) == 0;
   67   return $formula;
   68 }
   69 
   70 #
   71 #  Return the complex type
   72 #
   73 sub typeRef {return $Value::Type{complex}}
   74 
   75 ##################################################
   76 
   77 #
   78 #  Return a complex if it already is one, otherwise make it one
   79 #    (Guarantees that we have both parts in an array ref)
   80 #
   81 sub promote {
   82   my $x = shift;
   83   return $x if (ref($x) eq $pkg && scalar(@_) == 0);
   84   return $pkg->new($x,@_);
   85 }
   86 #
   87 #  Get the data from the promoted item
   88 #    (guarantees that we have an array with two elements)
   89 #
   90 sub promoteData {@{(promote(shift))->data}}
   91 
   92 ##################################################
   93 #
   94 #  Binary operations
   95 #
   96 
   97 sub add {
   98   my ($l,$r,$flag) = @_;
   99   if ($l->promotePrecedence($r)) {return $r->add($l,!$flag)}
  100   my ($a,$b) = (@{$l->data});
  101   my ($c,$d) = promoteData($r);
  102   return $pkg->make($a + $c, $b + $d);
  103 }
  104 
  105 sub sub {
  106   my ($l,$r,$flag) = @_;
  107   if ($l->promotePrecedence($r)) {return $r->sub($l,!$flag)}
  108   $r = promote($r);
  109   if ($flag) {my $tmp = $l; $l = $r; $r = $tmp}
  110   my ($a,$b) = (@{$l->data});
  111   my ($c,$d) = (@{$r->data});
  112   return $pkg->make($a - $c, $b - $d);
  113 }
  114 
  115 sub mult {
  116   my ($l,$r,$flag) = @_;
  117   if ($l->promotePrecedence($r)) {return $r->mult($l,!$flag)}
  118   my ($a,$b) = (@{$l->data});
  119   my ($c,$d) = promoteData($r);
  120   return $pkg->make($a*$c - $b*$d, $b*$c + $a*$d);
  121 }
  122 
  123 sub div {
  124   my ($l,$r,$flag) = @_;
  125   if ($l->promotePrecedence($r)) {return $r->div($l,!$flag)}
  126   $r = promote($r);
  127   if ($flag) {my $tmp = $l; $l = $r; $r = $tmp}
  128   my ($a,$b) = (@{$l->data});
  129   my ($c,$d) = (@{$r->data});
  130   my $x = $c*$c + $d*$d;
  131   Value::Error("Division by zero") if $x == 0;
  132   return $pkg->make(($a*$c + $b*$d)/$x,($b*$c - $a*$d)/$x);
  133 }
  134 
  135 sub power {
  136   my ($l,$r,$flag) = @_;
  137   if ($l->promotePrecedence($r)) {return $r->power($l,!$flag)}
  138   $r = promote($r);
  139   if ($flag) {my $tmp = $l; $l = $r; $r = $tmp}
  140   my ($a,$b) = (@{$l->data});
  141   my ($c,$d) = (@{$r->data});
  142   return Value::Real->make(1) if ($a eq '1' && $b == 0) || ($c == 0 && $d == 0);
  143   return Value::Real->make(0) if $c > 0 && ($a == 0 && $b == 0);
  144   return exp($r * log($l))
  145  }
  146 
  147 sub equal {
  148   my ($l,$r,$flag) = @_;
  149   my ($a,$b) = (@{$l->data});
  150   my ($c,$d) = promoteData($r);
  151   return $a == $c && $b == $d;
  152 }
  153 
  154 sub compare {
  155   my ($l,$r,$flag) = @_;
  156   if ($l->promotePrecedence($r)) {return $r->power($l,!$flag)}
  157   $r = promote($r);
  158   if ($flag) {my $tmp = $l; $l = $r; $r = $tmp}
  159   my ($a,$b) = (@{$l->data});
  160   my ($c,$d) = (@{$r->data});
  161   return ($a <=> $c) if $a != $c;
  162   return ($b <=> $d);
  163 }
  164 
  165 ##################################################
  166 #
  167 #   Numeric functions
  168 #
  169 
  170 sub arg {
  171   my ($a,$b) = promoteData(@_);
  172   return CORE::atan2($b,$a);
  173 }
  174 
  175 sub mod {
  176   my ($a,$b) = promoteData(@_);
  177   return CORE::sqrt($a*$a+$b*$b);
  178 }
  179 
  180 sub Re {return (promote(@_))->data->[0]}
  181 sub Im {return (promote(@_))->data->[1]}
  182 
  183 sub abs {norm(@_)}
  184 sub norm {
  185   my ($a,$b) = promoteData(@_);
  186   return CORE::sqrt($a*$a+$b*$b);
  187 }
  188 
  189 sub neg {
  190   my ($a,$b) = promoteData(@_);
  191   return $pkg->make(-$a,-$b);
  192 }
  193 
  194 sub conj {
  195   my ($a,$b) = promoteData(@_);
  196   return $pkg->make($a,-$b);
  197 }
  198 
  199 sub exp {
  200   my ($a,$b) = promoteData(@_);
  201   my $e = CORE::exp($a);
  202   my ($c,$s) = (CORE::cos($b),CORE::sin($b));
  203   return $pkg->make($e*$c,$e*$s);
  204 }
  205 
  206 sub log {
  207   my $z = promote(@_);
  208   my ($r,$t) = ($z->mod,$z->arg);
  209   Value::Error("Can't compute log of zero") if ($r == 0);
  210   return $pkg->make(CORE::log($r),$t);
  211 }
  212 
  213 sub sqrt {
  214   my $z = promote(@_);
  215   $z->power(.5);
  216 }
  217 
  218 ##################################################
  219 #
  220 #   Trig functions
  221 #
  222 
  223 # sin(z) = (exp(iz) - exp(-iz))/(2i)
  224 sub sin {
  225   my ($a,$b) = promoteData(@_);
  226   my $e = CORE::exp($b); my $e1 = 1/$e;
  227   $pkg->make(CORE::sin($a)*($e+$e1)/2, CORE::cos($a)*($e-$e1)/2);
  228 }
  229 
  230 # cos(z) = (exp(iz) + exp(-iz))/2
  231 sub cos {
  232   my ($a,$b) = promoteData(@_);
  233   my $e = CORE::exp($b); my $e1 = 1/$e;
  234   $pkg->make(CORE::cos($a)*($e+$e1)/2, CORE::sin($a)*($e1-$e)/2);
  235 }
  236 
  237 # tan(z) = sin(z) / cos(z)
  238 sub tan {CORE::sin($_[0])/CORE::cos($_[0])}
  239 
  240 # sec(z) = 1 / cos(z)
  241 sub sec {1/CORE::cos($_[0])}
  242 
  243 # csc(z) = 1 / sin(z)
  244 sub csc {1/CORE::sin($_[0])}
  245 
  246 # cot(z) = cos(z) / sin(z)
  247 sub cot {CORE::cos($_[0])/CORE::sin($_[0])}
  248 
  249 # acos(z) = -i log(z + sqrt(z^2 - 1))
  250 sub acos {my $z = promote(@_); -$i * CORE::log($z + CORE::sqrt($z*$z - 1))}
  251 
  252 # asin(z) = -i log(iz + sqrt(1 - z^2))
  253 sub asin {my $z = promote(@_); -$i * CORE::log($i*$z + CORE::sqrt(1 - $z*$z))}
  254 
  255 # atan(z) = (i/2) log((i+z)/(i-z))
  256 sub atan {my $z = promote(@_); ($i/2)*CORE::log(($i+$z)/($i-$z))}
  257 
  258 # asec(z) = acos(1/z)
  259 sub asec {acos(1/$_[0])}
  260 
  261 # acsc(z) = asin(1/z)
  262 sub acsc {asin(1/$_[0])}
  263 
  264 # acot(z) = atan(1/z)
  265 sub acot {atan(1/$_[0])}
  266 
  267 # atan2(z1,z2) = atan(z1/z2)
  268 sub atan2 {
  269   my ($l,$r,$flag) = @_;
  270   if ($flag) {my $tmp = $l; $l = $r; $r = $l}
  271   my ($a,$b) = promoteData($l);
  272   my ($c,$d) = promoteData($r);
  273   if ($b == 0) {
  274     return CORE::atan2($a,$c) if $b == 0;
  275     return $pkg->make($pi/2,0) if $a == 0;
  276   }
  277   ($a,$b) = @{atan($l/$r)->data};
  278   $a += $pi if $c <0; $a -= 2*$pi if $a > $pi;
  279   return $pkg->make($a,$b);
  280 }
  281 
  282 ##################################################
  283 #
  284 #   Hyperbolic functions
  285 #
  286 
  287 # sinh(z) = (exp(z) - exp(-z))/2
  288 sub sinh {
  289   my ($a,$b) = promoteData(@_);
  290   my $e = CORE::exp($a); my $e1 = 1/$e;
  291   $pkg->make(CORE::cos($b)*($e-$e1)/2, CORE::sin($b)*($e+$e1)/2);
  292 }
  293 
  294 # cosh(z) = (exp(z) + exp(-z))/2
  295 sub cosh {
  296   my ($a,$b) = promoteData(@_);
  297   my $e = CORE::exp($a); my $e1 = 1/$e;
  298   $pkg->make(CORE::cos($b)*($e+$e1)/2, CORE::sin($b)*($e-$e1)/2);
  299 }
  300 
  301 # tanh(z) = sinh(z) / cosh(z)
  302 sub tanh {sinh($_[0])/cosh($_[0])}
  303 
  304 # sech(z) = 1 / cosh(z)
  305 sub sech {1/cosh($_[0])}
  306 
  307 # csch(z) = 1 / sinh(z)
  308 sub csch {1/sinh($_[0])}
  309 
  310 # coth(z) = cosh(z) / sinh(z)
  311 sub coth {cosh($_[0]) / sinh($_[0])}
  312 
  313 # asinh(z) = log(z + sqrt(z^2 + 1))
  314 sub asinh {my $z = promote(@_); CORE::log($z + CORE::sqrt($z*$z + 1))}
  315 
  316 # acosh(z) = log(z + sqrt(z^2 - 1))
  317 sub acosh {my $z = promote(@_); CORE::log($z + CORE::sqrt($z*$z - 1))}
  318 
  319 # atanh(z) = (1/2) log((1+z) / (1-z))
  320 sub atanh {my $z = promote(@_); CORE::log((1+$z)/(1-$z))/2}
  321 
  322 # asech(z) = acosh(1/z)
  323 sub asech {acosh(1/$_[0])}
  324 
  325 # acsch(z) = asinh(1/z)
  326 sub acsch {asinh(1/$_[0])}
  327 
  328 # acoth(z) = (1/2) log((1+z)/(z-1))
  329 sub acoth {my $z = promote(@_); CORE::log((1+$z)/($z-1))/2}
  330 
  331 ##################################################
  332 
  333 sub string {my $self = shift; Value::Complex::format(@{$self->data},'string',@_)}
  334 sub TeX {my $self = shift; Value::Complex::format(@{$self->data},'TeX',@_)}
  335 
  336 #
  337 #  Try to make a pretty version of the number
  338 #
  339 sub format {
  340   my ($a,$b) = (shift,shift);
  341   my $method = shift || 'string';
  342   my $equation = shift;
  343   $a = Value::Real->make($a) unless ref($a);
  344   $b = Value::Real->make($b) unless ref($b);
  345   my $bi = 'i';
  346   return $a->$method($equation) if $b == 0;
  347   $bi = abs($b)->$method($equation,1) . 'i' if abs($b) ne 1;
  348   $bi = '-' . $bi if $b < 0;
  349   return $bi if $a == 0;
  350   $bi = '+' . $bi if $b > 0;
  351   $a = $a->$method($equation); $a = "($a)" if $a =~ m/E/i;
  352   return $a.$bi;
  353 }
  354 
  355 #
  356 #  Values for i and pi
  357 #
  358 $i = $pkg->make(0,1);
  359 $pi = 4*atan2(1,1);
  360 
  361 #
  362 #  So that we can use 1+3*i rather than 1+3*$i, etc.
  363 #
  364 sub i () {return $i}
  365 sub pi () {return $pi}
  366 
  367 ###########################################################################
  368 
  369 1;
  370 

aubreyja at gmail dot com
ViewVC Help
Powered by ViewVC 1.0.9