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