[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 2625 - (download) (as text) (annotate)
Mon Aug 16 18:35:12 2004 UTC (15 years, 4 months ago) by dpvc
File size: 9363 byte(s)
Added string comparison to all Value object classes (to compare the
string value of an object to another string).

Overloaded perl '.' operator to do dot product when the operands are
formulas returning vectors.  (Part of the auto-generation of
formulas).

A few improvements to real and complex class output results.

Made Union class slightly more robust and removed need for makeUnion
method other than in the Union itself.

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

aubreyja at gmail dot com
ViewVC Help
Powered by ViewVC 1.0.9