View of /trunk/pg/lib/Complex1.pm

    1 #
    2 # Complex numbers and associated mathematical functions
    3 # -- Raphael Manfredi Since Sep 1996
    4 # -- Jarkko Hietaniemi  Since Mar 1997
    5 # -- Daniel S. Lewart Since Sep 1997
    6 #
    7 
    8 require Exporter;
    9 package Complex1;
   10 
   11 use strict;
   12 
   13 use vars qw($VERSION @ISA @EXPORT %EXPORT_TAGS);
   14 
   15 my ( $i, $ip2, %logn );
   16 
   17 $VERSION = sprintf("%s", q$Id$ =~ /(\d+\.\d+)/);
   18 
   19 @ISA = qw(Exporter);
   20 
   21 my @trig = qw(
   22         pi
   23         tan
   24         csc cosec sec cot cotan
   25         asin acos atan
   26         acsc acosec asec acot acotan
   27         sinh cosh tanh
   28         csch cosech sech coth cotanh
   29         asinh acosh atanh
   30         acsch acosech asech acoth acotanh
   31        );
   32 
   33 @EXPORT = (qw(
   34        i Re Im rho theta arg
   35        sqrt log ln
   36        log10 logn cbrt root
   37        cplx cplxe
   38        ),
   39      @trig);
   40 
   41 %EXPORT_TAGS = (
   42     'trig' => [@trig],
   43 );
   44 
   45 use overload
   46   '+' => \&plus,
   47   '-' => \&minus,
   48   '*' => \&multiply,
   49   '/' => \&divide,
   50   '**'  => \&power,
   51   '<=>' => \&spaceship,
   52   'neg' => \&negate,
   53   '~' => \&conjugate,
   54   'abs' => \&abs,
   55   'sqrt'  => \&sqrt,
   56   'exp' => \&exp,
   57   'log' => \&log,
   58   'sin' => \&sin,
   59   'cos' => \&cos,
   60   'tan' => \&tan,
   61   'atan2' => \&atan2,
   62   qw("" stringify);
   63 
   64 #
   65 # Package "privates"
   66 #
   67 
   68 #my $package = 'Math::Complex';   # Package name
   69 my $package = 'Complex1';
   70 my $display = 'cartesian';    # Default display format
   71 my $eps     = 1e-14;      # Epsilon
   72 
   73 #
   74 # Object attributes (internal):
   75 # cartesian [real, imaginary] -- cartesian form
   76 # polar   [rho, theta] -- polar form
   77 # c_dirty   cartesian form not up-to-date
   78 # p_dirty   polar form not up-to-date
   79 # display   display format (package's global when not set)
   80 #
   81 
   82 # Die on bad *make() arguments.
   83 
   84 sub _cannot_make {
   85     die "@{[(caller(1))[3]]}: Cannot take $_[0] of $_[1].\n";
   86 }
   87 
   88 #
   89 # ->make
   90 #
   91 # Create a new complex number (cartesian form)
   92 #
   93 sub make {
   94   my $self = bless {}, shift;
   95   my ($re, $im) = @_;
   96   my $rre = ref $re;
   97   if ( $rre ) {
   98       if ( $rre eq ref $self ) {
   99     $re = Re($re);
  100       } else {
  101     _cannot_make("real part", $rre);
  102       }
  103   }
  104   my $rim = ref $im;
  105   if ( $rim ) {
  106       if ( $rim eq ref $self ) {
  107     $im = Im($im);
  108       } else {
  109     _cannot_make("imaginary part", $rim);
  110       }
  111   }
  112   $self->{'cartesian'} = [ $re, $im ];
  113   $self->{c_dirty} = 0;
  114   $self->{p_dirty} = 1;
  115   $self->display_format('cartesian');
  116   return $self;
  117 }
  118 
  119 #
  120 # ->emake
  121 #
  122 # Create a new complex number (exponential form)
  123 #
  124 sub emake {
  125   my $self = bless {}, shift;
  126   my ($rho, $theta) = @_;
  127   my $rrh = ref $rho;
  128   if ( $rrh ) {
  129       if ( $rrh eq ref $self ) {
  130     $rho = rho($rho);
  131       } else {
  132     _cannot_make("rho", $rrh);
  133       }
  134   }
  135   my $rth = ref $theta;
  136   if ( $rth ) {
  137       if ( $rth eq ref $self ) {
  138     $theta = theta($theta);
  139       } else {
  140     _cannot_make("theta", $rth);
  141       }
  142   }
  143   if ($rho < 0) {
  144       $rho   = -$rho;
  145       $theta = ($theta <= 0) ? $theta + pi() : $theta - pi();
  146   }
  147   $self->{'polar'} = [$rho, $theta];
  148   $self->{p_dirty} = 0;
  149   $self->{c_dirty} = 1;
  150   $self->display_format('polar');
  151   return $self;
  152 }
  153 
  154 sub new { &make }   # For backward compatibility only.
  155 
  156 #
  157 # cplx
  158 #
  159 # Creates a complex number from a (re, im) tuple.
  160 # This avoids the burden of writing Math::Complex->make(re, im).
  161 #
  162 sub cplx {
  163   my ($re, $im) = @_;
  164   return $package->make(defined $re ? $re : 0, defined $im ? $im : 0);
  165 }
  166 # cplx adn cplxe changed by MEG
  167 #
  168 # cplxe
  169 #
  170 # Creates a complex number from a (rho, theta) tuple.
  171 # This avoids the burden of writing Math::Complex->emake(rho, theta).
  172 #
  173 sub cplxe {
  174   my ($rho, $theta) = @_;
  175   return $package->emake(defined $rho ? $rho : 0, defined $theta ? $theta : 0);
  176 }
  177 
  178 #
  179 # pi
  180 #
  181 # The number defined as pi = 180 degrees
  182 #
  183 use constant pi => 4 * CORE::atan2(1, 1);
  184 
  185 #
  186 # pit2
  187 #
  188 # The full circle
  189 #
  190 use constant pit2 => 2 * pi;
  191 
  192 #
  193 # pip2
  194 #
  195 # The quarter circle
  196 #
  197 use constant pip2 => pi / 2;
  198 
  199 #
  200 # deg1
  201 #
  202 # One degree in radians, used in stringify_polar.
  203 #
  204 
  205 use constant deg1 => pi / 180;
  206 
  207 #
  208 # uplog10
  209 #
  210 # Used in log10().
  211 #
  212 use constant uplog10 => 1 / CORE::log(10);
  213 
  214 #
  215 # i
  216 #
  217 # The number defined as i*i = -1;
  218 #
  219 sub i () {
  220         return $i if ($i);
  221   $i = bless {};
  222   $i->{'cartesian'} = [0, 1];
  223   $i->{'polar'}     = [1, pip2];
  224   $i->{c_dirty} = 0;
  225   $i->{p_dirty} = 0;
  226   return $i;
  227 }
  228 
  229 #
  230 # Attribute access/set routines
  231 #
  232 
  233 sub cartesian {$_[0]->{c_dirty} ?
  234        $_[0]->update_cartesian : $_[0]->{'cartesian'}}
  235 sub polar     {$_[0]->{p_dirty} ?
  236        $_[0]->update_polar : $_[0]->{'polar'}}
  237 
  238 sub set_cartesian { $_[0]->{p_dirty}++; $_[0]->{'cartesian'} = $_[1] }
  239 sub set_polar     { $_[0]->{c_dirty}++; $_[0]->{'polar'} = $_[1] }
  240 
  241 #
  242 # ->update_cartesian
  243 #
  244 # Recompute and return the cartesian form, given accurate polar form.
  245 #
  246 sub update_cartesian {
  247   my $self = shift;
  248   my ($r, $t) = @{$self->{'polar'}};
  249   $self->{c_dirty} = 0;
  250   return $self->{'cartesian'} = [$r * CORE::cos($t), $r * CORE::sin($t)];
  251 }
  252 
  253 #
  254 #
  255 # ->update_polar
  256 #
  257 # Recompute and return the polar form, given accurate cartesian form.
  258 #
  259 sub update_polar {
  260   my $self = shift;
  261   my ($x, $y) = @{$self->{'cartesian'}};
  262   $self->{p_dirty} = 0;
  263   return $self->{'polar'} = [0, 0] if $x == 0 && $y == 0;
  264   return $self->{'polar'} = [CORE::sqrt($x*$x + $y*$y), CORE::atan2($y, $x)];
  265 }
  266 
  267 #
  268 # (plus)
  269 #
  270 # Computes z1+z2.
  271 #
  272 sub plus {
  273   my ($z1, $z2, $regular) = @_;
  274   my ($re1, $im1) = @{$z1->cartesian};
  275   $z2 = cplx($z2) unless ref $z2;
  276   my ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
  277   unless (defined $regular) {
  278     $z1->set_cartesian([$re1 + $re2, $im1 + $im2]);
  279     return $z1;
  280   }
  281   return (ref $z1)->make($re1 + $re2, $im1 + $im2);
  282 }
  283 
  284 #
  285 # (minus)
  286 #
  287 # Computes z1-z2.
  288 #
  289 sub minus {
  290   my ($z1, $z2, $inverted) = @_;
  291   my ($re1, $im1) = @{$z1->cartesian};
  292   $z2 = cplx($z2) unless ref $z2;
  293   my ($re2, $im2) = @{$z2->cartesian};
  294   unless (defined $inverted) {
  295     $z1->set_cartesian([$re1 - $re2, $im1 - $im2]);
  296     return $z1;
  297   }
  298 
  299    return ( $inverted) ?
  300     (ref $z1)->make($re2 - $re1, $im2 - $im1) :
  301     (ref $z1)->make($re1 - $re2, $im1 - $im2);
  302 
  303 }
  304 
  305 #
  306 # (multiply)
  307 #
  308 # Computes z1*z2.
  309 #
  310 sub multiply {
  311         my ($z1, $z2, $regular) = @_;
  312   if ($z1->{p_dirty} == 0 and ref $z2 and $z2->{p_dirty} == 0) {
  313       # if both polar better use polar to avoid rounding errors
  314       my ($r1, $t1) = @{$z1->polar};
  315       my ($r2, $t2) = @{$z2->polar};
  316       my $t = $t1 + $t2;
  317       if    ($t >   pi()) { $t -= pit2 }
  318       elsif ($t <= -pi()) { $t += pit2 }
  319       unless (defined $regular) {
  320     $z1->set_polar([$r1 * $r2, $t]);
  321     return $z1;
  322       }
  323       return (ref $z1)->emake($r1 * $r2, $t);
  324   } else {
  325       my ($x1, $y1) = @{$z1->cartesian};
  326       if (ref $z2) {
  327     my ($x2, $y2) = @{$z2->cartesian};
  328     return (ref $z1)->make($x1*$x2-$y1*$y2, $x1*$y2+$y1*$x2);
  329       } else {
  330     return (ref $z1)->make($x1*$z2, $y1*$z2);
  331       }
  332   }
  333 }
  334 
  335 #
  336 # _divbyzero
  337 #
  338 # Die on division by zero.
  339 #
  340 sub _divbyzero {
  341     my $mess = "$_[0]: Division by zero.\n";
  342 
  343     if (defined $_[1]) {
  344   $mess .= "(Because in the definition of $_[0], the divisor ";
  345   $mess .= "$_[1] " unless ($_[1] eq '0');
  346   $mess .= "is 0)\n";
  347     }
  348 
  349     my @up = caller(1);
  350 
  351     $mess .= "Died at $up[1] line $up[2].\n";
  352 
  353     die $mess;
  354 }
  355 
  356 #
  357 # (divide)
  358 #
  359 # Computes z1/z2.
  360 #
  361 sub divide {
  362   my ($z1, $z2, $inverted) = @_;
  363   if ($z1->{p_dirty} == 0 and ref $z2 and $z2->{p_dirty} == 0) {
  364       # if both polar better use polar to avoid rounding errors
  365       my ($r1, $t1) = @{$z1->polar};
  366       my ($r2, $t2) = @{$z2->polar};
  367       my $t;
  368       if ($inverted) {
  369     _divbyzero "$z2/0" if ($r1 == 0);
  370     $t = $t2 - $t1;
  371     if    ($t >   pi()) { $t -= pit2 }
  372     elsif ($t <= -pi()) { $t += pit2 }
  373     return (ref $z1)->emake($r2 / $r1, $t);
  374       } else {
  375     _divbyzero "$z1/0" if ($r2 == 0);
  376     $t = $t1 - $t2;
  377     if    ($t >   pi()) { $t -= pit2 }
  378     elsif ($t <= -pi()) { $t += pit2 }
  379     return (ref $z1)->emake($r1 / $r2, $t);
  380       }
  381   } else {
  382       my ($d, $x2, $y2);
  383       if ($inverted) {
  384     ($x2, $y2) = @{$z1->cartesian};
  385     $d = $x2*$x2 + $y2*$y2;
  386     _divbyzero "$z2/0" if $d == 0;
  387     return (ref $z1)->make(($x2*$z2)/$d, -($y2*$z2)/$d);
  388       } else {
  389     my ($x1, $y1) = @{$z1->cartesian};
  390     if (ref $z2) {
  391         ($x2, $y2) = @{$z2->cartesian};
  392         $d = $x2*$x2 + $y2*$y2;
  393         _divbyzero "$z1/0" if $d == 0;
  394         my $u = ($x1*$x2 + $y1*$y2)/$d;
  395         my $v = ($y1*$x2 - $x1*$y2)/$d;
  396         return (ref $z1)->make($u, $v);
  397     } else {
  398         _divbyzero "$z1/0" if $z2 == 0;
  399         return (ref $z1)->make($x1/$z2, $y1/$z2);
  400     }
  401       }
  402   }
  403 }
  404 
  405 #
  406 # (power)
  407 #
  408 # Computes z1**z2 = exp(z2 * log z1)).
  409 #
  410 sub power {
  411   my ($z1, $z2, $inverted) = @_;
  412   if ($inverted) {
  413       return 1 if $z1 == 0 || $z2 == 1;
  414       return 0 if $z2 == 0 && Re($z1) > 0;
  415   } else {
  416       return 1 if $z2 == 0 || $z1 == 1;
  417       return 0 if $z1 == 0 && Re($z2) > 0;
  418   }
  419   my $w = $inverted ? CORE::exp($z1 * CORE::log($z2))
  420                     : CORE::exp($z2 * CORE::log($z1));
  421   # If both arguments cartesian, return cartesian, else polar.
  422   return $z1->{c_dirty} == 0 &&
  423          (not ref $z2 or $z2->{c_dirty} == 0) ?
  424          cplx(@{$w->cartesian}) : $w;
  425 }
  426 
  427 #
  428 # (spaceship)
  429 #
  430 # Computes z1 <=> z2.
  431 # Sorts on the real part first, then on the imaginary part. Thus 2-4i < 3+8i.
  432 #
  433 sub spaceship {
  434   my ($z1, $z2, $inverted) = @_;
  435   my ($re1, $im1) = ref $z1 ? @{$z1->cartesian} : ($z1, 0);
  436   my ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
  437   my $sgn = $inverted ? -1 : 1;
  438   return $sgn * ($re1 <=> $re2) if $re1 != $re2;
  439   return $sgn * ($im1 <=> $im2);
  440 }
  441 
  442 #
  443 # (negate)
  444 #
  445 # Computes -z.
  446 #
  447 sub negate {
  448   my ($z) = @_;
  449   if ($z->{c_dirty}) {
  450     my ($r, $t) = @{$z->polar};
  451     $t = ($t <= 0) ? $t + pi : $t - pi;
  452     return (ref $z)->emake($r, $t);
  453   }
  454   my ($re, $im) = @{$z->cartesian};
  455   return (ref $z)->make(-$re, -$im);
  456 }
  457 
  458 #
  459 # (conjugate)
  460 #
  461 # Compute complex's conjugate.
  462 #
  463 sub conjugate {
  464   my ($z) = @_;
  465   if ($z->{c_dirty}) {
  466     my ($r, $t) = @{$z->polar};
  467     return (ref $z)->emake($r, -$t);
  468   }
  469   my ($re, $im) = @{$z->cartesian};
  470   return (ref $z)->make($re, -$im);
  471 }
  472 
  473 #
  474 # (abs)
  475 #
  476 # Compute or set complex's norm (rho).
  477 #
  478 sub abs {
  479   my ($z, $rho) = @_;
  480   return $z unless ref $z;
  481   if (defined $rho) {
  482       $z->{'polar'} = [ $rho, ${$z->polar}[1] ];
  483       $z->{p_dirty} = 0;
  484       $z->{c_dirty} = 1;
  485       return $rho;
  486   } else {
  487       return ${$z->polar}[0];
  488   }
  489 }
  490 
  491 sub _theta {
  492     my $theta = $_[0];
  493 
  494     if    ($$theta >   pi()) { $$theta -= pit2 }
  495     elsif ($$theta <= -pi()) { $$theta += pit2 }
  496 }
  497 
  498 #
  499 # arg
  500 #
  501 # Compute or set complex's argument (theta).
  502 #
  503 sub arg {
  504   my ($z, $theta) = @_;
  505   return $z unless ref $z;
  506   if (defined $theta) {
  507       _theta(\$theta);
  508       $z->{'polar'} = [ ${$z->polar}[0], $theta ];
  509       $z->{p_dirty} = 0;
  510       $z->{c_dirty} = 1;
  511   } else {
  512       $theta = ${$z->polar}[1];
  513       _theta(\$theta);
  514   }
  515   return $theta;
  516 }
  517 
  518 #
  519 # (sqrt)
  520 #
  521 # Compute sqrt(z).
  522 #
  523 # It is quite tempting to use wantarray here so that in list context
  524 # sqrt() would return the two solutions.  This, however, would
  525 # break things like
  526 #
  527 # print "sqrt(z) = ", sqrt($z), "\n";
  528 #
  529 # The two values would be printed side by side without no intervening
  530 # whitespace, quite confusing.
  531 # Therefore if you want the two solutions use the root().
  532 #
  533 sub sqrt {
  534   my ($z) = @_;
  535   my ($re, $im) = ref $z ? @{$z->cartesian} : ($z, 0);
  536   return $re < 0 ? cplx(0, CORE::sqrt(-$re)) : CORE::sqrt($re) if $im == 0;
  537   my ($r, $t) = @{$z->polar};
  538   return (ref $z)->emake(CORE::sqrt($r), $t/2);
  539 }
  540 
  541 #
  542 # cbrt
  543 #
  544 # Compute cbrt(z) (cubic root).
  545 #
  546 # Why are we not returning three values?  The same answer as for sqrt().
  547 #
  548 sub cbrt {
  549   my ($z) = @_;
  550   return $z < 0 ? -CORE::exp(CORE::log(-$z)/3) : ($z > 0 ? CORE::exp(CORE::log($z)/3): 0)
  551       unless ref $z;
  552   my ($r, $t) = @{$z->polar};
  553   return (ref $z)->emake(CORE::exp(CORE::log($r)/3), $t/3);
  554 }
  555 
  556 #
  557 # _rootbad
  558 #
  559 # Die on bad root.
  560 #
  561 sub _rootbad {
  562     my $mess = "Root $_[0] not defined, root must be positive integer.\n";
  563 
  564     my @up = caller(1);
  565 
  566     $mess .= "Died at $up[1] line $up[2].\n";
  567 
  568     die $mess;
  569 }
  570 
  571 #
  572 # root
  573 #
  574 # Computes all nth root for z, returning an array whose size is n.
  575 # `n' must be a positive integer.
  576 #
  577 # The roots are given by (for k = 0..n-1):
  578 #
  579 # z^(1/n) = r^(1/n) (cos ((t+2 k pi)/n) + i sin ((t+2 k pi)/n))
  580 #
  581 sub root {
  582   my ($z, $n) = @_;
  583   _rootbad($n) if ($n < 1 or int($n) != $n);
  584   my ($r, $t) = ref $z ? @{$z->polar} : (CORE::abs($z), $z >= 0 ? 0 : pi);
  585   my @root;
  586   my $k;
  587   my $theta_inc = pit2 / $n;
  588   my $rho = $r ** (1/$n);
  589   my $theta;
  590   my $cartesian = ref $z && $z->{c_dirty} == 0;
  591   for ($k = 0, $theta = $t / $n; $k < $n; $k++, $theta += $theta_inc) {
  592       my $w = cplxe($rho, $theta);
  593       # Yes, $cartesian is loop invariant.
  594       push @root, $cartesian ? cplx(@{$w->cartesian}) : $w;
  595   }
  596   return @root;
  597 }
  598 
  599 #
  600 # Re
  601 #
  602 # Return or set Re(z).
  603 #
  604 sub Re {
  605   my ($z, $Re) = @_;
  606   return $z unless ref $z;
  607   if (defined $Re) {
  608       $z->{'cartesian'} = [ $Re, ${$z->cartesian}[1] ];
  609       $z->{c_dirty} = 0;
  610       $z->{p_dirty} = 1;
  611   } else {
  612       return ${$z->cartesian}[0];
  613   }
  614 }
  615 
  616 #
  617 # Im
  618 #
  619 # Return or set Im(z).
  620 #
  621 sub Im {
  622   my ($z, $Im) = @_;
  623   return $z unless ref $z;
  624   if (defined $Im) {
  625       $z->{'cartesian'} = [ ${$z->cartesian}[0], $Im ];
  626       $z->{c_dirty} = 0;
  627       $z->{p_dirty} = 1;
  628   } else {
  629       return ${$z->cartesian}[1];
  630   }
  631 }
  632 
  633 #
  634 # rho
  635 #
  636 # Return or set rho(w).
  637 #
  638 sub rho {
  639     Complex1::abs(@_);
  640 }
  641 
  642 #
  643 # theta
  644 #
  645 # Return or set theta(w).
  646 #
  647 sub theta {
  648     Complex1::arg(@_);
  649 }
  650 
  651 #
  652 # (exp)
  653 #
  654 # Computes exp(z).
  655 #
  656 sub exp {
  657   my ($z) = @_;
  658   my ($x, $y) = @{$z->cartesian};
  659   return (ref $z)->emake(CORE::exp($x), $y);
  660 }
  661 
  662 #
  663 # _logofzero
  664 #
  665 # Die on logarithm of zero.
  666 #
  667 sub _logofzero {
  668     my $mess = "$_[0]: Logarithm of zero.\n";
  669 
  670     if (defined $_[1]) {
  671   $mess .= "(Because in the definition of $_[0], the argument ";
  672   $mess .= "$_[1] " unless ($_[1] eq '0');
  673   $mess .= "is 0)\n";
  674     }
  675 
  676     my @up = caller(1);
  677 
  678     $mess .= "Died at $up[1] line $up[2].\n";
  679 
  680     die $mess;
  681 }
  682 
  683 #
  684 # (log)
  685 #
  686 # Compute log(z).
  687 #
  688 sub log {
  689   my ($z) = @_;
  690   unless (ref $z) {
  691       _logofzero("log") if $z == 0;
  692       return $z > 0 ? CORE::log($z) : cplx(CORE::log(-$z), pi);
  693   }
  694   my ($r, $t) = @{$z->polar};
  695   _logofzero("log") if $r == 0;
  696   if    ($t >   pi()) { $t -= pit2 }
  697   elsif ($t <= -pi()) { $t += pit2 }
  698   return (ref $z)->make(CORE::log($r), $t);
  699 }
  700 
  701 #
  702 # ln
  703 #
  704 # Alias for log().
  705 #
  706 sub ln { Complex1::log(@_) }
  707 
  708 #
  709 # log10
  710 #
  711 # Compute log10(z).
  712 #
  713 
  714 sub log10 {
  715   return Complex1::log($_[0]) * uplog10;
  716 }
  717 
  718 #
  719 # logn
  720 #
  721 # Compute logn(z,n) = log(z) / log(n)
  722 #
  723 sub logn {
  724   my ($z, $n) = @_;
  725   $z = cplx($z, 0) unless ref $z;
  726   my $logn = $logn{$n};
  727   $logn = $logn{$n} = CORE::log($n) unless defined $logn; # Cache log(n)
  728   return CORE::log($z) / $logn;
  729 }
  730 
  731 #
  732 # (cos)
  733 #
  734 # Compute cos(z) = (exp(iz) + exp(-iz))/2.
  735 #
  736 sub cos {
  737   my ($z) = @_;
  738   my ($x, $y) = @{$z->cartesian};
  739   my $ey = CORE::exp($y);
  740   my $ey_1 = 1 / $ey;
  741   return (ref $z)->make(CORE::cos($x) * ($ey + $ey_1)/2,
  742             CORE::sin($x) * ($ey_1 - $ey)/2);
  743 }
  744 
  745 #
  746 # (sin)
  747 #
  748 # Compute sin(z) = (exp(iz) - exp(-iz))/2.
  749 #
  750 sub sin {
  751   my ($z) = @_;
  752   my ($x, $y) = @{$z->cartesian};
  753   my $ey = CORE::exp($y);
  754   my $ey_1 = 1 / $ey;
  755   return (ref $z)->make(CORE::sin($x) * ($ey + $ey_1)/2,
  756             CORE::cos($x) * ($ey - $ey_1)/2);
  757 }
  758 
  759 #
  760 # tan
  761 #
  762 # Compute tan(z) = sin(z) / cos(z).
  763 #
  764 sub tan {
  765   my ($z) = @_;
  766   my $cz = CORE::cos($z);
  767   _divbyzero "tan($z)", "cos($z)" if (CORE::abs($cz) < $eps);
  768   return CORE::sin($z) / $cz;
  769 }
  770 
  771 #
  772 # sec
  773 #
  774 # Computes the secant sec(z) = 1 / cos(z).
  775 #
  776 sub sec {
  777   my ($z) = @_;
  778   my $cz = CORE::cos($z);
  779   _divbyzero "sec($z)", "cos($z)" if ($cz == 0);
  780   return 1 / $cz;
  781 }
  782 
  783 #
  784 # csc
  785 #
  786 # Computes the cosecant csc(z) = 1 / sin(z).
  787 #
  788 sub csc {
  789   my ($z) = @_;
  790   my $sz = CORE::sin($z);
  791   _divbyzero "csc($z)", "sin($z)" if ($sz == 0);
  792   return 1 / $sz;
  793 }
  794 
  795 #
  796 # cosec
  797 #
  798 # Alias for csc().
  799 #
  800 sub cosec { Complex1::csc(@_) }
  801 
  802 #
  803 # cot
  804 #
  805 # Computes cot(z) = cos(z) / sin(z).
  806 #
  807 sub cot {
  808   my ($z) = @_;
  809   my $sz = CORE::sin($z);
  810   _divbyzero "cot($z)", "sin($z)" if ($sz == 0);
  811   return CORE::cos($z) / $sz;
  812 }
  813 
  814 #
  815 # cotan
  816 #
  817 # Alias for cot().
  818 #
  819 sub cotan { Complex1::cot(@_) }
  820 
  821 #
  822 # acos
  823 #
  824 # Computes the arc cosine acos(z) = -i log(z + sqrt(z*z-1)).
  825 #
  826 sub acos {
  827   my $z = $_[0];
  828   return CORE::atan2(CORE::sqrt(1-$z*$z), $z) if (! ref $z) && CORE::abs($z) <= 1;
  829   my ($x, $y) = ref $z ? @{$z->cartesian} : ($z, 0);
  830   my $t1 = CORE::sqrt(($x+1)*($x+1) + $y*$y);
  831   my $t2 = CORE::sqrt(($x-1)*($x-1) + $y*$y);
  832   my $alpha = ($t1 + $t2)/2;
  833   my $beta  = ($t1 - $t2)/2;
  834   $alpha = 1 if $alpha < 1;
  835   if    ($beta >  1) { $beta =  1 }
  836   elsif ($beta < -1) { $beta = -1 }
  837   my $u = CORE::atan2(CORE::sqrt(1-$beta*$beta), $beta);
  838   my $v = CORE::log($alpha + CORE::sqrt($alpha*$alpha-1));
  839   $v = -$v if $y > 0 || ($y == 0 && $x < -1);
  840   return $package->make($u, $v);
  841 }
  842 
  843 #
  844 # asin
  845 #
  846 # Computes the arc sine asin(z) = -i log(iz + sqrt(1-z*z)).
  847 #
  848 sub asin {
  849   my $z = $_[0];
  850   return CORE::atan2($z, CORE::sqrt(1-$z*$z)) if (! ref $z) && CORE::abs($z) <= 1;
  851   my ($x, $y) = ref $z ? @{$z->cartesian} : ($z, 0);
  852   my $t1 = CORE::sqrt(($x+1)*($x+1) + $y*$y);
  853   my $t2 = CORE::sqrt(($x-1)*($x-1) + $y*$y);
  854   my $alpha = ($t1 + $t2)/2;
  855   my $beta  = ($t1 - $t2)/2;
  856   $alpha = 1 if $alpha < 1;
  857   if    ($beta >  1) { $beta =  1 }
  858   elsif ($beta < -1) { $beta = -1 }
  859   my $u =  CORE::atan2($beta, CORE::sqrt(1-$beta*$beta));
  860   my $v = -CORE::log($alpha + CORE::sqrt($alpha*$alpha-1));
  861   $v = -$v if $y > 0 || ($y == 0 && $x < -1);
  862   return $package->make($u, $v);
  863 }
  864 
  865 #
  866 # atan
  867 #
  868 # Computes the arc tangent atan(z) = i/2 log((i+z) / (i-z)).
  869 #
  870 sub atan {
  871   my ($z) = @_;
  872   return CORE::atan2($z, 1) unless ref $z;
  873   _divbyzero "atan(i)"  if ( $z == i);
  874   _divbyzero "atan(-i)" if (-$z == i);
  875   my $log = CORE::log((i + $z) / (i - $z));
  876   $ip2 = 0.5 * i unless defined $ip2;
  877   return $ip2 * $log;
  878 }
  879 
  880 #
  881 # asec
  882 #
  883 # Computes the arc secant asec(z) = acos(1 / z).
  884 #
  885 sub asec {
  886   my ($z) = @_;
  887   _divbyzero "asec($z)", $z if ($z == 0);
  888   return acos(1 / $z);
  889 }
  890 
  891 #
  892 # acsc
  893 #
  894 # Computes the arc cosecant acsc(z) = asin(1 / z).
  895 #
  896 sub acsc {
  897   my ($z) = @_;
  898   _divbyzero "acsc($z)", $z if ($z == 0);
  899   return asin(1 / $z);
  900 }
  901 
  902 #
  903 # acosec
  904 #
  905 # Alias for acsc().
  906 #
  907 sub acosec { Complex1::acsc(@_) }
  908 
  909 #
  910 # acot
  911 #
  912 # Computes the arc cotangent acot(z) = atan(1 / z)
  913 #
  914 sub acot {
  915   my ($z) = @_;
  916   _divbyzero "acot(0)"  if (CORE::abs($z)     < $eps);
  917   return ($z >= 0) ? CORE::atan2(1, $z) : CORE::atan2(-1, -$z) unless ref $z;
  918   _divbyzero "acot(i)"  if (CORE::abs($z - i) < $eps);
  919   _logofzero "acot(-i)" if (CORE::abs($z + i) < $eps);
  920   return atan(1 / $z);
  921 }
  922 
  923 #
  924 # acotan
  925 #
  926 # Alias for acot().
  927 #
  928 sub acotan { Complex1::acot(@_) }
  929 
  930 #
  931 # cosh
  932 #
  933 # Computes the hyperbolic cosine cosh(z) = (exp(z) + exp(-z))/2.
  934 #
  935 sub cosh {
  936   my ($z) = @_;
  937   my $ex;
  938   unless (ref $z) {
  939       $ex = CORE::exp($z);
  940       return ($ex + 1/$ex)/2;
  941   }
  942   my ($x, $y) = @{$z->cartesian};
  943   $ex = CORE::exp($x);
  944   my $ex_1 = 1 / $ex;
  945   return (ref $z)->make(CORE::cos($y) * ($ex + $ex_1)/2,
  946             CORE::sin($y) * ($ex - $ex_1)/2);
  947 }
  948 
  949 #
  950 # sinh
  951 #
  952 # Computes the hyperbolic sine sinh(z) = (exp(z) - exp(-z))/2.
  953 #
  954 sub sinh {
  955   my ($z) = @_;
  956   my $ex;
  957   unless (ref $z) {
  958       $ex = CORE::exp($z);
  959       return ($ex - 1/$ex)/2;
  960   }
  961   my ($x, $y) = @{$z->cartesian};
  962   $ex = CORE::exp($x);
  963   my $ex_1 = 1 / $ex;
  964   return (ref $z)->make(CORE::cos($y) * ($ex - $ex_1)/2,
  965             CORE::sin($y) * ($ex + $ex_1)/2);
  966 }
  967 
  968 #
  969 # tanh
  970 #
  971 # Computes the hyperbolic tangent tanh(z) = sinh(z) / cosh(z).
  972 #
  973 sub tanh {
  974   my ($z) = @_;
  975   my $cz = cosh($z);
  976   _divbyzero "tanh($z)", "cosh($z)" if ($cz == 0);
  977   return sinh($z) / $cz;
  978 }
  979 
  980 #
  981 # sech
  982 #
  983 # Computes the hyperbolic secant sech(z) = 1 / cosh(z).
  984 #
  985 sub sech {
  986   my ($z) = @_;
  987   my $cz = cosh($z);
  988   _divbyzero "sech($z)", "cosh($z)" if ($cz == 0);
  989   return 1 / $cz;
  990 }
  991 
  992 #
  993 # csch
  994 #
  995 # Computes the hyperbolic cosecant csch(z) = 1 / sinh(z).
  996 #
  997 sub csch {
  998   my ($z) = @_;
  999   my $sz = sinh($z);
 1000   _divbyzero "csch($z)", "sinh($z)" if ($sz == 0);
 1001   return 1 / $sz;
 1002 }
 1003 
 1004 #
 1005 # cosech
 1006 #
 1007 # Alias for csch().
 1008 #
 1009 sub cosech { Complex1::csch(@_) }
 1010 
 1011 #
 1012 # coth
 1013 #
 1014 # Computes the hyperbolic cotangent coth(z) = cosh(z) / sinh(z).
 1015 #
 1016 sub coth {
 1017   my ($z) = @_;
 1018   my $sz = sinh($z);
 1019   _divbyzero "coth($z)", "sinh($z)" if ($sz == 0);
 1020   return cosh($z) / $sz;
 1021 }
 1022 
 1023 #
 1024 # cotanh
 1025 #
 1026 # Alias for coth().
 1027 #
 1028 sub cotanh { Complex1::coth(@_) }
 1029 
 1030 #
 1031 # acosh
 1032 #
 1033 # Computes the arc hyperbolic cosine acosh(z) = log(z + sqrt(z*z-1)).
 1034 #
 1035 sub acosh {
 1036   my ($z) = @_;
 1037   unless (ref $z) {
 1038       return CORE::log($z + CORE::sqrt($z*$z-1)) if $z >= 1;
 1039       $z = cplx($z, 0);
 1040   }
 1041   my ($re, $im) = @{$z->cartesian};
 1042   if ($im == 0) {
 1043       return cplx(CORE::log($re + CORE::sqrt($re*$re - 1)), 0) if $re >= 1;
 1044       return cplx(0, CORE::atan2(CORE::sqrt(1-$re*$re), $re)) if CORE::abs($re) <= 1;
 1045   }
 1046   return CORE::log($z + CORE::sqrt($z*$z - 1));
 1047 }
 1048 
 1049 #
 1050 # asinh
 1051 #
 1052 # Computes the arc hyperbolic sine asinh(z) = log(z + sqrt(z*z-1))
 1053 #
 1054 sub asinh {
 1055   my ($z) = @_;
 1056   return CORE::log($z + CORE::sqrt($z*$z + 1));
 1057 }
 1058 
 1059 #
 1060 # atanh
 1061 #
 1062 # Computes the arc hyperbolic tangent atanh(z) = 1/2 log((1+z) / (1-z)).
 1063 #
 1064 sub atanh {
 1065   my ($z) = @_;
 1066   unless (ref $z) {
 1067       return CORE::log((1 + $z)/(1 - $z))/2 if CORE::abs($z) < 1;
 1068       $z = cplx($z, 0);
 1069   }
 1070   _divbyzero 'atanh(1)',  "1 - $z" if ($z ==  1);
 1071   _logofzero 'atanh(-1)'           if ($z == -1);
 1072   return 0.5 * CORE::log((1 + $z) / (1 - $z));
 1073 }
 1074 
 1075 #
 1076 # asech
 1077 #
 1078 # Computes the hyperbolic arc secant asech(z) = acosh(1 / z).
 1079 #
 1080 sub asech {
 1081   my ($z) = @_;
 1082   _divbyzero 'asech(0)', $z if ($z == 0);
 1083   return acosh(1 / $z);
 1084 }
 1085 
 1086 #
 1087 # acsch
 1088 #
 1089 # Computes the hyperbolic arc cosecant acsch(z) = asinh(1 / z).
 1090 #
 1091 sub acsch {
 1092   my ($z) = @_;
 1093   _divbyzero 'acsch(0)', $z if ($z == 0);
 1094   return asinh(1 / $z);
 1095 }
 1096 
 1097 #
 1098 # acosech
 1099 #
 1100 # Alias for acosh().
 1101 #
 1102 sub acosech { Complex1::acsch(@_) }
 1103 
 1104 #
 1105 # acoth
 1106 #
 1107 # Computes the arc hyperbolic cotangent acoth(z) = 1/2 log((1+z) / (z-1)).
 1108 #
 1109 sub acoth {
 1110   my ($z) = @_;
 1111   _divbyzero 'acoth(0)'            if (CORE::abs($z)     < $eps);
 1112   unless (ref $z) {
 1113       return CORE::log(($z + 1)/($z - 1))/2 if CORE::abs($z) > 1;
 1114       $z = cplx($z, 0);
 1115   }
 1116   _divbyzero 'acoth(1)',  "$z - 1" if (CORE::abs($z - 1) < $eps);
 1117   _logofzero 'acoth(-1)', "1 / $z" if (CORE::abs($z + 1) < $eps);
 1118   return CORE::log((1 + $z) / ($z - 1)) / 2;
 1119 }
 1120 
 1121 #
 1122 # acotanh
 1123 #
 1124 # Alias for acot().
 1125 #
 1126 sub acotanh { Complex1::acoth(@_) }
 1127 
 1128 #
 1129 # (atan2)
 1130 #
 1131 # Compute atan(z1/z2).
 1132 #
 1133 sub atan2 {
 1134   my ($z1, $z2, $inverted) = @_;
 1135   my ($re1, $im1, $re2, $im2);
 1136   if ($inverted) {
 1137       ($re1, $im1) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
 1138       ($re2, $im2) = @{$z1->cartesian};
 1139   } else {
 1140       ($re1, $im1) = @{$z1->cartesian};
 1141       ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
 1142   }
 1143   if ($im2 == 0) {
 1144       return cplx(CORE::atan2($re1, $re2), 0) if $im1 == 0;
 1145       return cplx(($im1<=>0) * pip2, 0) if $re2 == 0;
 1146   }
 1147   my $w = atan($z1/$z2);
 1148   my ($u, $v) = ref $w ? @{$w->cartesian} : ($w, 0);
 1149   $u += pi   if $re2 < 0;
 1150   $u -= pit2 if $u > pi;
 1151   return cplx($u, $v);
 1152 }
 1153 
 1154 #
 1155 # display_format
 1156 # ->display_format
 1157 #
 1158 # Set (fetch if no argument) display format for all complex numbers that
 1159 # don't happen to have overridden it via ->display_format
 1160 #
 1161 # When called as a method, this actually sets the display format for
 1162 # the current object.
 1163 #
 1164 # Valid object formats are 'c' and 'p' for cartesian and polar. The first
 1165 # letter is used actually, so the type can be fully spelled out for clarity.
 1166 #
 1167 sub display_format {
 1168   my $self = shift;
 1169   my $format = undef;
 1170 
 1171   if (ref $self) {      # Called as a method
 1172     $format = shift;
 1173   } else {        # Regular procedure call
 1174     $format = $self;
 1175     undef $self;
 1176   }
 1177 
 1178   if (defined $self) {
 1179     return defined $self->{display} ? $self->{display} : $display
 1180       unless defined $format;
 1181     return $self->{display} = $format;
 1182   }
 1183 
 1184   return $display unless defined $format;
 1185   return $display = $format;
 1186 }
 1187 
 1188 #
 1189 # (stringify)
 1190 #
 1191 # Show nicely formatted complex number under its cartesian or polar form,
 1192 # depending on the current display format:
 1193 #
 1194 # . If a specific display format has been recorded for this object, use it.
 1195 # . Otherwise, use the generic current default for all complex numbers,
 1196 #   which is a package global variable.
 1197 #
 1198 sub stringify {
 1199   my ($z) = shift;
 1200   my $format;
 1201 
 1202   $format = $display;
 1203   $format = $z->{display} if defined $z->{display};
 1204 
 1205   return $z->stringify_polar if $format =~ /^p/i;
 1206   return $z->stringify_cartesian;
 1207 }
 1208 
 1209 #
 1210 # ->stringify_cartesian
 1211 #
 1212 # Stringify as a cartesian representation 'a+bi'.
 1213 #
 1214 sub stringify_cartesian {
 1215   my $z  = shift;
 1216   my ($x, $y) = @{$z->cartesian};
 1217   my ($re, $im);
 1218 
 1219   $x = int($x + ($x < 0 ? -1 : 1) * $eps)
 1220     if int(CORE::abs($x)) != int(CORE::abs($x) + $eps);
 1221   $y = int($y + ($y < 0 ? -1 : 1) * $eps)
 1222     if int(CORE::abs($y)) != int(CORE::abs($y) + $eps);
 1223 
 1224   $re = "$x" if CORE::abs($x) >= $eps;
 1225         if ($y == 1)                           { $im = 'i' }
 1226         elsif ($y == -1)                       { $im = '-i' }
 1227         elsif (CORE::abs($y) >= $eps)                { $im = $y . "i" }
 1228 
 1229   my $str = '';
 1230   $str = $re if defined $re;
 1231   $str .= "+$im" if defined $im;
 1232   $str =~ s/\+-/-/;
 1233   $str =~ s/^\+//;
 1234   $str =~ s/([-+])1i/$1i/; # Not redundant with the above 1/-1 tests.
 1235   $str = '0' unless $str;
 1236 
 1237   return $str;
 1238 }
 1239 
 1240 
 1241 # Helper for stringify_polar, a Greatest Common Divisor with a memory.
 1242 
 1243 sub _gcd {
 1244     my ($a, $b) = @_;
 1245 
 1246     use integer;
 1247 
 1248     # Loops forever if given negative inputs.
 1249 
 1250     if    ($b and $a > $b) { return gcd($a % $b, $b) }
 1251     elsif ($a and $b > $a) { return gcd($b % $a, $a) }
 1252     else                   { return $a ? $a : $b     }
 1253 }
 1254 
 1255 my %gcd;
 1256 
 1257 sub gcd {
 1258     my ($a, $b) = @_;
 1259 
 1260     my $id = "$a $b";
 1261 
 1262     unless (exists $gcd{$id}) {
 1263   $gcd{$id} = _gcd($a, $b);
 1264   $gcd{"$b $a"} = $gcd{$id};
 1265     }
 1266 
 1267     return $gcd{$id};
 1268 }
 1269 
 1270 #
 1271 # ->stringify_polar
 1272 #
 1273 # Stringify as a polar representation '[r,t]'.
 1274 #
 1275 sub stringify_polar {
 1276   my $z  = shift;
 1277   my ($r, $t) = @{$z->polar};
 1278   my $theta;
 1279 
 1280   return '[0,0]' if $r <= $eps;
 1281 
 1282   my $nt = $t / pit2;
 1283   $nt = ($nt - int($nt)) * pit2;
 1284   $nt += pit2 if $nt < 0;     # Range [0, 2pi]
 1285 
 1286   if (CORE::abs($nt) <= $eps)   { $theta = 0 }
 1287   elsif (CORE::abs(pi-$nt) <= $eps) { $theta = 'pi' }
 1288 
 1289   if (defined $theta) {
 1290     $r = int($r + ($r < 0 ? -1 : 1) * $eps)
 1291       if int(CORE::abs($r)) != int(CORE::abs($r) + $eps);
 1292     $theta = int($theta + ($theta < 0 ? -1 : 1) * $eps)
 1293       if ($theta ne 'pi' and
 1294           int(CORE::abs($theta)) != int(CORE::abs($theta) + $eps));
 1295     return "\[$r,$theta\]";
 1296   }
 1297 
 1298   #
 1299   # Okay, number is not a real. Try to identify pi/n and friends...
 1300   #
 1301 
 1302   $nt -= pit2 if $nt > pi;
 1303 
 1304   if (CORE::abs($nt) >= deg1) {
 1305       my ($n, $k, $kpi);
 1306 
 1307       for ($k = 1, $kpi = pi; $k < 10; $k++, $kpi += pi) {
 1308     $n = int($kpi / $nt + ($nt > 0 ? 1 : -1) * 0.5);
 1309     if (CORE::abs($kpi/$n - $nt) <= $eps) {
 1310         $n = CORE::abs($n);
 1311         my $gcd = gcd($k, $n);
 1312         if ($gcd > 1) {
 1313       $k /= $gcd;
 1314       $n /= $gcd;
 1315         }
 1316         next if $n > 360;
 1317         $theta = ($nt < 0 ? '-':'').
 1318            ($k == 1 ? 'pi':"${k}pi");
 1319         $theta .= '/'.$n if $n > 1;
 1320         last;
 1321     }
 1322       }
 1323   }
 1324 
 1325   $theta = $nt unless defined $theta;
 1326 
 1327   $r = int($r + ($r < 0 ? -1 : 1) * $eps)
 1328     if int(CORE::abs($r)) != int(CORE::abs($r) + $eps);
 1329   $theta = int($theta + ($theta < 0 ? -1 : 1) * $eps)
 1330     if ($theta !~ m(^-?\d*pi/\d+$) and
 1331         int(CORE::abs($theta)) != int(CORE::abs($theta) + $eps));
 1332 
 1333   return "\[$r,$theta\]";
 1334 }
 1335 
 1336 1;
 1337 __END__
 1338 
 1339 =head1 NAME
 1340 
 1341 Math::Complex - complex numbers and associated mathematical functions
 1342 
 1343 =head1 SYNOPSIS
 1344 
 1345   use Math::Complex;
 1346 
 1347   $z = Math::Complex->make(5, 6);
 1348   $t = 4 - 3*i + $z;
 1349   $j = cplxe(1, 2*pi/3);
 1350 
 1351 =head1 DESCRIPTION
 1352 
 1353 This package lets you create and manipulate complex numbers. By default,
 1354 I<Perl> limits itself to real numbers, but an extra C<use> statement brings
 1355 full complex support, along with a full set of mathematical functions
 1356 typically associated with and/or extended to complex numbers.
 1357 
 1358 If you wonder what complex numbers are, they were invented to be able to solve
 1359 the following equation:
 1360 
 1361   x*x = -1
 1362 
 1363 and by definition, the solution is noted I<i> (engineers use I<j> instead since
 1364 I<i> usually denotes an intensity, but the name does not matter). The number
 1365 I<i> is a pure I<imaginary> number.
 1366 
 1367 The arithmetics with pure imaginary numbers works just like you would expect
 1368 it with real numbers... you just have to remember that
 1369 
 1370   i*i = -1
 1371 
 1372 so you have:
 1373 
 1374   5i + 7i = i * (5 + 7) = 12i
 1375   4i - 3i = i * (4 - 3) = i
 1376   4i * 2i = -8
 1377   6i / 2i = 3
 1378   1 / i = -i
 1379 
 1380 Complex numbers are numbers that have both a real part and an imaginary
 1381 part, and are usually noted:
 1382 
 1383   a + bi
 1384 
 1385 where C<a> is the I<real> part and C<b> is the I<imaginary> part. The
 1386 arithmetic with complex numbers is straightforward. You have to
 1387 keep track of the real and the imaginary parts, but otherwise the
 1388 rules used for real numbers just apply:
 1389 
 1390   (4 + 3i) + (5 - 2i) = (4 + 5) + i(3 - 2) = 9 + i
 1391   (2 + i) * (4 - i) = 2*4 + 4i -2i -i*i = 8 + 2i + 1 = 9 + 2i
 1392 
 1393 A graphical representation of complex numbers is possible in a plane
 1394 (also called the I<complex plane>, but it's really a 2D plane).
 1395 The number
 1396 
 1397   z = a + bi
 1398 
 1399 is the point whose coordinates are (a, b). Actually, it would
 1400 be the vector originating from (0, 0) to (a, b). It follows that the addition
 1401 of two complex numbers is a vectorial addition.
 1402 
 1403 Since there is a bijection between a point in the 2D plane and a complex
 1404 number (i.e. the mapping is unique and reciprocal), a complex number
 1405 can also be uniquely identified with polar coordinates:
 1406 
 1407   [rho, theta]
 1408 
 1409 where C<rho> is the distance to the origin, and C<theta> the angle between
 1410 the vector and the I<x> axis. There is a notation for this using the
 1411 exponential form, which is:
 1412 
 1413   rho * exp(i * theta)
 1414 
 1415 where I<i> is the famous imaginary number introduced above. Conversion
 1416 between this form and the cartesian form C<a + bi> is immediate:
 1417 
 1418   a = rho * cos(theta)
 1419   b = rho * sin(theta)
 1420 
 1421 which is also expressed by this formula:
 1422 
 1423   z = rho * exp(i * theta) = rho * (cos theta + i * sin theta)
 1424 
 1425 In other words, it's the projection of the vector onto the I<x> and I<y>
 1426 axes. Mathematicians call I<rho> the I<norm> or I<modulus> and I<theta>
 1427 the I<argument> of the complex number. The I<norm> of C<z> will be
 1428 noted C<abs(z)>.
 1429 
 1430 The polar notation (also known as the trigonometric
 1431 representation) is much more handy for performing multiplications and
 1432 divisions of complex numbers, whilst the cartesian notation is better
 1433 suited for additions and subtractions. Real numbers are on the I<x>
 1434 axis, and therefore I<theta> is zero or I<pi>.
 1435 
 1436 All the common operations that can be performed on a real number have
 1437 been defined to work on complex numbers as well, and are merely
 1438 I<extensions> of the operations defined on real numbers. This means
 1439 they keep their natural meaning when there is no imaginary part, provided
 1440 the number is within their definition set.
 1441 
 1442 For instance, the C<sqrt> routine which computes the square root of
 1443 its argument is only defined for non-negative real numbers and yields a
 1444 non-negative real number (it is an application from B<R+> to B<R+>).
 1445 If we allow it to return a complex number, then it can be extended to
 1446 negative real numbers to become an application from B<R> to B<C> (the
 1447 set of complex numbers):
 1448 
 1449   sqrt(x) = x >= 0 ? sqrt(x) : sqrt(-x)*i
 1450 
 1451 It can also be extended to be an application from B<C> to B<C>,
 1452 whilst its restriction to B<R> behaves as defined above by using
 1453 the following definition:
 1454 
 1455   sqrt(z = [r,t]) = sqrt(r) * exp(i * t/2)
 1456 
 1457 Indeed, a negative real number can be noted C<[x,pi]> (the modulus
 1458 I<x> is always non-negative, so C<[x,pi]> is really C<-x>, a negative
 1459 number) and the above definition states that
 1460 
 1461   sqrt([x,pi]) = sqrt(x) * exp(i*pi/2) = [sqrt(x),pi/2] = sqrt(x)*i
 1462 
 1463 which is exactly what we had defined for negative real numbers above.
 1464 The C<sqrt> returns only one of the solutions: if you want the both,
 1465 use the C<root> function.
 1466 
 1467 All the common mathematical functions defined on real numbers that
 1468 are extended to complex numbers share that same property of working
 1469 I<as usual> when the imaginary part is zero (otherwise, it would not
 1470 be called an extension, would it?).
 1471 
 1472 A I<new> operation possible on a complex number that is
 1473 the identity for real numbers is called the I<conjugate>, and is noted
 1474 with an horizontal bar above the number, or C<~z> here.
 1475 
 1476    z = a + bi
 1477   ~z = a - bi
 1478 
 1479 Simple... Now look:
 1480 
 1481   z * ~z = (a + bi) * (a - bi) = a*a + b*b
 1482 
 1483 We saw that the norm of C<z> was noted C<abs(z)> and was defined as the
 1484 distance to the origin, also known as:
 1485 
 1486   rho = abs(z) = sqrt(a*a + b*b)
 1487 
 1488 so
 1489 
 1490   z * ~z = abs(z) ** 2
 1491 
 1492 If z is a pure real number (i.e. C<b == 0>), then the above yields:
 1493 
 1494   a * a = abs(a) ** 2
 1495 
 1496 which is true (C<abs> has the regular meaning for real number, i.e. stands
 1497 for the absolute value). This example explains why the norm of C<z> is
 1498 noted C<abs(z)>: it extends the C<abs> function to complex numbers, yet
 1499 is the regular C<abs> we know when the complex number actually has no
 1500 imaginary part... This justifies I<a posteriori> our use of the C<abs>
 1501 notation for the norm.
 1502 
 1503 =head1 OPERATIONS
 1504 
 1505 Given the following notations:
 1506 
 1507   z1 = a + bi = r1 * exp(i * t1)
 1508   z2 = c + di = r2 * exp(i * t2)
 1509   z = <any complex or real number>
 1510 
 1511 the following (overloaded) operations are supported on complex numbers:
 1512 
 1513   z1 + z2 = (a + c) + i(b + d)
 1514   z1 - z2 = (a - c) + i(b - d)
 1515   z1 * z2 = (r1 * r2) * exp(i * (t1 + t2))
 1516   z1 / z2 = (r1 / r2) * exp(i * (t1 - t2))
 1517   z1 ** z2 = exp(z2 * log z1)
 1518   ~z = a - bi
 1519   abs(z) = r1 = sqrt(a*a + b*b)
 1520   sqrt(z) = sqrt(r1) * exp(i * t/2)
 1521   exp(z) = exp(a) * exp(i * b)
 1522   log(z) = log(r1) + i*t
 1523   sin(z) = 1/2i (exp(i * z1) - exp(-i * z))
 1524   cos(z) = 1/2 (exp(i * z1) + exp(-i * z))
 1525   atan2(z1, z2) = atan(z1/z2)
 1526 
 1527 The following extra operations are supported on both real and complex
 1528 numbers:
 1529 
 1530   Re(z) = a
 1531   Im(z) = b
 1532   arg(z) = t
 1533   abs(z) = r
 1534 
 1535   cbrt(z) = z ** (1/3)
 1536   log10(z) = log(z) / log(10)
 1537   logn(z, n) = log(z) / log(n)
 1538 
 1539   tan(z) = sin(z) / cos(z)
 1540 
 1541   csc(z) = 1 / sin(z)
 1542   sec(z) = 1 / cos(z)
 1543   cot(z) = 1 / tan(z)
 1544 
 1545   asin(z) = -i * log(i*z + sqrt(1-z*z))
 1546   acos(z) = -i * log(z + i*sqrt(1-z*z))
 1547   atan(z) = i/2 * log((i+z) / (i-z))
 1548 
 1549   acsc(z) = asin(1 / z)
 1550   asec(z) = acos(1 / z)
 1551   acot(z) = atan(1 / z) = -i/2 * log((i+z) / (z-i))
 1552 
 1553   sinh(z) = 1/2 (exp(z) - exp(-z))
 1554   cosh(z) = 1/2 (exp(z) + exp(-z))
 1555   tanh(z) = sinh(z) / cosh(z) = (exp(z) - exp(-z)) / (exp(z) + exp(-z))
 1556 
 1557   csch(z) = 1 / sinh(z)
 1558   sech(z) = 1 / cosh(z)
 1559   coth(z) = 1 / tanh(z)
 1560 
 1561   asinh(z) = log(z + sqrt(z*z+1))
 1562   acosh(z) = log(z + sqrt(z*z-1))
 1563   atanh(z) = 1/2 * log((1+z) / (1-z))
 1564 
 1565   acsch(z) = asinh(1 / z)
 1566   asech(z) = acosh(1 / z)
 1567   acoth(z) = atanh(1 / z) = 1/2 * log((1+z) / (z-1))
 1568 
 1569 I<arg>, I<abs>, I<log>, I<csc>, I<cot>, I<acsc>, I<acot>, I<csch>,
 1570 I<coth>, I<acosech>, I<acotanh>, have aliases I<rho>, I<theta>, I<ln>,
 1571 I<cosec>, I<cotan>, I<acosec>, I<acotan>, I<cosech>, I<cotanh>,
 1572 I<acosech>, I<acotanh>, respectively.  C<Re>, C<Im>, C<arg>, C<abs>,
 1573 C<rho>, and C<theta> can be used also also mutators.  The C<cbrt>
 1574 returns only one of the solutions: if you want all three, use the
 1575 C<root> function.
 1576 
 1577 The I<root> function is available to compute all the I<n>
 1578 roots of some complex, where I<n> is a strictly positive integer.
 1579 There are exactly I<n> such roots, returned as a list. Getting the
 1580 number mathematicians call C<j> such that:
 1581 
 1582   1 + j + j*j = 0;
 1583 
 1584 is a simple matter of writing:
 1585 
 1586   $j = ((root(1, 3))[1];
 1587 
 1588 The I<k>th root for C<z = [r,t]> is given by:
 1589 
 1590   (root(z, n))[k] = r**(1/n) * exp(i * (t + 2*k*pi)/n)
 1591 
 1592 The I<spaceship> comparison operator, E<lt>=E<gt>, is also defined. In
 1593 order to ensure its restriction to real numbers is conform to what you
 1594 would expect, the comparison is run on the real part of the complex
 1595 number first, and imaginary parts are compared only when the real
 1596 parts match.
 1597 
 1598 =head1 CREATION
 1599 
 1600 To create a complex number, use either:
 1601 
 1602   $z = Math::Complex->make(3, 4);
 1603   $z = cplx(3, 4);
 1604 
 1605 if you know the cartesian form of the number, or
 1606 
 1607   $z = 3 + 4*i;
 1608 
 1609 if you like. To create a number using the polar form, use either:
 1610 
 1611   $z = Math::Complex->emake(5, pi/3);
 1612   $x = cplxe(5, pi/3);
 1613 
 1614 instead. The first argument is the modulus, the second is the angle
 1615 (in radians, the full circle is 2*pi).  (Mnemonic: C<e> is used as a
 1616 notation for complex numbers in the polar form).
 1617 
 1618 It is possible to write:
 1619 
 1620   $x = cplxe(-3, pi/4);
 1621 
 1622 but that will be silently converted into C<[3,-3pi/4]>, since the modulus
 1623 must be non-negative (it represents the distance to the origin in the complex
 1624 plane).
 1625 
 1626 It is also possible to have a complex number as either argument of
 1627 either the C<make> or C<emake>: the appropriate component of
 1628 the argument will be used.
 1629 
 1630   $z1 = cplx(-2,  1);
 1631   $z2 = cplx($z1, 4);
 1632 
 1633 =head1 STRINGIFICATION
 1634 
 1635 When printed, a complex number is usually shown under its cartesian
 1636 form I<a+bi>, but there are legitimate cases where the polar format
 1637 I<[r,t]> is more appropriate.
 1638 
 1639 By calling the routine C<Complex1::display_format> and supplying either
 1640 C<"polar"> or C<"cartesian">, you override the default display format,
 1641 which is C<"cartesian">. Not supplying any argument returns the current
 1642 setting.
 1643 
 1644 This default can be overridden on a per-number basis by calling the
 1645 C<display_format> method instead. As before, not supplying any argument
 1646 returns the current display format for this number. Otherwise whatever you
 1647 specify will be the new display format for I<this> particular number.
 1648 
 1649 For instance:
 1650 
 1651   use Math::Complex;
 1652 
 1653   Complex1::display_format('polar');
 1654   $j = ((root(1, 3))[1];
 1655   print "j = $j\n";   # Prints "j = [1,2pi/3]
 1656   $j->display_format('cartesian');
 1657   print "j = $j\n";   # Prints "j = -0.5+0.866025403784439i"
 1658 
 1659 The polar format attempts to emphasize arguments like I<k*pi/n>
 1660 (where I<n> is a positive integer and I<k> an integer within [-9,+9]).
 1661 
 1662 =head1 USAGE
 1663 
 1664 Thanks to overloading, the handling of arithmetics with complex numbers
 1665 is simple and almost transparent.
 1666 
 1667 Here are some examples:
 1668 
 1669   use Math::Complex;
 1670 
 1671   $j = cplxe(1, 2*pi/3);  # $j ** 3 == 1
 1672   print "j = $j, j**3 = ", $j ** 3, "\n";
 1673   print "1 + j + j**2 = ", 1 + $j + $j**2, "\n";
 1674 
 1675   $z = -16 + 0*i;     # Force it to be a complex
 1676   print "sqrt($z) = ", sqrt($z), "\n";
 1677 
 1678   $k = exp(i * 2*pi/3);
 1679   print "$j - $k = ", $j - $k, "\n";
 1680 
 1681   $z->Re(3);      # Re, Im, arg, abs,
 1682   $j->arg(2);     # (the last two aka rho, theta)
 1683           # can be used also as mutators.
 1684 
 1685 =head1 ERRORS DUE TO DIVISION BY ZERO OR LOGARITHM OF ZERO
 1686 
 1687 The division (/) and the following functions
 1688 
 1689   log ln  log10 logn
 1690   tan sec csc cot
 1691   atan  asec  acsc  acot
 1692   tanh  sech  csch  coth
 1693   atanh asech acsch acoth
 1694 
 1695 cannot be computed for all arguments because that would mean dividing
 1696 by zero or taking logarithm of zero. These situations cause fatal
 1697 runtime errors looking like this
 1698 
 1699   cot(0): Division by zero.
 1700   (Because in the definition of cot(0), the divisor sin(0) is 0)
 1701   Died at ...
 1702 
 1703 or
 1704 
 1705   atanh(-1): Logarithm of zero.
 1706   Died at...
 1707 
 1708 For the C<csc>, C<cot>, C<asec>, C<acsc>, C<acot>, C<csch>, C<coth>,
 1709 C<asech>, C<acsch>, the argument cannot be C<0> (zero).  For the the
 1710 logarithmic functions and the C<atanh>, C<acoth>, the argument cannot
 1711 be C<1> (one).  For the C<atanh>, C<acoth>, the argument cannot be
 1712 C<-1> (minus one).  For the C<atan>, C<acot>, the argument cannot be
 1713 C<i> (the imaginary unit).  For the C<atan>, C<acoth>, the argument
 1714 cannot be C<-i> (the negative imaginary unit).  For the C<tan>,
 1715 C<sec>, C<tanh>, the argument cannot be I<pi/2 + k * pi>, where I<k>
 1716 is any integer.
 1717 
 1718 Note that because we are operating on approximations of real numbers,
 1719 these errors can happen when merely `too close' to the singularities
 1720 listed above.  For example C<tan(2*atan2(1,1)+1e-15)> will die of
 1721 division by zero.
 1722 
 1723 =head1 ERRORS DUE TO INDIGESTIBLE ARGUMENTS
 1724 
 1725 The C<make> and C<emake> accept both real and complex arguments.
 1726 When they cannot recognize the arguments they will die with error
 1727 messages like the following
 1728 
 1729     Complex1::make: Cannot take real part of ...
 1730     Complex1::make: Cannot take real part of ...
 1731     Complex1::emake: Cannot take rho of ...
 1732     Complex1::emake: Cannot take theta of ...
 1733 
 1734 =head1 BUGS
 1735 
 1736 Saying C<use Math::Complex;> exports many mathematical routines in the
 1737 caller environment and even overrides some (C<sqrt>, C<log>).
 1738 This is construed as a feature by the Authors, actually... ;-)
 1739 
 1740 All routines expect to be given real or complex numbers. Don't attempt to
 1741 use BigFloat, since Perl has currently no rule to disambiguate a '+'
 1742 operation (for instance) between two overloaded entities.
 1743 
 1744 In Cray UNICOS there is some strange numerical instability that results
 1745 in root(), cos(), sin(), cosh(), sinh(), losing accuracy fast.  Beware.
 1746 The bug may be in UNICOS math libs, in UNICOS C compiler, in Math::Complex.
 1747 Whatever it is, it does not manifest itself anywhere else where Perl runs.
 1748 
 1749 =head1 AUTHORS
 1750 
 1751 Raphael Manfredi <F<Raphael_Manfredi@grenoble.hp.com>> and
 1752 Jarkko Hietaniemi <F<jhi@iki.fi>>.
 1753 
 1754 Extensive patches by Daniel S. Lewart <F<d-lewart@uiuc.edu>>.
 1755 
 1756 =cut
 1757 
 1758 1;
 1759 
 1760 # eof