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