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