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