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