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

aubreyja at gmail dot com
ViewVC Help
Powered by ViewVC 1.0.9