[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 5696 - (download) (as text) (annotate)
Sat Jun 14 12:14:21 2008 UTC (11 years, 8 months ago) by dpvc
File size: 10238 byte(s)
Added

	 no strict "refs"

to try to avoid new error checking in Perl 5.10.

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

aubreyja at gmail dot com
ViewVC Help
Powered by ViewVC 1.0.9