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