"Fossies" - the Fresh Open Source Software Archive

Member "install-tl-20200916/tlpkg/tlperl/lib/Math/Complex.pm" (8 Mar 2018, 50065 Bytes) of package /windows/misc/install-tl.zip:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) Perl source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file.

    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 package Math::Complex;
    9 
   10 { use 5.006; }
   11 use strict;
   12 
   13 our $VERSION = 1.59_01;
   14 
   15 use Config;
   16 
   17 our ($Inf, $ExpInf);
   18 our ($vax_float, $has_inf, $has_nan);
   19 
   20 BEGIN {
   21     $vax_float = (pack("d",1) =~ /^[\x80\x10]\x40/);
   22     $has_inf   = !$vax_float;
   23     $has_nan   = !$vax_float;
   24 
   25     unless ($has_inf) {
   26       # For example in vax, there is no Inf,
   27       # and just mentioning the DBL_MAX (1.70141183460469229e+38)
   28       # causes SIGFPE.
   29 
   30       # These are pretty useless without a real infinity,
   31       # but setting them makes for less warnings about their
   32       # undefined values.
   33       $Inf = "Inf";
   34       $ExpInf = "Inf";
   35       return;
   36     }
   37 
   38     my %DBL_MAX =  # These are IEEE 754 maxima.
   39     (
   40       4  => '1.70141183460469229e+38',
   41       8  => '1.7976931348623157e+308',
   42      # AFAICT the 10, 12, and 16-byte long doubles
   43      # all have the same maximum.
   44      10 => '1.1897314953572317650857593266280070162E+4932',
   45      12 => '1.1897314953572317650857593266280070162E+4932',
   46      16 => '1.1897314953572317650857593266280070162E+4932',
   47     );
   48 
   49     my $nvsize = $Config{nvsize} ||
   50             ($Config{uselongdouble} && $Config{longdblsize}) ||
   51                  $Config{doublesize};
   52     die "Math::Complex: Could not figure out nvsize\n"
   53     unless defined $nvsize;
   54     die "Math::Complex: Cannot not figure out max nv (nvsize = $nvsize)\n"
   55     unless defined $DBL_MAX{$nvsize};
   56     my $DBL_MAX = eval $DBL_MAX{$nvsize};
   57     die "Math::Complex: Could not figure out max nv (nvsize = $nvsize)\n"
   58     unless defined $DBL_MAX;
   59     my $BIGGER_THAN_THIS = 1e30;  # Must find something bigger than this.
   60     if ($^O eq 'unicosmk') {
   61     $Inf = $DBL_MAX;
   62     } else {
   63     local $SIG{FPE} = sub { };
   64         local $!;
   65     # We do want an arithmetic overflow, Inf INF inf Infinity.
   66     for my $t (
   67         'exp(99999)',  # Enough even with 128-bit long doubles.
   68         'inf',
   69         'Inf',
   70         'INF',
   71         'infinity',
   72         'Infinity',
   73         'INFINITY',
   74         '1e99999',
   75         ) {
   76         local $^W = 0;
   77         my $i = eval "$t+1.0";
   78         if (defined $i && $i > $BIGGER_THAN_THIS) {
   79         $Inf = $i;
   80         last;
   81         }
   82           }
   83     $Inf = $DBL_MAX unless defined $Inf;  # Oh well, close enough.
   84     die "Math::Complex: Could not get Infinity"
   85         unless $Inf > $BIGGER_THAN_THIS;
   86     $ExpInf = eval 'exp(99999)';
   87       }
   88     # print "# On this machine, Inf = '$Inf'\n";
   89 }
   90 
   91 use Scalar::Util qw(set_prototype);
   92 
   93 use warnings;
   94 no warnings 'syntax';  # To avoid the (_) warnings.
   95 
   96 BEGIN {
   97     # For certain functions that we override, in 5.10 or better
   98     # we can set a smarter prototype that will handle the lexical $_
   99     # (also a 5.10+ feature).
  100     if ($] >= 5.010000) {
  101         set_prototype \&abs, '_';
  102         set_prototype \&cos, '_';
  103         set_prototype \&exp, '_';
  104         set_prototype \&log, '_';
  105         set_prototype \&sin, '_';
  106         set_prototype \&sqrt, '_';
  107     }
  108 }
  109 
  110 my $i;
  111 my %LOGN;
  112 
  113 # Regular expression for floating point numbers.
  114 # These days we could use Scalar::Util::lln(), I guess.
  115 my $gre = qr'\s*([\+\-]?(?:(?:(?:\d+(?:_\d+)*(?:\.\d*(?:_\d+)*)?|\.\d+(?:_\d+)*)(?:[eE][\+\-]?\d+(?:_\d+)*)?))|inf)'i;
  116 
  117 require Exporter;
  118 
  119 our @ISA = qw(Exporter);
  120 
  121 my @trig = qw(
  122           pi
  123           tan
  124           csc cosec sec cot cotan
  125           asin acos atan
  126           acsc acosec asec acot acotan
  127           sinh cosh tanh
  128           csch cosech sech coth cotanh
  129           asinh acosh atanh
  130           acsch acosech asech acoth acotanh
  131          );
  132 
  133 our @EXPORT = (qw(
  134          i Re Im rho theta arg
  135          sqrt log ln
  136          log10 logn cbrt root
  137          cplx cplxe
  138          atan2
  139          ),
  140        @trig);
  141 
  142 my @pi = qw(pi pi2 pi4 pip2 pip4 Inf);
  143 
  144 our @EXPORT_OK = @pi;
  145 
  146 our %EXPORT_TAGS = (
  147     'trig' => [@trig],
  148     'pi' => [@pi],
  149 );
  150 
  151 use overload
  152     '=' => \&_copy,
  153     '+='    => \&_plus,
  154     '+' => \&_plus,
  155     '-='    => \&_minus,
  156     '-' => \&_minus,
  157     '*='    => \&_multiply,
  158     '*' => \&_multiply,
  159     '/='    => \&_divide,
  160     '/' => \&_divide,
  161     '**='   => \&_power,
  162     '**'    => \&_power,
  163     '=='    => \&_numeq,
  164     '<=>'   => \&_spaceship,
  165     'neg'   => \&_negate,
  166     '~' => \&_conjugate,
  167     'abs'   => \&abs,
  168     'sqrt'  => \&sqrt,
  169     'exp'   => \&exp,
  170     'log'   => \&log,
  171     'sin'   => \&sin,
  172     'cos'   => \&cos,
  173     'atan2' => \&atan2,
  174         '""'    => \&_stringify;
  175 
  176 #
  177 # Package "privates"
  178 #
  179 
  180 my %DISPLAY_FORMAT = ('style' => 'cartesian',
  181               'polar_pretty_print' => 1);
  182 my $eps            = 1e-14;     # Epsilon
  183 
  184 #
  185 # Object attributes (internal):
  186 #   cartesian   [real, imaginary] -- cartesian form
  187 #   polar       [rho, theta] -- polar form
  188 #   c_dirty     cartesian form not up-to-date
  189 #   p_dirty     polar form not up-to-date
  190 #   display     display format (package's global when not set)
  191 #
  192 
  193 # Die on bad *make() arguments.
  194 
  195 sub _cannot_make {
  196     die "@{[(caller(1))[3]]}: Cannot take $_[0] of '$_[1]'.\n";
  197 }
  198 
  199 sub _make {
  200     my $arg = shift;
  201     my ($p, $q);
  202 
  203     if ($arg =~ /^$gre$/) {
  204     ($p, $q) = ($1, 0);
  205     } elsif ($arg =~ /^(?:$gre)?$gre\s*i\s*$/) {
  206     ($p, $q) = ($1 || 0, $2);
  207     } elsif ($arg =~ /^\s*\(\s*$gre\s*(?:,\s*$gre\s*)?\)\s*$/) {
  208     ($p, $q) = ($1, $2 || 0);
  209     }
  210 
  211     if (defined $p) {
  212     $p =~ s/^\+//;
  213     $p =~ s/^(-?)inf$/"${1}9**9**9"/e if $has_inf;
  214     $q =~ s/^\+//;
  215     $q =~ s/^(-?)inf$/"${1}9**9**9"/e if $has_inf;
  216     }
  217 
  218     return ($p, $q);
  219 }
  220 
  221 sub _emake {
  222     my $arg = shift;
  223     my ($p, $q);
  224 
  225     if ($arg =~ /^\s*\[\s*$gre\s*(?:,\s*$gre\s*)?\]\s*$/) {
  226     ($p, $q) = ($1, $2 || 0);
  227     } elsif ($arg =~ m!^\s*\[\s*$gre\s*(?:,\s*([-+]?\d*\s*)?pi(?:/\s*(\d+))?\s*)?\]\s*$!) {
  228     ($p, $q) = ($1, ($2 eq '-' ? -1 : ($2 || 1)) * pi() / ($3 || 1));
  229     } elsif ($arg =~ /^\s*\[\s*$gre\s*\]\s*$/) {
  230     ($p, $q) = ($1, 0);
  231     } elsif ($arg =~ /^\s*$gre\s*$/) {
  232     ($p, $q) = ($1, 0);
  233     }
  234 
  235     if (defined $p) {
  236     $p =~ s/^\+//;
  237     $q =~ s/^\+//;
  238     $p =~ s/^(-?)inf$/"${1}9**9**9"/e if $has_inf;
  239     $q =~ s/^(-?)inf$/"${1}9**9**9"/e if $has_inf;
  240     }
  241 
  242     return ($p, $q);
  243 }
  244 
  245 sub _copy {
  246     my $self = shift;
  247     my $clone = {%$self};
  248     if ($self->{'cartesian'}) {
  249     $clone->{'cartesian'} = [@{$self->{'cartesian'}}];
  250     }
  251     if ($self->{'polar'}) {
  252     $clone->{'polar'} = [@{$self->{'polar'}}];
  253     }
  254     bless $clone,__PACKAGE__;
  255     return $clone;
  256 }
  257 
  258 #
  259 # ->make
  260 #
  261 # Create a new complex number (cartesian form)
  262 #
  263 sub make {
  264     my $self = bless {}, shift;
  265     my ($re, $im);
  266     if (@_ == 0) {
  267     ($re, $im) = (0, 0);
  268     } elsif (@_ == 1) {
  269     return (ref $self)->emake($_[0])
  270         if ($_[0] =~ /^\s*\[/);
  271     ($re, $im) = _make($_[0]);
  272     } elsif (@_ == 2) {
  273     ($re, $im) = @_;
  274     }
  275     if (defined $re) {
  276     _cannot_make("real part",      $re) unless $re =~ /^$gre$/;
  277     }
  278     $im ||= 0;
  279     _cannot_make("imaginary part", $im) unless $im =~ /^$gre$/;
  280     $self->_set_cartesian([$re, $im ]);
  281     $self->display_format('cartesian');
  282 
  283     return $self;
  284 }
  285 
  286 #
  287 # ->emake
  288 #
  289 # Create a new complex number (exponential form)
  290 #
  291 sub emake {
  292     my $self = bless {}, shift;
  293     my ($rho, $theta);
  294     if (@_ == 0) {
  295     ($rho, $theta) = (0, 0);
  296     } elsif (@_ == 1) {
  297     return (ref $self)->make($_[0])
  298         if ($_[0] =~ /^\s*\(/ || $_[0] =~ /i\s*$/);
  299     ($rho, $theta) = _emake($_[0]);
  300     } elsif (@_ == 2) {
  301     ($rho, $theta) = @_;
  302     }
  303     if (defined $rho && defined $theta) {
  304     if ($rho < 0) {
  305         $rho   = -$rho;
  306         $theta = ($theta <= 0) ? $theta + pi() : $theta - pi();
  307     }
  308     }
  309     if (defined $rho) {
  310     _cannot_make("rho",   $rho)   unless $rho   =~ /^$gre$/;
  311     }
  312     $theta ||= 0;
  313     _cannot_make("theta", $theta) unless $theta =~ /^$gre$/;
  314     $self->_set_polar([$rho, $theta]);
  315     $self->display_format('polar');
  316 
  317     return $self;
  318 }
  319 
  320 sub new { &make }       # For backward compatibility only.
  321 
  322 #
  323 # cplx
  324 #
  325 # Creates a complex number from a (re, im) tuple.
  326 # This avoids the burden of writing Math::Complex->make(re, im).
  327 #
  328 sub cplx {
  329     return __PACKAGE__->make(@_);
  330 }
  331 
  332 #
  333 # cplxe
  334 #
  335 # Creates a complex number from a (rho, theta) tuple.
  336 # This avoids the burden of writing Math::Complex->emake(rho, theta).
  337 #
  338 sub cplxe {
  339     return __PACKAGE__->emake(@_);
  340 }
  341 
  342 #
  343 # pi
  344 #
  345 # The number defined as pi = 180 degrees
  346 #
  347 sub pi () { 4 * CORE::atan2(1, 1) }
  348 
  349 #
  350 # pi2
  351 #
  352 # The full circle
  353 #
  354 sub pi2 () { 2 * pi }
  355 
  356 #
  357 # pi4
  358 #
  359 # The full circle twice.
  360 #
  361 sub pi4 () { 4 * pi }
  362 
  363 #
  364 # pip2
  365 #
  366 # The quarter circle
  367 #
  368 sub pip2 () { pi / 2 }
  369 
  370 #
  371 # pip4
  372 #
  373 # The eighth circle.
  374 #
  375 sub pip4 () { pi / 4 }
  376 
  377 #
  378 # _uplog10
  379 #
  380 # Used in log10().
  381 #
  382 sub _uplog10 () { 1 / CORE::log(10) }
  383 
  384 #
  385 # i
  386 #
  387 # The number defined as i*i = -1;
  388 #
  389 sub i () {
  390         return $i if ($i);
  391     $i = bless {};
  392     $i->{'cartesian'} = [0, 1];
  393     $i->{'polar'}     = [1, pip2];
  394     $i->{c_dirty} = 0;
  395     $i->{p_dirty} = 0;
  396     return $i;
  397 }
  398 
  399 #
  400 # _ip2
  401 #
  402 # Half of i.
  403 #
  404 sub _ip2 () { i / 2 }
  405 
  406 #
  407 # Attribute access/set routines
  408 #
  409 
  410 sub _cartesian {$_[0]->{c_dirty} ?
  411            $_[0]->_update_cartesian : $_[0]->{'cartesian'}}
  412 sub _polar     {$_[0]->{p_dirty} ?
  413            $_[0]->_update_polar : $_[0]->{'polar'}}
  414 
  415 sub _set_cartesian { $_[0]->{p_dirty}++; $_[0]->{c_dirty} = 0;
  416              $_[0]->{'cartesian'} = $_[1] }
  417 sub _set_polar     { $_[0]->{c_dirty}++; $_[0]->{p_dirty} = 0;
  418              $_[0]->{'polar'} = $_[1] }
  419 
  420 #
  421 # ->_update_cartesian
  422 #
  423 # Recompute and return the cartesian form, given accurate polar form.
  424 #
  425 sub _update_cartesian {
  426     my $self = shift;
  427     my ($r, $t) = @{$self->{'polar'}};
  428     $self->{c_dirty} = 0;
  429     return $self->{'cartesian'} = [$r * CORE::cos($t), $r * CORE::sin($t)];
  430 }
  431 
  432 #
  433 #
  434 # ->_update_polar
  435 #
  436 # Recompute and return the polar form, given accurate cartesian form.
  437 #
  438 sub _update_polar {
  439     my $self = shift;
  440     my ($x, $y) = @{$self->{'cartesian'}};
  441     $self->{p_dirty} = 0;
  442     return $self->{'polar'} = [0, 0] if $x == 0 && $y == 0;
  443     return $self->{'polar'} = [CORE::sqrt($x*$x + $y*$y),
  444                    CORE::atan2($y, $x)];
  445 }
  446 
  447 #
  448 # (_plus)
  449 #
  450 # Computes z1+z2.
  451 #
  452 sub _plus {
  453     my ($z1, $z2, $regular) = @_;
  454     my ($re1, $im1) = @{$z1->_cartesian};
  455     $z2 = cplx($z2) unless ref $z2;
  456     my ($re2, $im2) = ref $z2 ? @{$z2->_cartesian} : ($z2, 0);
  457     unless (defined $regular) {
  458         $z1->_set_cartesian([$re1 + $re2, $im1 + $im2]);
  459         return $z1;
  460     }
  461     return (ref $z1)->make($re1 + $re2, $im1 + $im2);
  462 }
  463 
  464 #
  465 # (_minus)
  466 #
  467 # Computes z1-z2.
  468 #
  469 sub _minus {
  470     my ($z1, $z2, $inverted) = @_;
  471     my ($re1, $im1) = @{$z1->_cartesian};
  472     $z2 = cplx($z2) unless ref $z2;
  473     my ($re2, $im2) = @{$z2->_cartesian};
  474     unless (defined $inverted) {
  475         $z1->_set_cartesian([$re1 - $re2, $im1 - $im2]);
  476         return $z1;
  477     }
  478     return $inverted ?
  479         (ref $z1)->make($re2 - $re1, $im2 - $im1) :
  480         (ref $z1)->make($re1 - $re2, $im1 - $im2);
  481 
  482 }
  483 
  484 #
  485 # (_multiply)
  486 #
  487 # Computes z1*z2.
  488 #
  489 sub _multiply {
  490         my ($z1, $z2, $regular) = @_;
  491     if ($z1->{p_dirty} == 0 and ref $z2 and $z2->{p_dirty} == 0) {
  492         # if both polar better use polar to avoid rounding errors
  493         my ($r1, $t1) = @{$z1->_polar};
  494         my ($r2, $t2) = @{$z2->_polar};
  495         my $t = $t1 + $t2;
  496         if    ($t >   pi()) { $t -= pi2 }
  497         elsif ($t <= -pi()) { $t += pi2 }
  498         unless (defined $regular) {
  499         $z1->_set_polar([$r1 * $r2, $t]);
  500         return $z1;
  501         }
  502         return (ref $z1)->emake($r1 * $r2, $t);
  503     } else {
  504         my ($x1, $y1) = @{$z1->_cartesian};
  505         if (ref $z2) {
  506         my ($x2, $y2) = @{$z2->_cartesian};
  507         return (ref $z1)->make($x1*$x2-$y1*$y2, $x1*$y2+$y1*$x2);
  508         } else {
  509         return (ref $z1)->make($x1*$z2, $y1*$z2);
  510         }
  511     }
  512 }
  513 
  514 #
  515 # _divbyzero
  516 #
  517 # Die on division by zero.
  518 #
  519 sub _divbyzero {
  520     my $mess = "$_[0]: Division by zero.\n";
  521 
  522     if (defined $_[1]) {
  523     $mess .= "(Because in the definition of $_[0], the divisor ";
  524     $mess .= "$_[1] " unless ("$_[1]" eq '0');
  525     $mess .= "is 0)\n";
  526     }
  527 
  528     my @up = caller(1);
  529 
  530     $mess .= "Died at $up[1] line $up[2].\n";
  531 
  532     die $mess;
  533 }
  534 
  535 #
  536 # (_divide)
  537 #
  538 # Computes z1/z2.
  539 #
  540 sub _divide {
  541     my ($z1, $z2, $inverted) = @_;
  542     if ($z1->{p_dirty} == 0 and ref $z2 and $z2->{p_dirty} == 0) {
  543         # if both polar better use polar to avoid rounding errors
  544         my ($r1, $t1) = @{$z1->_polar};
  545         my ($r2, $t2) = @{$z2->_polar};
  546         my $t;
  547         if ($inverted) {
  548         _divbyzero "$z2/0" if ($r1 == 0);
  549         $t = $t2 - $t1;
  550         if    ($t >   pi()) { $t -= pi2 }
  551         elsif ($t <= -pi()) { $t += pi2 }
  552         return (ref $z1)->emake($r2 / $r1, $t);
  553         } else {
  554         _divbyzero "$z1/0" if ($r2 == 0);
  555         $t = $t1 - $t2;
  556         if    ($t >   pi()) { $t -= pi2 }
  557         elsif ($t <= -pi()) { $t += pi2 }
  558         return (ref $z1)->emake($r1 / $r2, $t);
  559         }
  560     } else {
  561         my ($d, $x2, $y2);
  562         if ($inverted) {
  563         ($x2, $y2) = @{$z1->_cartesian};
  564         $d = $x2*$x2 + $y2*$y2;
  565         _divbyzero "$z2/0" if $d == 0;
  566         return (ref $z1)->make(($x2*$z2)/$d, -($y2*$z2)/$d);
  567         } else {
  568         my ($x1, $y1) = @{$z1->_cartesian};
  569         if (ref $z2) {
  570             ($x2, $y2) = @{$z2->_cartesian};
  571             $d = $x2*$x2 + $y2*$y2;
  572             _divbyzero "$z1/0" if $d == 0;
  573             my $u = ($x1*$x2 + $y1*$y2)/$d;
  574             my $v = ($y1*$x2 - $x1*$y2)/$d;
  575             return (ref $z1)->make($u, $v);
  576         } else {
  577             _divbyzero "$z1/0" if $z2 == 0;
  578             return (ref $z1)->make($x1/$z2, $y1/$z2);
  579         }
  580         }
  581     }
  582 }
  583 
  584 #
  585 # (_power)
  586 #
  587 # Computes z1**z2 = exp(z2 * log z1)).
  588 #
  589 sub _power {
  590     my ($z1, $z2, $inverted) = @_;
  591     if ($inverted) {
  592         return 1 if $z1 == 0 || $z2 == 1;
  593         return 0 if $z2 == 0 && Re($z1) > 0;
  594     } else {
  595         return 1 if $z2 == 0 || $z1 == 1;
  596         return 0 if $z1 == 0 && Re($z2) > 0;
  597     }
  598     my $w = $inverted ? &exp($z1 * &log($z2))
  599                       : &exp($z2 * &log($z1));
  600     # If both arguments cartesian, return cartesian, else polar.
  601     return $z1->{c_dirty} == 0 &&
  602            (not ref $z2 or $z2->{c_dirty} == 0) ?
  603            cplx(@{$w->_cartesian}) : $w;
  604 }
  605 
  606 #
  607 # (_spaceship)
  608 #
  609 # Computes z1 <=> z2.
  610 # Sorts on the real part first, then on the imaginary part. Thus 2-4i < 3+8i.
  611 #
  612 sub _spaceship {
  613     my ($z1, $z2, $inverted) = @_;
  614     my ($re1, $im1) = ref $z1 ? @{$z1->_cartesian} : ($z1, 0);
  615     my ($re2, $im2) = ref $z2 ? @{$z2->_cartesian} : ($z2, 0);
  616     my $sgn = $inverted ? -1 : 1;
  617     return $sgn * ($re1 <=> $re2) if $re1 != $re2;
  618     return $sgn * ($im1 <=> $im2);
  619 }
  620 
  621 #
  622 # (_numeq)
  623 #
  624 # Computes z1 == z2.
  625 #
  626 # (Required in addition to _spaceship() because of NaNs.)
  627 sub _numeq {
  628     my ($z1, $z2, $inverted) = @_;
  629     my ($re1, $im1) = ref $z1 ? @{$z1->_cartesian} : ($z1, 0);
  630     my ($re2, $im2) = ref $z2 ? @{$z2->_cartesian} : ($z2, 0);
  631     return $re1 == $re2 && $im1 == $im2 ? 1 : 0;
  632 }
  633 
  634 #
  635 # (_negate)
  636 #
  637 # Computes -z.
  638 #
  639 sub _negate {
  640     my ($z) = @_;
  641     if ($z->{c_dirty}) {
  642         my ($r, $t) = @{$z->_polar};
  643         $t = ($t <= 0) ? $t + pi : $t - pi;
  644         return (ref $z)->emake($r, $t);
  645     }
  646     my ($re, $im) = @{$z->_cartesian};
  647     return (ref $z)->make(-$re, -$im);
  648 }
  649 
  650 #
  651 # (_conjugate)
  652 #
  653 # Compute complex's _conjugate.
  654 #
  655 sub _conjugate {
  656     my ($z) = @_;
  657     if ($z->{c_dirty}) {
  658         my ($r, $t) = @{$z->_polar};
  659         return (ref $z)->emake($r, -$t);
  660     }
  661     my ($re, $im) = @{$z->_cartesian};
  662     return (ref $z)->make($re, -$im);
  663 }
  664 
  665 #
  666 # (abs)
  667 #
  668 # Compute or set complex's norm (rho).
  669 #
  670 sub abs {
  671     my ($z, $rho) = @_ ? @_ : $_;
  672     unless (ref $z) {
  673         if (@_ == 2) {
  674         $_[0] = $_[1];
  675         } else {
  676         return CORE::abs($z);
  677         }
  678     }
  679     if (defined $rho) {
  680         $z->{'polar'} = [ $rho, ${$z->_polar}[1] ];
  681         $z->{p_dirty} = 0;
  682         $z->{c_dirty} = 1;
  683         return $rho;
  684     } else {
  685         return ${$z->_polar}[0];
  686     }
  687 }
  688 
  689 sub _theta {
  690     my $theta = $_[0];
  691 
  692     if    ($$theta >   pi()) { $$theta -= pi2 }
  693     elsif ($$theta <= -pi()) { $$theta += pi2 }
  694 }
  695 
  696 #
  697 # arg
  698 #
  699 # Compute or set complex's argument (theta).
  700 #
  701 sub arg {
  702     my ($z, $theta) = @_;
  703     return $z unless ref $z;
  704     if (defined $theta) {
  705         _theta(\$theta);
  706         $z->{'polar'} = [ ${$z->_polar}[0], $theta ];
  707         $z->{p_dirty} = 0;
  708         $z->{c_dirty} = 1;
  709     } else {
  710         $theta = ${$z->_polar}[1];
  711         _theta(\$theta);
  712     }
  713     return $theta;
  714 }
  715 
  716 #
  717 # (sqrt)
  718 #
  719 # Compute sqrt(z).
  720 #
  721 # It is quite tempting to use wantarray here so that in list context
  722 # sqrt() would return the two solutions.  This, however, would
  723 # break things like
  724 #
  725 #   print "sqrt(z) = ", sqrt($z), "\n";
  726 #
  727 # The two values would be printed side by side without no intervening
  728 # whitespace, quite confusing.
  729 # Therefore if you want the two solutions use the root().
  730 #
  731 sub sqrt {
  732     my ($z) = @_ ? $_[0] : $_;
  733     my ($re, $im) = ref $z ? @{$z->_cartesian} : ($z, 0);
  734     return $re < 0 ? cplx(0, CORE::sqrt(-$re)) : CORE::sqrt($re)
  735         if $im == 0;
  736     my ($r, $t) = @{$z->_polar};
  737     return (ref $z)->emake(CORE::sqrt($r), $t/2);
  738 }
  739 
  740 #
  741 # cbrt
  742 #
  743 # Compute cbrt(z) (cubic root).
  744 #
  745 # Why are we not returning three values?  The same answer as for sqrt().
  746 #
  747 sub cbrt {
  748     my ($z) = @_;
  749     return $z < 0 ?
  750         -CORE::exp(CORE::log(-$z)/3) :
  751         ($z > 0 ? CORE::exp(CORE::log($z)/3): 0)
  752         unless ref $z;
  753     my ($r, $t) = @{$z->_polar};
  754     return 0 if $r == 0;
  755     return (ref $z)->emake(CORE::exp(CORE::log($r)/3), $t/3);
  756 }
  757 
  758 #
  759 # _rootbad
  760 #
  761 # Die on bad root.
  762 #
  763 sub _rootbad {
  764     my $mess = "Root '$_[0]' illegal, root rank must be positive integer.\n";
  765 
  766     my @up = caller(1);
  767 
  768     $mess .= "Died at $up[1] line $up[2].\n";
  769 
  770     die $mess;
  771 }
  772 
  773 #
  774 # root
  775 #
  776 # Computes all nth root for z, returning an array whose size is n.
  777 # `n' must be a positive integer.
  778 #
  779 # The roots are given by (for k = 0..n-1):
  780 #
  781 # z^(1/n) = r^(1/n) (cos ((t+2 k pi)/n) + i sin ((t+2 k pi)/n))
  782 #
  783 sub root {
  784     my ($z, $n, $k) = @_;
  785     _rootbad($n) if ($n < 1 or int($n) != $n);
  786     my ($r, $t) = ref $z ?
  787         @{$z->_polar} : (CORE::abs($z), $z >= 0 ? 0 : pi);
  788     my $theta_inc = pi2 / $n;
  789     my $rho = $r ** (1/$n);
  790     my $cartesian = ref $z && $z->{c_dirty} == 0;
  791     if (@_ == 2) {
  792         my @root;
  793         for (my $i = 0, my $theta = $t / $n;
  794          $i < $n;
  795          $i++, $theta += $theta_inc) {
  796         my $w = cplxe($rho, $theta);
  797         # Yes, $cartesian is loop invariant.
  798         push @root, $cartesian ? cplx(@{$w->_cartesian}) : $w;
  799         }
  800         return @root;
  801     } elsif (@_ == 3) {
  802         my $w = cplxe($rho, $t / $n + $k * $theta_inc);
  803         return $cartesian ? cplx(@{$w->_cartesian}) : $w;
  804     }
  805 }
  806 
  807 #
  808 # Re
  809 #
  810 # Return or set Re(z).
  811 #
  812 sub Re {
  813     my ($z, $Re) = @_;
  814     return $z unless ref $z;
  815     if (defined $Re) {
  816         $z->{'cartesian'} = [ $Re, ${$z->_cartesian}[1] ];
  817         $z->{c_dirty} = 0;
  818         $z->{p_dirty} = 1;
  819     } else {
  820         return ${$z->_cartesian}[0];
  821     }
  822 }
  823 
  824 #
  825 # Im
  826 #
  827 # Return or set Im(z).
  828 #
  829 sub Im {
  830     my ($z, $Im) = @_;
  831     return 0 unless ref $z;
  832     if (defined $Im) {
  833         $z->{'cartesian'} = [ ${$z->_cartesian}[0], $Im ];
  834         $z->{c_dirty} = 0;
  835         $z->{p_dirty} = 1;
  836     } else {
  837         return ${$z->_cartesian}[1];
  838     }
  839 }
  840 
  841 #
  842 # rho
  843 #
  844 # Return or set rho(w).
  845 #
  846 sub rho {
  847     Math::Complex::abs(@_);
  848 }
  849 
  850 #
  851 # theta
  852 #
  853 # Return or set theta(w).
  854 #
  855 sub theta {
  856     Math::Complex::arg(@_);
  857 }
  858 
  859 #
  860 # (exp)
  861 #
  862 # Computes exp(z).
  863 #
  864 sub exp {
  865     my ($z) = @_ ? @_ : $_;
  866     return CORE::exp($z) unless ref $z;
  867     my ($x, $y) = @{$z->_cartesian};
  868     return (ref $z)->emake(CORE::exp($x), $y);
  869 }
  870 
  871 #
  872 # _logofzero
  873 #
  874 # Die on logarithm of zero.
  875 #
  876 sub _logofzero {
  877     my $mess = "$_[0]: Logarithm of zero.\n";
  878 
  879     if (defined $_[1]) {
  880     $mess .= "(Because in the definition of $_[0], the argument ";
  881     $mess .= "$_[1] " unless ($_[1] eq '0');
  882     $mess .= "is 0)\n";
  883     }
  884 
  885     my @up = caller(1);
  886 
  887     $mess .= "Died at $up[1] line $up[2].\n";
  888 
  889     die $mess;
  890 }
  891 
  892 #
  893 # (log)
  894 #
  895 # Compute log(z).
  896 #
  897 sub log {
  898     my ($z) = @_ ? @_ : $_;
  899     unless (ref $z) {
  900         _logofzero("log") if $z == 0;
  901         return $z > 0 ? CORE::log($z) : cplx(CORE::log(-$z), pi);
  902     }
  903     my ($r, $t) = @{$z->_polar};
  904     _logofzero("log") if $r == 0;
  905     if    ($t >   pi()) { $t -= pi2 }
  906     elsif ($t <= -pi()) { $t += pi2 }
  907     return (ref $z)->make(CORE::log($r), $t);
  908 }
  909 
  910 #
  911 # ln
  912 #
  913 # Alias for log().
  914 #
  915 sub ln { Math::Complex::log(@_) }
  916 
  917 #
  918 # log10
  919 #
  920 # Compute log10(z).
  921 #
  922 
  923 sub log10 {
  924     return Math::Complex::log($_[0]) * _uplog10;
  925 }
  926 
  927 #
  928 # logn
  929 #
  930 # Compute logn(z,n) = log(z) / log(n)
  931 #
  932 sub logn {
  933     my ($z, $n) = @_;
  934     $z = cplx($z, 0) unless ref $z;
  935     my $logn = $LOGN{$n};
  936     $logn = $LOGN{$n} = CORE::log($n) unless defined $logn; # Cache log(n)
  937     return &log($z) / $logn;
  938 }
  939 
  940 #
  941 # (cos)
  942 #
  943 # Compute cos(z) = (exp(iz) + exp(-iz))/2.
  944 #
  945 sub cos {
  946     my ($z) = @_ ? @_ : $_;
  947     return CORE::cos($z) unless ref $z;
  948     my ($x, $y) = @{$z->_cartesian};
  949     my $ey = CORE::exp($y);
  950     my $sx = CORE::sin($x);
  951     my $cx = CORE::cos($x);
  952     my $ey_1 = $ey ? 1 / $ey : Inf();
  953     return (ref $z)->make($cx * ($ey + $ey_1)/2,
  954                   $sx * ($ey_1 - $ey)/2);
  955 }
  956 
  957 #
  958 # (sin)
  959 #
  960 # Compute sin(z) = (exp(iz) - exp(-iz))/2.
  961 #
  962 sub sin {
  963     my ($z) = @_ ? @_ : $_;
  964     return CORE::sin($z) unless ref $z;
  965     my ($x, $y) = @{$z->_cartesian};
  966     my $ey = CORE::exp($y);
  967     my $sx = CORE::sin($x);
  968     my $cx = CORE::cos($x);
  969     my $ey_1 = $ey ? 1 / $ey : Inf();
  970     return (ref $z)->make($sx * ($ey + $ey_1)/2,
  971                   $cx * ($ey - $ey_1)/2);
  972 }
  973 
  974 #
  975 # tan
  976 #
  977 # Compute tan(z) = sin(z) / cos(z).
  978 #
  979 sub tan {
  980     my ($z) = @_;
  981     my $cz = &cos($z);
  982     _divbyzero "tan($z)", "cos($z)" if $cz == 0;
  983     return &sin($z) / $cz;
  984 }
  985 
  986 #
  987 # sec
  988 #
  989 # Computes the secant sec(z) = 1 / cos(z).
  990 #
  991 sub sec {
  992     my ($z) = @_;
  993     my $cz = &cos($z);
  994     _divbyzero "sec($z)", "cos($z)" if ($cz == 0);
  995     return 1 / $cz;
  996 }
  997 
  998 #
  999 # csc
 1000 #
 1001 # Computes the cosecant csc(z) = 1 / sin(z).
 1002 #
 1003 sub csc {
 1004     my ($z) = @_;
 1005     my $sz = &sin($z);
 1006     _divbyzero "csc($z)", "sin($z)" if ($sz == 0);
 1007     return 1 / $sz;
 1008 }
 1009 
 1010 #
 1011 # cosec
 1012 #
 1013 # Alias for csc().
 1014 #
 1015 sub cosec { Math::Complex::csc(@_) }
 1016 
 1017 #
 1018 # cot
 1019 #
 1020 # Computes cot(z) = cos(z) / sin(z).
 1021 #
 1022 sub cot {
 1023     my ($z) = @_;
 1024     my $sz = &sin($z);
 1025     _divbyzero "cot($z)", "sin($z)" if ($sz == 0);
 1026     return &cos($z) / $sz;
 1027 }
 1028 
 1029 #
 1030 # cotan
 1031 #
 1032 # Alias for cot().
 1033 #
 1034 sub cotan { Math::Complex::cot(@_) }
 1035 
 1036 #
 1037 # acos
 1038 #
 1039 # Computes the arc cosine acos(z) = -i log(z + sqrt(z*z-1)).
 1040 #
 1041 sub acos {
 1042     my $z = $_[0];
 1043     return CORE::atan2(CORE::sqrt(1-$z*$z), $z)
 1044         if (! ref $z) && CORE::abs($z) <= 1;
 1045     $z = cplx($z, 0) unless ref $z;
 1046     my ($x, $y) = @{$z->_cartesian};
 1047     return 0 if $x == 1 && $y == 0;
 1048     my $t1 = CORE::sqrt(($x+1)*($x+1) + $y*$y);
 1049     my $t2 = CORE::sqrt(($x-1)*($x-1) + $y*$y);
 1050     my $alpha = ($t1 + $t2)/2;
 1051     my $beta  = ($t1 - $t2)/2;
 1052     $alpha = 1 if $alpha < 1;
 1053     if    ($beta >  1) { $beta =  1 }
 1054     elsif ($beta < -1) { $beta = -1 }
 1055     my $u = CORE::atan2(CORE::sqrt(1-$beta*$beta), $beta);
 1056     my $v = CORE::log($alpha + CORE::sqrt($alpha*$alpha-1));
 1057     $v = -$v if $y > 0 || ($y == 0 && $x < -1);
 1058     return (ref $z)->make($u, $v);
 1059 }
 1060 
 1061 #
 1062 # asin
 1063 #
 1064 # Computes the arc sine asin(z) = -i log(iz + sqrt(1-z*z)).
 1065 #
 1066 sub asin {
 1067     my $z = $_[0];
 1068     return CORE::atan2($z, CORE::sqrt(1-$z*$z))
 1069         if (! ref $z) && CORE::abs($z) <= 1;
 1070     $z = cplx($z, 0) unless ref $z;
 1071     my ($x, $y) = @{$z->_cartesian};
 1072     return 0 if $x == 0 && $y == 0;
 1073     my $t1 = CORE::sqrt(($x+1)*($x+1) + $y*$y);
 1074     my $t2 = CORE::sqrt(($x-1)*($x-1) + $y*$y);
 1075     my $alpha = ($t1 + $t2)/2;
 1076     my $beta  = ($t1 - $t2)/2;
 1077     $alpha = 1 if $alpha < 1;
 1078     if    ($beta >  1) { $beta =  1 }
 1079     elsif ($beta < -1) { $beta = -1 }
 1080     my $u =  CORE::atan2($beta, CORE::sqrt(1-$beta*$beta));
 1081     my $v = -CORE::log($alpha + CORE::sqrt($alpha*$alpha-1));
 1082     $v = -$v if $y > 0 || ($y == 0 && $x < -1);
 1083     return (ref $z)->make($u, $v);
 1084 }
 1085 
 1086 #
 1087 # atan
 1088 #
 1089 # Computes the arc tangent atan(z) = i/2 log((i+z) / (i-z)).
 1090 #
 1091 sub atan {
 1092     my ($z) = @_;
 1093     return CORE::atan2($z, 1) unless ref $z;
 1094     my ($x, $y) = ref $z ? @{$z->_cartesian} : ($z, 0);
 1095     return 0 if $x == 0 && $y == 0;
 1096     _divbyzero "atan(i)"  if ( $z == i);
 1097     _logofzero "atan(-i)" if (-$z == i); # -i is a bad file test...
 1098     my $log = &log((i + $z) / (i - $z));
 1099     return _ip2 * $log;
 1100 }
 1101 
 1102 #
 1103 # asec
 1104 #
 1105 # Computes the arc secant asec(z) = acos(1 / z).
 1106 #
 1107 sub asec {
 1108     my ($z) = @_;
 1109     _divbyzero "asec($z)", $z if ($z == 0);
 1110     return acos(1 / $z);
 1111 }
 1112 
 1113 #
 1114 # acsc
 1115 #
 1116 # Computes the arc cosecant acsc(z) = asin(1 / z).
 1117 #
 1118 sub acsc {
 1119     my ($z) = @_;
 1120     _divbyzero "acsc($z)", $z if ($z == 0);
 1121     return asin(1 / $z);
 1122 }
 1123 
 1124 #
 1125 # acosec
 1126 #
 1127 # Alias for acsc().
 1128 #
 1129 sub acosec { Math::Complex::acsc(@_) }
 1130 
 1131 #
 1132 # acot
 1133 #
 1134 # Computes the arc cotangent acot(z) = atan(1 / z)
 1135 #
 1136 sub acot {
 1137     my ($z) = @_;
 1138     _divbyzero "acot(0)"  if $z == 0;
 1139     return ($z >= 0) ? CORE::atan2(1, $z) : CORE::atan2(-1, -$z)
 1140         unless ref $z;
 1141     _divbyzero "acot(i)"  if ($z - i == 0);
 1142     _logofzero "acot(-i)" if ($z + i == 0);
 1143     return atan(1 / $z);
 1144 }
 1145 
 1146 #
 1147 # acotan
 1148 #
 1149 # Alias for acot().
 1150 #
 1151 sub acotan { Math::Complex::acot(@_) }
 1152 
 1153 #
 1154 # cosh
 1155 #
 1156 # Computes the hyperbolic cosine cosh(z) = (exp(z) + exp(-z))/2.
 1157 #
 1158 sub cosh {
 1159     my ($z) = @_;
 1160     my $ex;
 1161     unless (ref $z) {
 1162         $ex = CORE::exp($z);
 1163             return $ex ? ($ex == $ExpInf ? Inf() : ($ex + 1/$ex)/2) : Inf();
 1164     }
 1165     my ($x, $y) = @{$z->_cartesian};
 1166     $ex = CORE::exp($x);
 1167     my $ex_1 = $ex ? 1 / $ex : Inf();
 1168     return (ref $z)->make(CORE::cos($y) * ($ex + $ex_1)/2,
 1169                   CORE::sin($y) * ($ex - $ex_1)/2);
 1170 }
 1171 
 1172 #
 1173 # sinh
 1174 #
 1175 # Computes the hyperbolic sine sinh(z) = (exp(z) - exp(-z))/2.
 1176 #
 1177 sub sinh {
 1178     my ($z) = @_;
 1179     my $ex;
 1180     unless (ref $z) {
 1181         return 0 if $z == 0;
 1182         $ex = CORE::exp($z);
 1183             return $ex ? ($ex == $ExpInf ? Inf() : ($ex - 1/$ex)/2) : -Inf();
 1184     }
 1185     my ($x, $y) = @{$z->_cartesian};
 1186     my $cy = CORE::cos($y);
 1187     my $sy = CORE::sin($y);
 1188     $ex = CORE::exp($x);
 1189     my $ex_1 = $ex ? 1 / $ex : Inf();
 1190     return (ref $z)->make(CORE::cos($y) * ($ex - $ex_1)/2,
 1191                   CORE::sin($y) * ($ex + $ex_1)/2);
 1192 }
 1193 
 1194 #
 1195 # tanh
 1196 #
 1197 # Computes the hyperbolic tangent tanh(z) = sinh(z) / cosh(z).
 1198 #
 1199 sub tanh {
 1200     my ($z) = @_;
 1201     my $cz = cosh($z);
 1202     _divbyzero "tanh($z)", "cosh($z)" if ($cz == 0);
 1203     my $sz = sinh($z);
 1204     return  1 if $cz ==  $sz;
 1205     return -1 if $cz == -$sz;
 1206     return $sz / $cz;
 1207 }
 1208 
 1209 #
 1210 # sech
 1211 #
 1212 # Computes the hyperbolic secant sech(z) = 1 / cosh(z).
 1213 #
 1214 sub sech {
 1215     my ($z) = @_;
 1216     my $cz = cosh($z);
 1217     _divbyzero "sech($z)", "cosh($z)" if ($cz == 0);
 1218     return 1 / $cz;
 1219 }
 1220 
 1221 #
 1222 # csch
 1223 #
 1224 # Computes the hyperbolic cosecant csch(z) = 1 / sinh(z).
 1225 #
 1226 sub csch {
 1227     my ($z) = @_;
 1228     my $sz = sinh($z);
 1229     _divbyzero "csch($z)", "sinh($z)" if ($sz == 0);
 1230     return 1 / $sz;
 1231 }
 1232 
 1233 #
 1234 # cosech
 1235 #
 1236 # Alias for csch().
 1237 #
 1238 sub cosech { Math::Complex::csch(@_) }
 1239 
 1240 #
 1241 # coth
 1242 #
 1243 # Computes the hyperbolic cotangent coth(z) = cosh(z) / sinh(z).
 1244 #
 1245 sub coth {
 1246     my ($z) = @_;
 1247     my $sz = sinh($z);
 1248     _divbyzero "coth($z)", "sinh($z)" if $sz == 0;
 1249     my $cz = cosh($z);
 1250     return  1 if $cz ==  $sz;
 1251     return -1 if $cz == -$sz;
 1252     return $cz / $sz;
 1253 }
 1254 
 1255 #
 1256 # cotanh
 1257 #
 1258 # Alias for coth().
 1259 #
 1260 sub cotanh { Math::Complex::coth(@_) }
 1261 
 1262 #
 1263 # acosh
 1264 #
 1265 # Computes the area/inverse hyperbolic cosine acosh(z) = log(z + sqrt(z*z-1)).
 1266 #
 1267 sub acosh {
 1268     my ($z) = @_;
 1269     unless (ref $z) {
 1270         $z = cplx($z, 0);
 1271     }
 1272     my ($re, $im) = @{$z->_cartesian};
 1273     if ($im == 0) {
 1274         return CORE::log($re + CORE::sqrt($re*$re - 1))
 1275         if $re >= 1;
 1276         return cplx(0, CORE::atan2(CORE::sqrt(1 - $re*$re), $re))
 1277         if CORE::abs($re) < 1;
 1278     }
 1279     my $t = &sqrt($z * $z - 1) + $z;
 1280     # Try Taylor if looking bad (this usually means that
 1281     # $z was large negative, therefore the sqrt is really
 1282     # close to abs(z), summing that with z...)
 1283     $t = 1/(2 * $z) - 1/(8 * $z**3) + 1/(16 * $z**5) - 5/(128 * $z**7)
 1284         if $t == 0;
 1285     my $u = &log($t);
 1286     $u->Im(-$u->Im) if $re < 0 && $im == 0;
 1287     return $re < 0 ? -$u : $u;
 1288 }
 1289 
 1290 #
 1291 # asinh
 1292 #
 1293 # Computes the area/inverse hyperbolic sine asinh(z) = log(z + sqrt(z*z+1))
 1294 #
 1295 sub asinh {
 1296     my ($z) = @_;
 1297     unless (ref $z) {
 1298         my $t = $z + CORE::sqrt($z*$z + 1);
 1299         return CORE::log($t) if $t;
 1300     }
 1301     my $t = &sqrt($z * $z + 1) + $z;
 1302     # Try Taylor if looking bad (this usually means that
 1303     # $z was large negative, therefore the sqrt is really
 1304     # close to abs(z), summing that with z...)
 1305     $t = 1/(2 * $z) - 1/(8 * $z**3) + 1/(16 * $z**5) - 5/(128 * $z**7)
 1306         if $t == 0;
 1307     return &log($t);
 1308 }
 1309 
 1310 #
 1311 # atanh
 1312 #
 1313 # Computes the area/inverse hyperbolic tangent atanh(z) = 1/2 log((1+z) / (1-z)).
 1314 #
 1315 sub atanh {
 1316     my ($z) = @_;
 1317     unless (ref $z) {
 1318         return CORE::log((1 + $z)/(1 - $z))/2 if CORE::abs($z) < 1;
 1319         $z = cplx($z, 0);
 1320     }
 1321     _divbyzero 'atanh(1)',  "1 - $z" if (1 - $z == 0);
 1322     _logofzero 'atanh(-1)'           if (1 + $z == 0);
 1323     return 0.5 * &log((1 + $z) / (1 - $z));
 1324 }
 1325 
 1326 #
 1327 # asech
 1328 #
 1329 # Computes the area/inverse hyperbolic secant asech(z) = acosh(1 / z).
 1330 #
 1331 sub asech {
 1332     my ($z) = @_;
 1333     _divbyzero 'asech(0)', "$z" if ($z == 0);
 1334     return acosh(1 / $z);
 1335 }
 1336 
 1337 #
 1338 # acsch
 1339 #
 1340 # Computes the area/inverse hyperbolic cosecant acsch(z) = asinh(1 / z).
 1341 #
 1342 sub acsch {
 1343     my ($z) = @_;
 1344     _divbyzero 'acsch(0)', $z if ($z == 0);
 1345     return asinh(1 / $z);
 1346 }
 1347 
 1348 #
 1349 # acosech
 1350 #
 1351 # Alias for acosh().
 1352 #
 1353 sub acosech { Math::Complex::acsch(@_) }
 1354 
 1355 #
 1356 # acoth
 1357 #
 1358 # Computes the area/inverse hyperbolic cotangent acoth(z) = 1/2 log((1+z) / (z-1)).
 1359 #
 1360 sub acoth {
 1361     my ($z) = @_;
 1362     _divbyzero 'acoth(0)'            if ($z == 0);
 1363     unless (ref $z) {
 1364         return CORE::log(($z + 1)/($z - 1))/2 if CORE::abs($z) > 1;
 1365         $z = cplx($z, 0);
 1366     }
 1367     _divbyzero 'acoth(1)',  "$z - 1" if ($z - 1 == 0);
 1368     _logofzero 'acoth(-1)', "1 + $z" if (1 + $z == 0);
 1369     return &log((1 + $z) / ($z - 1)) / 2;
 1370 }
 1371 
 1372 #
 1373 # acotanh
 1374 #
 1375 # Alias for acot().
 1376 #
 1377 sub acotanh { Math::Complex::acoth(@_) }
 1378 
 1379 #
 1380 # (atan2)
 1381 #
 1382 # Compute atan(z1/z2), minding the right quadrant.
 1383 #
 1384 sub atan2 {
 1385     my ($z1, $z2, $inverted) = @_;
 1386     my ($re1, $im1, $re2, $im2);
 1387     if ($inverted) {
 1388         ($re1, $im1) = ref $z2 ? @{$z2->_cartesian} : ($z2, 0);
 1389         ($re2, $im2) = ref $z1 ? @{$z1->_cartesian} : ($z1, 0);
 1390     } else {
 1391         ($re1, $im1) = ref $z1 ? @{$z1->_cartesian} : ($z1, 0);
 1392         ($re2, $im2) = ref $z2 ? @{$z2->_cartesian} : ($z2, 0);
 1393     }
 1394     if ($im1 || $im2) {
 1395         # In MATLAB the imaginary parts are ignored.
 1396         # warn "atan2: Imaginary parts ignored";
 1397         # http://documents.wolfram.com/mathematica/functions/ArcTan
 1398         # NOTE: Mathematica ArcTan[x,y] while atan2(y,x)
 1399         my $s = $z1 * $z1 + $z2 * $z2;
 1400         _divbyzero("atan2") if $s == 0;
 1401         my $i = &i;
 1402         my $r = $z2 + $z1 * $i;
 1403         return -$i * &log($r / &sqrt( $s ));
 1404     }
 1405     return CORE::atan2($re1, $re2);
 1406 }
 1407 
 1408 #
 1409 # display_format
 1410 # ->display_format
 1411 #
 1412 # Set (get if no argument) the display format for all complex numbers that
 1413 # don't happen to have overridden it via ->display_format
 1414 #
 1415 # When called as an object method, this actually sets the display format for
 1416 # the current object.
 1417 #
 1418 # Valid object formats are 'c' and 'p' for cartesian and polar. The first
 1419 # letter is used actually, so the type can be fully spelled out for clarity.
 1420 #
 1421 sub display_format {
 1422     my $self  = shift;
 1423     my %display_format = %DISPLAY_FORMAT;
 1424 
 1425     if (ref $self) {            # Called as an object method
 1426         if (exists $self->{display_format}) {
 1427         my %obj = %{$self->{display_format}};
 1428         @display_format{keys %obj} = values %obj;
 1429         }
 1430     }
 1431     if (@_ == 1) {
 1432         $display_format{style} = shift;
 1433     } else {
 1434         my %new = @_;
 1435         @display_format{keys %new} = values %new;
 1436     }
 1437 
 1438     if (ref $self) { # Called as an object method
 1439         $self->{display_format} = { %display_format };
 1440         return
 1441         wantarray ?
 1442             %{$self->{display_format}} :
 1443             $self->{display_format}->{style};
 1444     }
 1445 
 1446         # Called as a class method
 1447     %DISPLAY_FORMAT = %display_format;
 1448     return
 1449         wantarray ?
 1450         %DISPLAY_FORMAT :
 1451             $DISPLAY_FORMAT{style};
 1452 }
 1453 
 1454 #
 1455 # (_stringify)
 1456 #
 1457 # Show nicely formatted complex number under its cartesian or polar form,
 1458 # depending on the current display format:
 1459 #
 1460 # . If a specific display format has been recorded for this object, use it.
 1461 # . Otherwise, use the generic current default for all complex numbers,
 1462 #   which is a package global variable.
 1463 #
 1464 sub _stringify {
 1465     my ($z) = shift;
 1466 
 1467     my $style = $z->display_format;
 1468 
 1469     $style = $DISPLAY_FORMAT{style} unless defined $style;
 1470 
 1471     return $z->_stringify_polar if $style =~ /^p/i;
 1472     return $z->_stringify_cartesian;
 1473 }
 1474 
 1475 #
 1476 # ->_stringify_cartesian
 1477 #
 1478 # Stringify as a cartesian representation 'a+bi'.
 1479 #
 1480 sub _stringify_cartesian {
 1481     my $z  = shift;
 1482     my ($x, $y) = @{$z->_cartesian};
 1483     my ($re, $im);
 1484 
 1485     my %format = $z->display_format;
 1486     my $format = $format{format};
 1487 
 1488     if ($x) {
 1489         if ($x =~ /^NaN[QS]?$/i) {
 1490         $re = $x;
 1491         } else {
 1492         if ($x =~ /^-?\Q$Inf\E$/oi) {
 1493             $re = $x;
 1494         } else {
 1495             $re = defined $format ? sprintf($format, $x) : $x;
 1496         }
 1497         }
 1498     } else {
 1499         undef $re;
 1500     }
 1501 
 1502     if ($y) {
 1503         if ($y =~ /^(NaN[QS]?)$/i) {
 1504         $im = $y;
 1505         } else {
 1506         if ($y =~ /^-?\Q$Inf\E$/oi) {
 1507             $im = $y;
 1508         } else {
 1509             $im =
 1510             defined $format ?
 1511                 sprintf($format, $y) :
 1512                 ($y == 1 ? "" : ($y == -1 ? "-" : $y));
 1513         }
 1514         }
 1515         $im .= "i";
 1516     } else {
 1517         undef $im;
 1518     }
 1519 
 1520     my $str = $re;
 1521 
 1522     if (defined $im) {
 1523         if ($y < 0) {
 1524         $str .= $im;
 1525         } elsif ($y > 0 || $im =~ /^NaN[QS]?i$/i)  {
 1526         $str .= "+" if defined $re;
 1527         $str .= $im;
 1528         }
 1529     } elsif (!defined $re) {
 1530         $str = "0";
 1531     }
 1532 
 1533     return $str;
 1534 }
 1535 
 1536 
 1537 #
 1538 # ->_stringify_polar
 1539 #
 1540 # Stringify as a polar representation '[r,t]'.
 1541 #
 1542 sub _stringify_polar {
 1543     my $z  = shift;
 1544     my ($r, $t) = @{$z->_polar};
 1545     my $theta;
 1546 
 1547     my %format = $z->display_format;
 1548     my $format = $format{format};
 1549 
 1550     if ($t =~ /^NaN[QS]?$/i || $t =~ /^-?\Q$Inf\E$/oi) {
 1551         $theta = $t; 
 1552     } elsif ($t == pi) {
 1553         $theta = "pi";
 1554     } elsif ($r == 0 || $t == 0) {
 1555         $theta = defined $format ? sprintf($format, $t) : $t;
 1556     }
 1557 
 1558     return "[$r,$theta]" if defined $theta;
 1559 
 1560     #
 1561     # Try to identify pi/n and friends.
 1562     #
 1563 
 1564     $t -= int(CORE::abs($t) / pi2) * pi2;
 1565 
 1566     if ($format{polar_pretty_print} && $t) {
 1567         my ($a, $b);
 1568         for $a (2..9) {
 1569         $b = $t * $a / pi;
 1570         if ($b =~ /^-?\d+$/) {
 1571             $b = $b < 0 ? "-" : "" if CORE::abs($b) == 1;
 1572             $theta = "${b}pi/$a";
 1573             last;
 1574         }
 1575         }
 1576     }
 1577 
 1578         if (defined $format) {
 1579         $r     = sprintf($format, $r);
 1580         $theta = sprintf($format, $t) unless defined $theta;
 1581     } else {
 1582         $theta = $t unless defined $theta;
 1583     }
 1584 
 1585     return "[$r,$theta]";
 1586 }
 1587 
 1588 sub Inf {
 1589     return $Inf;
 1590 }
 1591 
 1592 1;
 1593 __END__
 1594 
 1595 =pod
 1596 
 1597 =head1 NAME
 1598 
 1599 Math::Complex - complex numbers and associated mathematical functions
 1600 
 1601 =head1 SYNOPSIS
 1602 
 1603     use Math::Complex;
 1604 
 1605     $z = Math::Complex->make(5, 6);
 1606     $t = 4 - 3*i + $z;
 1607     $j = cplxe(1, 2*pi/3);
 1608 
 1609 =head1 DESCRIPTION
 1610 
 1611 This package lets you create and manipulate complex numbers. By default,
 1612 I<Perl> limits itself to real numbers, but an extra C<use> statement brings
 1613 full complex support, along with a full set of mathematical functions
 1614 typically associated with and/or extended to complex numbers.
 1615 
 1616 If you wonder what complex numbers are, they were invented to be able to solve
 1617 the following equation:
 1618 
 1619     x*x = -1
 1620 
 1621 and by definition, the solution is noted I<i> (engineers use I<j> instead since
 1622 I<i> usually denotes an intensity, but the name does not matter). The number
 1623 I<i> is a pure I<imaginary> number.
 1624 
 1625 The arithmetics with pure imaginary numbers works just like you would expect
 1626 it with real numbers... you just have to remember that
 1627 
 1628     i*i = -1
 1629 
 1630 so you have:
 1631 
 1632     5i + 7i = i * (5 + 7) = 12i
 1633     4i - 3i = i * (4 - 3) = i
 1634     4i * 2i = -8
 1635     6i / 2i = 3
 1636     1 / i = -i
 1637 
 1638 Complex numbers are numbers that have both a real part and an imaginary
 1639 part, and are usually noted:
 1640 
 1641     a + bi
 1642 
 1643 where C<a> is the I<real> part and C<b> is the I<imaginary> part. The
 1644 arithmetic with complex numbers is straightforward. You have to
 1645 keep track of the real and the imaginary parts, but otherwise the
 1646 rules used for real numbers just apply:
 1647 
 1648     (4 + 3i) + (5 - 2i) = (4 + 5) + i(3 - 2) = 9 + i
 1649     (2 + i) * (4 - i) = 2*4 + 4i -2i -i*i = 8 + 2i + 1 = 9 + 2i
 1650 
 1651 A graphical representation of complex numbers is possible in a plane
 1652 (also called the I<complex plane>, but it's really a 2D plane).
 1653 The number
 1654 
 1655     z = a + bi
 1656 
 1657 is the point whose coordinates are (a, b). Actually, it would
 1658 be the vector originating from (0, 0) to (a, b). It follows that the addition
 1659 of two complex numbers is a vectorial addition.
 1660 
 1661 Since there is a bijection between a point in the 2D plane and a complex
 1662 number (i.e. the mapping is unique and reciprocal), a complex number
 1663 can also be uniquely identified with polar coordinates:
 1664 
 1665     [rho, theta]
 1666 
 1667 where C<rho> is the distance to the origin, and C<theta> the angle between
 1668 the vector and the I<x> axis. There is a notation for this using the
 1669 exponential form, which is:
 1670 
 1671     rho * exp(i * theta)
 1672 
 1673 where I<i> is the famous imaginary number introduced above. Conversion
 1674 between this form and the cartesian form C<a + bi> is immediate:
 1675 
 1676     a = rho * cos(theta)
 1677     b = rho * sin(theta)
 1678 
 1679 which is also expressed by this formula:
 1680 
 1681     z = rho * exp(i * theta) = rho * (cos theta + i * sin theta)
 1682 
 1683 In other words, it's the projection of the vector onto the I<x> and I<y>
 1684 axes. Mathematicians call I<rho> the I<norm> or I<modulus> and I<theta>
 1685 the I<argument> of the complex number. The I<norm> of C<z> is
 1686 marked here as C<abs(z)>.
 1687 
 1688 The polar notation (also known as the trigonometric representation) is
 1689 much more handy for performing multiplications and divisions of
 1690 complex numbers, whilst the cartesian notation is better suited for
 1691 additions and subtractions. Real numbers are on the I<x> axis, and
 1692 therefore I<y> or I<theta> is zero or I<pi>.
 1693 
 1694 All the common operations that can be performed on a real number have
 1695 been defined to work on complex numbers as well, and are merely
 1696 I<extensions> of the operations defined on real numbers. This means
 1697 they keep their natural meaning when there is no imaginary part, provided
 1698 the number is within their definition set.
 1699 
 1700 For instance, the C<sqrt> routine which computes the square root of
 1701 its argument is only defined for non-negative real numbers and yields a
 1702 non-negative real number (it is an application from B<R+> to B<R+>).
 1703 If we allow it to return a complex number, then it can be extended to
 1704 negative real numbers to become an application from B<R> to B<C> (the
 1705 set of complex numbers):
 1706 
 1707     sqrt(x) = x >= 0 ? sqrt(x) : sqrt(-x)*i
 1708 
 1709 It can also be extended to be an application from B<C> to B<C>,
 1710 whilst its restriction to B<R> behaves as defined above by using
 1711 the following definition:
 1712 
 1713     sqrt(z = [r,t]) = sqrt(r) * exp(i * t/2)
 1714 
 1715 Indeed, a negative real number can be noted C<[x,pi]> (the modulus
 1716 I<x> is always non-negative, so C<[x,pi]> is really C<-x>, a negative
 1717 number) and the above definition states that
 1718 
 1719     sqrt([x,pi]) = sqrt(x) * exp(i*pi/2) = [sqrt(x),pi/2] = sqrt(x)*i
 1720 
 1721 which is exactly what we had defined for negative real numbers above.
 1722 The C<sqrt> returns only one of the solutions: if you want the both,
 1723 use the C<root> function.
 1724 
 1725 All the common mathematical functions defined on real numbers that
 1726 are extended to complex numbers share that same property of working
 1727 I<as usual> when the imaginary part is zero (otherwise, it would not
 1728 be called an extension, would it?).
 1729 
 1730 A I<new> operation possible on a complex number that is
 1731 the identity for real numbers is called the I<conjugate>, and is noted
 1732 with a horizontal bar above the number, or C<~z> here.
 1733 
 1734      z = a + bi
 1735     ~z = a - bi
 1736 
 1737 Simple... Now look:
 1738 
 1739     z * ~z = (a + bi) * (a - bi) = a*a + b*b
 1740 
 1741 We saw that the norm of C<z> was noted C<abs(z)> and was defined as the
 1742 distance to the origin, also known as:
 1743 
 1744     rho = abs(z) = sqrt(a*a + b*b)
 1745 
 1746 so
 1747 
 1748     z * ~z = abs(z) ** 2
 1749 
 1750 If z is a pure real number (i.e. C<b == 0>), then the above yields:
 1751 
 1752     a * a = abs(a) ** 2
 1753 
 1754 which is true (C<abs> has the regular meaning for real number, i.e. stands
 1755 for the absolute value). This example explains why the norm of C<z> is
 1756 noted C<abs(z)>: it extends the C<abs> function to complex numbers, yet
 1757 is the regular C<abs> we know when the complex number actually has no
 1758 imaginary part... This justifies I<a posteriori> our use of the C<abs>
 1759 notation for the norm.
 1760 
 1761 =head1 OPERATIONS
 1762 
 1763 Given the following notations:
 1764 
 1765     z1 = a + bi = r1 * exp(i * t1)
 1766     z2 = c + di = r2 * exp(i * t2)
 1767     z = <any complex or real number>
 1768 
 1769 the following (overloaded) operations are supported on complex numbers:
 1770 
 1771     z1 + z2 = (a + c) + i(b + d)
 1772     z1 - z2 = (a - c) + i(b - d)
 1773     z1 * z2 = (r1 * r2) * exp(i * (t1 + t2))
 1774     z1 / z2 = (r1 / r2) * exp(i * (t1 - t2))
 1775     z1 ** z2 = exp(z2 * log z1)
 1776     ~z = a - bi
 1777     abs(z) = r1 = sqrt(a*a + b*b)
 1778     sqrt(z) = sqrt(r1) * exp(i * t/2)
 1779     exp(z) = exp(a) * exp(i * b)
 1780     log(z) = log(r1) + i*t
 1781     sin(z) = 1/2i (exp(i * z1) - exp(-i * z))
 1782     cos(z) = 1/2 (exp(i * z1) + exp(-i * z))
 1783     atan2(y, x) = atan(y / x) # Minding the right quadrant, note the order.
 1784 
 1785 The definition used for complex arguments of atan2() is
 1786 
 1787        -i log((x + iy)/sqrt(x*x+y*y))
 1788 
 1789 Note that atan2(0, 0) is not well-defined.
 1790 
 1791 The following extra operations are supported on both real and complex
 1792 numbers:
 1793 
 1794     Re(z) = a
 1795     Im(z) = b
 1796     arg(z) = t
 1797     abs(z) = r
 1798 
 1799     cbrt(z) = z ** (1/3)
 1800     log10(z) = log(z) / log(10)
 1801     logn(z, n) = log(z) / log(n)
 1802 
 1803     tan(z) = sin(z) / cos(z)
 1804 
 1805     csc(z) = 1 / sin(z)
 1806     sec(z) = 1 / cos(z)
 1807     cot(z) = 1 / tan(z)
 1808 
 1809     asin(z) = -i * log(i*z + sqrt(1-z*z))
 1810     acos(z) = -i * log(z + i*sqrt(1-z*z))
 1811     atan(z) = i/2 * log((i+z) / (i-z))
 1812 
 1813     acsc(z) = asin(1 / z)
 1814     asec(z) = acos(1 / z)
 1815     acot(z) = atan(1 / z) = -i/2 * log((i+z) / (z-i))
 1816 
 1817     sinh(z) = 1/2 (exp(z) - exp(-z))
 1818     cosh(z) = 1/2 (exp(z) + exp(-z))
 1819     tanh(z) = sinh(z) / cosh(z) = (exp(z) - exp(-z)) / (exp(z) + exp(-z))
 1820 
 1821     csch(z) = 1 / sinh(z)
 1822     sech(z) = 1 / cosh(z)
 1823     coth(z) = 1 / tanh(z)
 1824 
 1825     asinh(z) = log(z + sqrt(z*z+1))
 1826     acosh(z) = log(z + sqrt(z*z-1))
 1827     atanh(z) = 1/2 * log((1+z) / (1-z))
 1828 
 1829     acsch(z) = asinh(1 / z)
 1830     asech(z) = acosh(1 / z)
 1831     acoth(z) = atanh(1 / z) = 1/2 * log((1+z) / (z-1))
 1832 
 1833 I<arg>, I<abs>, I<log>, I<csc>, I<cot>, I<acsc>, I<acot>, I<csch>,
 1834 I<coth>, I<acosech>, I<acotanh>, have aliases I<rho>, I<theta>, I<ln>,
 1835 I<cosec>, I<cotan>, I<acosec>, I<acotan>, I<cosech>, I<cotanh>,
 1836 I<acosech>, I<acotanh>, respectively.  C<Re>, C<Im>, C<arg>, C<abs>,
 1837 C<rho>, and C<theta> can be used also as mutators.  The C<cbrt>
 1838 returns only one of the solutions: if you want all three, use the
 1839 C<root> function.
 1840 
 1841 The I<root> function is available to compute all the I<n>
 1842 roots of some complex, where I<n> is a strictly positive integer.
 1843 There are exactly I<n> such roots, returned as a list. Getting the
 1844 number mathematicians call C<j> such that:
 1845 
 1846     1 + j + j*j = 0;
 1847 
 1848 is a simple matter of writing:
 1849 
 1850     $j = ((root(1, 3))[1];
 1851 
 1852 The I<k>th root for C<z = [r,t]> is given by:
 1853 
 1854     (root(z, n))[k] = r**(1/n) * exp(i * (t + 2*k*pi)/n)
 1855 
 1856 You can return the I<k>th root directly by C<root(z, n, k)>,
 1857 indexing starting from I<zero> and ending at I<n - 1>.
 1858 
 1859 The I<spaceship> numeric comparison operator, E<lt>=E<gt>, is also
 1860 defined. In order to ensure its restriction to real numbers is conform
 1861 to what you would expect, the comparison is run on the real part of
 1862 the complex number first, and imaginary parts are compared only when
 1863 the real parts match.
 1864 
 1865 =head1 CREATION
 1866 
 1867 To create a complex number, use either:
 1868 
 1869     $z = Math::Complex->make(3, 4);
 1870     $z = cplx(3, 4);
 1871 
 1872 if you know the cartesian form of the number, or
 1873 
 1874     $z = 3 + 4*i;
 1875 
 1876 if you like. To create a number using the polar form, use either:
 1877 
 1878     $z = Math::Complex->emake(5, pi/3);
 1879     $x = cplxe(5, pi/3);
 1880 
 1881 instead. The first argument is the modulus, the second is the angle
 1882 (in radians, the full circle is 2*pi).  (Mnemonic: C<e> is used as a
 1883 notation for complex numbers in the polar form).
 1884 
 1885 It is possible to write:
 1886 
 1887     $x = cplxe(-3, pi/4);
 1888 
 1889 but that will be silently converted into C<[3,-3pi/4]>, since the
 1890 modulus must be non-negative (it represents the distance to the origin
 1891 in the complex plane).
 1892 
 1893 It is also possible to have a complex number as either argument of the
 1894 C<make>, C<emake>, C<cplx>, and C<cplxe>: the appropriate component of
 1895 the argument will be used.
 1896 
 1897     $z1 = cplx(-2,  1);
 1898     $z2 = cplx($z1, 4);
 1899 
 1900 The C<new>, C<make>, C<emake>, C<cplx>, and C<cplxe> will also
 1901 understand a single (string) argument of the forms
 1902 
 1903         2-3i
 1904         -3i
 1905     [2,3]
 1906     [2,-3pi/4]
 1907     [2]
 1908 
 1909 in which case the appropriate cartesian and exponential components
 1910 will be parsed from the string and used to create new complex numbers.
 1911 The imaginary component and the theta, respectively, will default to zero.
 1912 
 1913 The C<new>, C<make>, C<emake>, C<cplx>, and C<cplxe> will also
 1914 understand the case of no arguments: this means plain zero or (0, 0).
 1915 
 1916 =head1 DISPLAYING
 1917 
 1918 When printed, a complex number is usually shown under its cartesian
 1919 style I<a+bi>, but there are legitimate cases where the polar style
 1920 I<[r,t]> is more appropriate.  The process of converting the complex
 1921 number into a string that can be displayed is known as I<stringification>.
 1922 
 1923 By calling the class method C<Math::Complex::display_format> and
 1924 supplying either C<"polar"> or C<"cartesian"> as an argument, you
 1925 override the default display style, which is C<"cartesian">. Not
 1926 supplying any argument returns the current settings.
 1927 
 1928 This default can be overridden on a per-number basis by calling the
 1929 C<display_format> method instead. As before, not supplying any argument
 1930 returns the current display style for this number. Otherwise whatever you
 1931 specify will be the new display style for I<this> particular number.
 1932 
 1933 For instance:
 1934 
 1935     use Math::Complex;
 1936 
 1937     Math::Complex::display_format('polar');
 1938     $j = (root(1, 3))[1];
 1939     print "j = $j\n";       # Prints "j = [1,2pi/3]"
 1940     $j->display_format('cartesian');
 1941     print "j = $j\n";       # Prints "j = -0.5+0.866025403784439i"
 1942 
 1943 The polar style attempts to emphasize arguments like I<k*pi/n>
 1944 (where I<n> is a positive integer and I<k> an integer within [-9, +9]),
 1945 this is called I<polar pretty-printing>.
 1946 
 1947 For the reverse of stringifying, see the C<make> and C<emake>.
 1948 
 1949 =head2 CHANGED IN PERL 5.6
 1950 
 1951 The C<display_format> class method and the corresponding
 1952 C<display_format> object method can now be called using
 1953 a parameter hash instead of just a one parameter.
 1954 
 1955 The old display format style, which can have values C<"cartesian"> or
 1956 C<"polar">, can be changed using the C<"style"> parameter.
 1957 
 1958     $j->display_format(style => "polar");
 1959 
 1960 The one parameter calling convention also still works.
 1961 
 1962     $j->display_format("polar");
 1963 
 1964 There are two new display parameters.
 1965 
 1966 The first one is C<"format">, which is a sprintf()-style format string
 1967 to be used for both numeric parts of the complex number(s).  The is
 1968 somewhat system-dependent but most often it corresponds to C<"%.15g">.
 1969 You can revert to the default by setting the C<format> to C<undef>.
 1970 
 1971     # the $j from the above example
 1972 
 1973     $j->display_format('format' => '%.5f');
 1974     print "j = $j\n";       # Prints "j = -0.50000+0.86603i"
 1975     $j->display_format('format' => undef);
 1976     print "j = $j\n";       # Prints "j = -0.5+0.86603i"
 1977 
 1978 Notice that this affects also the return values of the
 1979 C<display_format> methods: in list context the whole parameter hash
 1980 will be returned, as opposed to only the style parameter value.
 1981 This is a potential incompatibility with earlier versions if you
 1982 have been calling the C<display_format> method in list context.
 1983 
 1984 The second new display parameter is C<"polar_pretty_print">, which can
 1985 be set to true or false, the default being true.  See the previous
 1986 section for what this means.
 1987 
 1988 =head1 USAGE
 1989 
 1990 Thanks to overloading, the handling of arithmetics with complex numbers
 1991 is simple and almost transparent.
 1992 
 1993 Here are some examples:
 1994 
 1995     use Math::Complex;
 1996 
 1997     $j = cplxe(1, 2*pi/3);  # $j ** 3 == 1
 1998     print "j = $j, j**3 = ", $j ** 3, "\n";
 1999     print "1 + j + j**2 = ", 1 + $j + $j**2, "\n";
 2000 
 2001     $z = -16 + 0*i;         # Force it to be a complex
 2002     print "sqrt($z) = ", sqrt($z), "\n";
 2003 
 2004     $k = exp(i * 2*pi/3);
 2005     print "$j - $k = ", $j - $k, "\n";
 2006 
 2007     $z->Re(3);          # Re, Im, arg, abs,
 2008     $j->arg(2);         # (the last two aka rho, theta)
 2009                     # can be used also as mutators.
 2010 
 2011 =head1 CONSTANTS
 2012 
 2013 =head2 PI
 2014 
 2015 The constant C<pi> and some handy multiples of it (pi2, pi4,
 2016 and pip2 (pi/2) and pip4 (pi/4)) are also available if separately
 2017 exported:
 2018 
 2019     use Math::Complex ':pi'; 
 2020     $third_of_circle = pi2 / 3;
 2021 
 2022 =head2 Inf
 2023 
 2024 The floating point infinity can be exported as a subroutine Inf():
 2025 
 2026     use Math::Complex qw(Inf sinh);
 2027     my $AlsoInf = Inf() + 42;
 2028     my $AnotherInf = sinh(1e42);
 2029     print "$AlsoInf is $AnotherInf\n" if $AlsoInf == $AnotherInf;
 2030 
 2031 Note that the stringified form of infinity varies between platforms:
 2032 it can be for example any of
 2033 
 2034    inf
 2035    infinity
 2036    INF
 2037    1.#INF
 2038 
 2039 or it can be something else. 
 2040 
 2041 Also note that in some platforms trying to use the infinity in
 2042 arithmetic operations may result in Perl crashing because using
 2043 an infinity causes SIGFPE or its moral equivalent to be sent.
 2044 The way to ignore this is
 2045 
 2046   local $SIG{FPE} = sub { };
 2047 
 2048 =head1 ERRORS DUE TO DIVISION BY ZERO OR LOGARITHM OF ZERO
 2049 
 2050 The division (/) and the following functions
 2051 
 2052     log ln  log10   logn
 2053     tan sec csc cot
 2054     atan    asec    acsc    acot
 2055     tanh    sech    csch    coth
 2056     atanh   asech   acsch   acoth
 2057 
 2058 cannot be computed for all arguments because that would mean dividing
 2059 by zero or taking logarithm of zero. These situations cause fatal
 2060 runtime errors looking like this
 2061 
 2062     cot(0): Division by zero.
 2063     (Because in the definition of cot(0), the divisor sin(0) is 0)
 2064     Died at ...
 2065 
 2066 or
 2067 
 2068     atanh(-1): Logarithm of zero.
 2069     Died at...
 2070 
 2071 For the C<csc>, C<cot>, C<asec>, C<acsc>, C<acot>, C<csch>, C<coth>,
 2072 C<asech>, C<acsch>, the argument cannot be C<0> (zero).  For the
 2073 logarithmic functions and the C<atanh>, C<acoth>, the argument cannot
 2074 be C<1> (one).  For the C<atanh>, C<acoth>, the argument cannot be
 2075 C<-1> (minus one).  For the C<atan>, C<acot>, the argument cannot be
 2076 C<i> (the imaginary unit).  For the C<atan>, C<acoth>, the argument
 2077 cannot be C<-i> (the negative imaginary unit).  For the C<tan>,
 2078 C<sec>, C<tanh>, the argument cannot be I<pi/2 + k * pi>, where I<k>
 2079 is any integer.  atan2(0, 0) is undefined, and if the complex arguments
 2080 are used for atan2(), a division by zero will happen if z1**2+z2**2 == 0.
 2081 
 2082 Note that because we are operating on approximations of real numbers,
 2083 these errors can happen when merely `too close' to the singularities
 2084 listed above.
 2085 
 2086 =head1 ERRORS DUE TO INDIGESTIBLE ARGUMENTS
 2087 
 2088 The C<make> and C<emake> accept both real and complex arguments.
 2089 When they cannot recognize the arguments they will die with error
 2090 messages like the following
 2091 
 2092     Math::Complex::make: Cannot take real part of ...
 2093     Math::Complex::make: Cannot take real part of ...
 2094     Math::Complex::emake: Cannot take rho of ...
 2095     Math::Complex::emake: Cannot take theta of ...
 2096 
 2097 =head1 BUGS
 2098 
 2099 Saying C<use Math::Complex;> exports many mathematical routines in the
 2100 caller environment and even overrides some (C<sqrt>, C<log>, C<atan2>).
 2101 This is construed as a feature by the Authors, actually... ;-)
 2102 
 2103 All routines expect to be given real or complex numbers. Don't attempt to
 2104 use BigFloat, since Perl has currently no rule to disambiguate a '+'
 2105 operation (for instance) between two overloaded entities.
 2106 
 2107 In Cray UNICOS there is some strange numerical instability that results
 2108 in root(), cos(), sin(), cosh(), sinh(), losing accuracy fast.  Beware.
 2109 The bug may be in UNICOS math libs, in UNICOS C compiler, in Math::Complex.
 2110 Whatever it is, it does not manifest itself anywhere else where Perl runs.
 2111 
 2112 =head1 SEE ALSO
 2113 
 2114 L<Math::Trig>
 2115 
 2116 =head1 AUTHORS
 2117 
 2118 Daniel S. Lewart <F<lewart!at!uiuc.edu>>,
 2119 Jarkko Hietaniemi <F<jhi!at!iki.fi>>,
 2120 Raphael Manfredi <F<Raphael_Manfredi!at!pobox.com>>,
 2121 Zefram <zefram@fysh.org>
 2122 
 2123 =head1 LICENSE
 2124 
 2125 This library is free software; you can redistribute it and/or modify
 2126 it under the same terms as Perl itself. 
 2127 
 2128 =cut
 2129 
 2130 1;
 2131 
 2132 # eof