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