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