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