[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 3716 - (download) (as text) (annotate)
Sun Oct 16 03:37:17 2005 UTC (14 years, 1 month ago) by dpvc
File size: 9878 byte(s)
In the past, when Value objects were inserted into strings, they would
automatically include parentheses so that if you had $f equal to 1+x
and $g equal to 1-x, then Formula("$f/$g") would mean (1+x)/(1-x)
rather than 1+(x/1)-x, which is what would happen as a straing string
substitution.

The problem is that this would also happen for real numbers, vectors,
and everything else, even when it wasn't necessary.  So if $x=Real(3),
then "Let x = $x" would be "Let x = (3)".

I have changed the behavior of the string concatenation for Value
objects so that parentheses are only added in a few cases: for
Formulas, Complex numbers, and Unions.  This makes the other Value
objects work more like regular variables in strings, but might cause
some problems with strings that are used as formulas.  For example, if
$a = Real(-3), then "x + 2 $a" will become "x + 2 -3", or "x-1" rather
than the expected "x - 6".  (The old approach would have made it "x +
2 (-3)" which would have worked properly).  For the most part, it is
easier to use something like "x + 2*$a" or even "x" + 2*$a in this
case, so the extra trouble of having to avoid parentheses when you
really meant to substitute the value into a string didn't seem worth
it.

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

aubreyja at gmail dot com
ViewVC Help
Powered by ViewVC 1.0.9