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