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         $z = cplx($z, 0) unless ref $z;
 878         my ($x, $y) = @{$z->cartesian};
 879         return 0 if $x == 1 && $y == 0;
 880         my $t1 = CORE::sqrt(($x+1)*($x+1) + $y*$y);
 881         my $t2 = CORE::sqrt(($x-1)*($x-1) + $y*$y);
 882         my $alpha = ($t1 + $t2)/2;
 883         my $beta  = ($t1 - $t2)/2;
 884         $alpha = 1 if $alpha < 1;
 885         if    ($beta >  1) { $beta =  1 }
 886         elsif ($beta < -1) { $beta = -1 }
 887         my $u = CORE::atan2(CORE::sqrt(1-$beta*$beta), $beta);
 888         my $v = CORE::log($alpha + CORE::sqrt($alpha*$alpha-1));
 889         $v = -$v if $y > 0 || ($y == 0 && $x < -1);
 890         return (ref $z)->make($u, $v);
 891 }
 892
 893 #
 894 # asin
 895 #
 896 # Computes the arc sine asin(z) = -i log(iz + sqrt(1-z*z)).
 897 #
 898 sub asin {
 899         my $z = $_[0];
 900         return CORE::atan2($z, CORE::sqrt(1-$z*$z))
 901             if (! ref $z) && CORE::abs($z) <= 1;
 902         $z = cplx($z, 0) unless ref $z;
 903         my ($x, $y) = @{$z->cartesian};
 904         return 0 if $x == 0 && $y == 0;
 905         my $t1 = CORE::sqrt(($x+1)*($x+1) + $y*$y);
 906         my $t2 = CORE::sqrt(($x-1)*($x-1) + $y*$y);
 907         my $alpha = ($t1 + $t2)/2;
 908         my $beta  = ($t1 - $t2)/2;
 909         $alpha = 1 if $alpha < 1;
 910         if    ($beta >  1) { $beta =  1 }
 911         elsif ($beta < -1) { $beta = -1 }
 912         my $u =  CORE::atan2($beta, CORE::sqrt(1-$beta*$beta));
 913         my $v = -CORE::log($alpha + CORE::sqrt($alpha*$alpha-1));
 914         $v = -$v if $y > 0 || ($y == 0 && $x < -1);
 915         return (ref $z)->make($u, $v);
 916 }
 917
 918 #
 919 # atan
 920 #
 921 # Computes the arc tangent atan(z) = i/2 log((i+z) / (i-z)).
 922 #
 923 sub atan {
 924         my ($z) = @_;
 925         return CORE::atan2($z, 1) unless ref $z;
 926         my ($x, $y) = ref $z ? @{$z->cartesian} : ($z, 0);
 927         return 0 if $x == 0 && $y == 0;
 928         _divbyzero "atan(i)"  if ( $z == i);
 929         _logofzero "atan(-i)" if (-$z == i); # -i is a bad file test...
 930         my $log = &log((i + $z) / (i - $z));
 931         return ip2 * $log;
 932 }
 933
 934 #
 935 # asec
 936 #
 937 # Computes the arc secant asec(z) = acos(1 / z).
 938 #
 939 sub asec {
 940         my ($z) = @_;
 941         _divbyzero "asec($z)", $z if ($z == 0);
 942         return acos(1 / $z);
 943 }
 944
 945 #
 946 # acsc
 947 #
 948 # Computes the arc cosecant acsc(z) = asin(1 / z).
 949 #
 950 sub acsc {
 951         my ($z) = @_;
 952         _divbyzero "acsc($z)", $z if ($z == 0);
 953         return asin(1 / $z);
 954 }
 955
 956 #
 957 # acosec
 958 #
 959 # Alias for acsc().
 960 #
 961 sub acosec { Math::Complex::acsc(@_) }
 962
 963 #
 964 # acot
 965 #
 966 # Computes the arc cotangent acot(z) = atan(1 / z)
 967 #
 968 sub acot {
 969         my ($z) = @_;
 970         _divbyzero "acot(0)"  if $z == 0;
 971         return ($z >= 0) ? CORE::atan2(1, $z) : CORE::atan2(-1, -$z)
 972             unless ref $z;
 973         _divbyzero "acot(i)"  if ($z - i == 0);
 974         _logofzero "acot(-i)" if ($z + i == 0);
 975         return atan(1 / $z);
 976 }
 977
 978 #
 979 # acotan
 980 #
 981 # Alias for acot().
 982 #
 983 sub acotan { Math::Complex::acot(@_) }
 984
 985 #
 986 # cosh
 987 #
 988 # Computes the hyperbolic cosine cosh(z) = (exp(z) + exp(-z))/2.
 989 #
 990 sub cosh {
 991         my ($z) = @_;
 992         my $ex;
 993         unless (ref $z) {
 994             $ex = CORE::exp($z);
 995             return $ex ? ($ex + 1/$ex)/2 : $Inf;
 996         }
 997         my ($x, $y) = @{$z->cartesian};
 998         $ex = CORE::exp($x);
 999         my $ex_1 = $ex ? 1 / $ex : $Inf;
1000         return (ref $z)->make(CORE::cos($y) * ($ex + $ex_1)/2,
1001                               CORE::sin($y) * ($ex - $ex_1)/2);
1002 }
1003
1004 #
1005 # sinh
1006 #
1007 # Computes the hyperbolic sine sinh(z) = (exp(z) - exp(-z))/2.
1008 #
1009 sub sinh {
1010         my ($z) = @_;
1011         my $ex;
1012         unless (ref $z) {
1013             return 0 if $z == 0;
1014             $ex = CORE::exp($z);
1015             return $ex ? ($ex - 1/$ex)/2 : "-$Inf";
1016         }
1017         my ($x, $y) = @{$z->cartesian};
1018         my $cy = CORE::cos($y);
1019         my $sy = CORE::sin($y);
1020         $ex = CORE::exp($x);
1021         my $ex_1 = $ex ? 1 / $ex : $Inf;
1022         return (ref $z)->make(CORE::cos($y) * ($ex - $ex_1)/2,
1023                               CORE::sin($y) * ($ex + $ex_1)/2);
1024 }
1025
1026 #
1027 # tanh
1028 #
1029 # Computes the hyperbolic tangent tanh(z) = sinh(z) / cosh(z).
1030 #
1031 sub tanh {
1032         my ($z) = @_;
1033         my $cz = cosh($z);
1034         _divbyzero "tanh($z)", "cosh($z)" if ($cz == 0);
1035         return sinh($z) / $cz;
1036 }
1037
1038 #
1039 # sech
1040 #
1041 # Computes the hyperbolic secant sech(z) = 1 / cosh(z).
1042 #
1043 sub sech {
1044         my ($z) = @_;
1045         my $cz = cosh($z);
1046         _divbyzero "sech($z)", "cosh($z)" if ($cz == 0);
1047         return 1 / $cz;
1048 }
1049
1050 #
1051 # csch
1052 #
1053 # Computes the hyperbolic cosecant csch(z) = 1 / sinh(z).
1054 #
1055 sub csch {
1056         my ($z) = @_;
1057         my $sz = sinh($z);
1058         _divbyzero "csch($z)", "sinh($z)" if ($sz == 0);
1059         return 1 / $sz;
1060 }
1061
1062 #
1063 # cosech
1064 #
1065 # Alias for csch().
1066 #
1067 sub cosech { Math::Complex::csch(@_) }
1068
1069 #
1070 # coth
1071 #
1072 # Computes the hyperbolic cotangent coth(z) = cosh(z) / sinh(z).
1073 #
1074 sub coth {
1075         my ($z) = @_;
1076         my $sz = sinh($z);
1077         _divbyzero "coth($z)", "sinh($z)" if $sz == 0;
1078         return cosh($z) / $sz;
1079 }
1080
1081 #
1082 # cotanh
1083 #
1084 # Alias for coth().
1085 #
1086 sub cotanh { Math::Complex::coth(@_) }
1087
1088 #
1089 # acosh
1090 #
1091 # Computes the arc hyperbolic cosine acosh(z) = log(z + sqrt(z*z-1)).
1092 #
1093 sub acosh {
1094         my ($z) = @_;
1095         unless (ref $z) {
1096             $z = cplx($z, 0);
1097         }
1098         my ($re, $im) = @{$z->cartesian};
1099         if ($im == 0) {
1100             return CORE::log($re + CORE::sqrt($re*$re - 1))
1101                 if $re >= 1;
1102             return cplx(0, CORE::atan2(CORE::sqrt(1 - $re*$re), $re))
1103                 if CORE::abs($re) < 1;
1104         }
1105         my $t = &sqrt($z * $z - 1) + $z;
1106         # Try Taylor if looking bad (this usually means that
1107         # $z was large negative, therefore the sqrt is really
1108         # close to abs(z), summing that with z...)
1109         $t = 1/(2 * $z) - 1/(8 * $z**3) + 1/(16 * $z**5) - 5/(128 * $z**7)
1110             if $t == 0;
1111         my $u = &log($t);
1112         $u->Im(-$u->Im) if $re < 0 && $im == 0;
1113         return $re < 0 ? -$u : $u;
1114 }
1115
1116 #
1117 # asinh
1118 #
1119 # Computes the arc hyperbolic sine asinh(z) = log(z + sqrt(z*z+1))
1120 #
1121 sub asinh {
1122         my ($z) = @_;
1123         unless (ref $z) {
1124             my $t = $z + CORE::sqrt($z*$z + 1);
1125             return CORE::log($t) if $t;
1126         }
1127         my $t = &sqrt($z * $z + 1) + $z;
1128         # Try Taylor if looking bad (this usually means that
1129         # $z was large negative, therefore the sqrt is really
1130         # close to abs(z), summing that with z...)
1131         $t = 1/(2 * $z) - 1/(8 * $z**3) + 1/(16 * $z**5) - 5/(128 * $z**7)
1132             if $t == 0;
1133         return &log($t);
1134 }
1135
1136 #
1137 # atanh
1138 #
1139 # Computes the arc hyperbolic tangent atanh(z) = 1/2 log((1+z) / (1-z)).
1140 #
1141 sub atanh {
1142         my ($z) = @_;
1143         unless (ref $z) {
1144             return CORE::log((1 + $z)/(1 - $z))/2 if CORE::abs($z) < 1;
1145             $z = cplx($z, 0);
1146         }
1147         _divbyzero 'atanh(1)',  "1 - $z" if (1 - $z == 0);
1148         _logofzero 'atanh(-1)'           if (1 + $z == 0);
1149         return 0.5 * &log((1 + $z) / (1 - $z));
1150 }
1151
1152 #
1153 # asech
1154 #
1155 # Computes the hyperbolic arc secant asech(z) = acosh(1 / z).
1156 #
1157 sub asech {
1158         my ($z) = @_;
1159         _divbyzero 'asech(0)', "$z" if ($z == 0);
1160         return acosh(1 / $z);
1161 }
1162
1163 #
1164 # acsch
1165 #
1166 # Computes the hyperbolic arc cosecant acsch(z) = asinh(1 / z).
1167 #
1168 sub acsch {
1169         my ($z) = @_;
1170         _divbyzero 'acsch(0)', $z if ($z == 0);
1171         return asinh(1 / $z);
1172 }
1173
1174 #
1175 # acosech
1176 #
1177 # Alias for acosh().
1178 #
1179 sub acosech { Math::Complex::acsch(@_) }
1180
1181 #
1182 # acoth
1183 #
1184 # Computes the arc hyperbolic cotangent acoth(z) = 1/2 log((1+z) / (z-1)).
1185 #
1186 sub acoth {
1187         my ($z) = @_;
1188         _divbyzero 'acoth(0)'            if ($z == 0);
1189         unless (ref $z) {
1190             return CORE::log(($z + 1)/($z - 1))/2 if CORE::abs($z) > 1;
1191             $z = cplx($z, 0);
1192         }
1193         _divbyzero 'acoth(1)',  "$z - 1" if ($z - 1 == 0);
1194         _logofzero 'acoth(-1)', "1 + $z" if (1 + $z == 0);
1195         return &log((1 + $z) / ($z - 1)) / 2;
1196 }
1197
1198 #
1199 # acotanh
1200 #
1201 # Alias for acot().
1202 #
1203 sub acotanh { Math::Complex::acoth(@_) }
1204
1205 #
1206 # (atan2)
1207 #
1208 # Compute atan(z1/z2).
1209 #
1210 sub atan2 {
1211         my ($z1, $z2, $inverted) = @_;
1212         my ($re1, $im1, $re2, $im2);
1213         if ($inverted) {
1214             ($re1, $im1) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
1215             ($re2, $im2) = @{$z1->cartesian};
1216         } else {
1217             ($re1, $im1) = @{$z1->cartesian};
1218             ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
1219         }
1220         if ($im2 == 0) {
1221             return CORE::atan2($re1, $re2) if $im1 == 0;
1222             return ($im1<=>0) * pip2 if $re2 == 0;
1223         }
1224         my $w = atan($z1/$z2);
1225         my ($u, $v) = ref $w ? @{$w->cartesian} : ($w, 0);
1226         $u += pi   if $re2 < 0;
1227         $u -= pit2 if $u > pi;
1228         return cplx($u, $v);
1229 }
1230
1231 #
1232 # display_format
1233 # ->display_format
1234 #
1235 # Set (get if no argument) the display format for all complex numbers that
1236 # don't happen to have overridden it via ->display_format
1237 #
1238 # When called as an object method, this actually sets the display format for
1239 # the current object.
1240 #
1241 # Valid object formats are 'c' and 'p' for cartesian and polar. The first
1242 # letter is used actually, so the type can be fully spelled out for clarity.
1243 #
1244 sub display_format {
1245         my $self  = shift;
1246         my %display_format = %DISPLAY_FORMAT;
1247
1248         if (ref $self) {                        # Called as an object method
1249             if (exists $self->{display_format}) {
1250                 my %obj = %{$self->{display_format}};
1251                 @display_format{keys %obj} = values %obj;
1252             }
1253             if (@_ == 1) {
1254                 $display_format{style} = shift;
1255             } else {
1256                 my %new = @_;
1257                 @display_format{keys %new} = values %new;
1258             }
1259         } else {                                # Called as a class method
1260             if (@_ = 1) {
1261                 $display_format{style} = $self;
1262             } else {
1263                 my %new = @_;
1264                 @display_format{keys %new} = values %new;
1265             }
1266             undef $self;
1267         }
1268
1269         if (defined $self) {
1270             $self->{display_format} = { %display_format };
1271             return
1272                 wantarray ?
1273                     %{$self->{display_format}} :
1274                     $self->{display_format}->{style};
1275         }
1276
1277         %DISPLAY_FORMAT = %display_format;
1278         return
1279             wantarray ?
1280                 %DISPLAY_FORMAT :
1281                     $DISPLAY_FORMAT{style};
1282 }
1283
1284 #
1285 # (stringify)
1286 #
1287 # Show nicely formatted complex number under its cartesian or polar form,
1288 # depending on the current display format:
1289 #
1290 # . If a specific display format has been recorded for this object, use it.
1291 # . Otherwise, use the generic current default for all complex numbers,
1292 #   which is a package global variable.
1293 #
1294 sub stringify {
1295         my ($z) = shift;
1296
1297         my $style = $z->display_format;
1298
1299         $style = $DISPLAY_FORMAT{style} unless defined $style;
1300
1301         return $z->stringify_polar if $style =~ /^p/i;
1302         return $z->stringify_cartesian;
1303 }
1304
1305 #
1306 # ->stringify_cartesian
1307 #
1308 # Stringify as a cartesian representation 'a+bi'.
1309 #
1310 sub stringify_cartesian {
1311         my $z  = shift;
1312         my ($x, $y) = @{$z->cartesian};
1313         my ($re, $im);
1314
1315         my %format = $z->display_format;
1316         my $format = $format{format};
1317
1318         if ($x) {
1319             if ($x =~ /^NaN[QS]?$/i) {
1320                 $re = $x;
1321             } else {
1322                 if ($x =~ /^-?$Inf$/oi) {
1323                     $re = $x;
1324                 } else {
1325                     $re = defined $format ? sprintf($format, $x) : $x;
1326                 }
1327             }
1328         } else {
1329             undef $re;
1330         }
1331
1332         if ($y) {
1333             if ($y =~ /^(NaN[QS]?)$/i) {
1334                 $im = $y;
1335             } else {
1336                 if ($y =~ /^-?$Inf$/oi) {
1337                     $im = $y;
1338                 } else {
1339                     $im =
1340                         defined $format ?
1341                             sprintf($format, $y) :
1342                             ($y == 1 ? "" : ($y == -1 ? "-" : $y));
1343                 }
1344             }
1345             $im .= "i";
1346         } else {
1347             undef $im;
1348         }
1349
1350         my $str = $re;
1351
1352         if (defined $im) {
1353             if ($y < 0) {
1354                 $str .= $im;
1355             } elsif ($y > 0 || $im =~ /^NaN[QS]?i$/i)  {
1356                 $str .= "+" if defined $re;
1357                 $str .= $im;
1358             }
1359         } elsif (!defined $re) {
1360             $str = "0";
1361         }
1362
1363         return $str;
1364 }
1365
1366
1367 #
1368 # ->stringify_polar
1369 #
1370 # Stringify as a polar representation '[r,t]'.
1371 #
1372 sub stringify_polar {
1373         my $z  = shift;
1374         my ($r, $t) = @{$z->polar};
1375         my $theta;
1376
1377         my %format = $z->display_format;
1378         my $format = $format{format};
1379
1380         if ($t =~ /^NaN[QS]?$/i || $t =~ /^-?$Inf$/oi) {
1381             $theta = $t;
1382         } elsif ($t == pi) {
1383             $theta = "pi";
1384         } elsif ($r == 0 || $t == 0) {
1385             $theta = defined $format ? sprintf($format, $t) : $t;
1386         }
1387
1388         return "[$r,$theta]" if defined $theta;
1389
1390         #
1391         # Try to identify pi/n and friends.
1392         #
1393
1394         $t -= int(CORE::abs($t) / pit2) * pit2;
1395
1396         if ($format{polar_pretty_print}) {
1397             my ($a, $b);
1398             for $a (2..9) {
1399                 $b = $t * $a / pi;
1400                 if (int($b) == $b) {
1401                     $b = $b < 0 ? "-" : "" if CORE::abs($b) == 1;
1402                     $theta = "${b}pi/$a";
1403                     last;
1404                 }
1405             }
1406         }
1407
1408         if (defined $format) {
1409             $r     = sprintf($format, $r);
1410             $theta = sprintf($format, $theta) unless defined $theta;
1411         } else {
1412             $theta = $t unless defined $theta;
1413         }
1414
1415         return "[$r,$theta]";
1416 }
1417
1418 1;
1419 __END__
1420
1421 =head1 NAME
1422
1423 Math::Complex - complex numbers and associated mathematical functions
1424
1425 =head1 SYNOPSIS
1426
1427         use Math::Complex;
1428
1429         $z = Math::Complex->make(5, 6);
1430         $t = 4 - 3*i + $z;
1431         $j = cplxe(1, 2*pi/3);
1432
1433 =head1 DESCRIPTION
1434
1435 This package lets you create and manipulate complex numbers. By default,
1436 I<Perl> limits itself to real numbers, but an extra C<use> statement brings
1437 full complex support, along with a full set of mathematical functions
1438 typically associated with and/or extended to complex numbers.
1439
1440 If you wonder what complex numbers are, they were invented to be able to solve
1441 the following equation:
1442
1443         x*x = -1
1444
1445 and by definition, the solution is noted I<i> (engineers use I<j> instead since
1446 I<i> usually denotes an intensity, but the name does not matter). The number
1447 I<i> is a pure I<imaginary> number.
1448
1449 The arithmetics with pure imaginary numbers works just like you would expect
1450 it with real numbers... you just have to remember that
1451
1452         i*i = -1
1453
1454 so you have:
1455
1456         5i + 7i = i * (5 + 7) = 12i
1457         4i - 3i = i * (4 - 3) = i
1458         4i * 2i = -8
1459         6i / 2i = 3
1460         1 / i = -i
1461
1462 Complex numbers are numbers that have both a real part and an imaginary
1463 part, and are usually noted:
1464
1465         a + bi
1466
1467 where C<a> is the I<real> part and C<b> is the I<imaginary> part. The
1468 arithmetic with complex numbers is straightforward. You have to
1469 keep track of the real and the imaginary parts, but otherwise the
1470 rules used for real numbers just apply:
1471
1472         (4 + 3i) + (5 - 2i) = (4 + 5) + i(3 - 2) = 9 + i
1473         (2 + i) * (4 - i) = 2*4 + 4i -2i -i*i = 8 + 2i + 1 = 9 + 2i
1474
1475 A graphical representation of complex numbers is possible in a plane
1476 (also called the I<complex plane>, but it's really a 2D plane).
1477 The number
1478
1479         z = a + bi
1480
1481 is the point whose coordinates are (a, b). Actually, it would
1482 be the vector originating from (0, 0) to (a, b). It follows that the addition
1483 of two complex numbers is a vectorial addition.
1484
1485 Since there is a bijection between a point in the 2D plane and a complex
1486 number (i.e. the mapping is unique and reciprocal), a complex number
1487 can also be uniquely identified with polar coordinates:
1488
1489         [rho, theta]
1490
1491 where C<rho> is the distance to the origin, and C<theta> the angle between
1492 the vector and the I<x> axis. There is a notation for this using the
1493 exponential form, which is:
1494
1495         rho * exp(i * theta)
1496
1497 where I<i> is the famous imaginary number introduced above. Conversion
1498 between this form and the cartesian form C<a + bi> is immediate:
1499
1500         a = rho * cos(theta)
1501         b = rho * sin(theta)
1502
1503 which is also expressed by this formula:
1504
1505         z = rho * exp(i * theta) = rho * (cos theta + i * sin theta)
1506
1507 In other words, it's the projection of the vector onto the I<x> and I<y>
1508 axes. Mathematicians call I<rho> the I<norm> or I<modulus> and I<theta>
1509 the I<argument> of the complex number. The I<norm> of C<z> will be
1510 noted C<abs(z)>.
1511
1512 The polar notation (also known as the trigonometric
1513 representation) is much more handy for performing multiplications and
1514 divisions of complex numbers, whilst the cartesian notation is better
1515 suited for additions and subtractions. Real numbers are on the I<x>
1516 axis, and therefore I<theta> is zero or I<pi>.
1517
1518 All the common operations that can be performed on a real number have
1519 been defined to work on complex numbers as well, and are merely
1520 I<extensions> of the operations defined on real numbers. This means
1521 they keep their natural meaning when there is no imaginary part, provided
1522 the number is within their definition set.
1523
1524 For instance, the C<sqrt> routine which computes the square root of
1525 its argument is only defined for non-negative real numbers and yields a
1526 non-negative real number (it is an application from B<R+> to B<R+>).
1527 If we allow it to return a complex number, then it can be extended to
1528 negative real numbers to become an application from B<R> to B<C> (the
1529 set of complex numbers):
1530
1531         sqrt(x) = x >= 0 ? sqrt(x) : sqrt(-x)*i
1532
1533 It can also be extended to be an application from B<C> to B<C>,
1534 whilst its restriction to B<R> behaves as defined above by using
1535 the following definition:
1536
1537         sqrt(z = [r,t]) = sqrt(r) * exp(i * t/2)
1538
1539 Indeed, a negative real number can be noted C<[x,pi]> (the modulus
1540 I<x> is always non-negative, so C<[x,pi]> is really C<-x>, a negative
1541 number) and the above definition states that
1542
1543         sqrt([x,pi]) = sqrt(x) * exp(i*pi/2) = [sqrt(x),pi/2] = sqrt(x)*i
1544
1545 which is exactly what we had defined for negative real numbers above.
1546 The C<sqrt> returns only one of the solutions: if you want the both,
1547 use the C<root> function.
1548
1549 All the common mathematical functions defined on real numbers that
1550 are extended to complex numbers share that same property of working
1551 I<as usual> when the imaginary part is zero (otherwise, it would not
1552 be called an extension, would it?).
1553
1554 A I<new> operation possible on a complex number that is
1555 the identity for real numbers is called the I<conjugate>, and is noted
1556 with an horizontal bar above the number, or C<~z> here.
1557
1558          z = a + bi
1559         ~z = a - bi
1560
1561 Simple... Now look:
1562
1563         z * ~z = (a + bi) * (a - bi) = a*a + b*b
1564
1565 We saw that the norm of C<z> was noted C<abs(z)> and was defined as the
1566 distance to the origin, also known as:
1567
1568         rho = abs(z) = sqrt(a*a + b*b)
1569
1570 so
1571
1572         z * ~z = abs(z) ** 2
1573
1574 If z is a pure real number (i.e. C<b == 0>), then the above yields:
1575
1576         a * a = abs(a) ** 2
1577
1578 which is true (C<abs> has the regular meaning for real number, i.e. stands
1579 for the absolute value). This example explains why the norm of C<z> is
1580 noted C<abs(z)>: it extends the C<abs> function to complex numbers, yet
1581 is the regular C<abs> we know when the complex number actually has no
1582 imaginary part... This justifies I<a posteriori> our use of the C<abs>
1583 notation for the norm.
1584
1585 =head1 OPERATIONS
1586
1587 Given the following notations:
1588
1589         z1 = a + bi = r1 * exp(i * t1)
1590         z2 = c + di = r2 * exp(i * t2)
1591         z = <any complex or real number>
1592
1593 the following (overloaded) operations are supported on complex numbers:
1594
1595         z1 + z2 = (a + c) + i(b + d)
1596         z1 - z2 = (a - c) + i(b - d)
1597         z1 * z2 = (r1 * r2) * exp(i * (t1 + t2))
1598         z1 / z2 = (r1 / r2) * exp(i * (t1 - t2))
1599         z1 ** z2 = exp(z2 * log z1)
1600         ~z = a - bi
1601         abs(z) = r1 = sqrt(a*a + b*b)
1602         sqrt(z) = sqrt(r1) * exp(i * t/2)
1603         exp(z) = exp(a) * exp(i * b)
1604         log(z) = log(r1) + i*t
1605         sin(z) = 1/2i (exp(i * z1) - exp(-i * z))
1606         cos(z) = 1/2 (exp(i * z1) + exp(-i * z))
1607         atan2(z1, z2) = atan(z1/z2)
1608
1609 The following extra operations are supported on both real and complex
1610 numbers:
1611
1612         Re(z) = a
1613         Im(z) = b
1614         arg(z) = t
1615         abs(z) = r
1616
1617         cbrt(z) = z ** (1/3)
1618         log10(z) = log(z) / log(10)
1619         logn(z, n) = log(z) / log(n)
1620
1621         tan(z) = sin(z) / cos(z)
1622
1623         csc(z) = 1 / sin(z)
1624         sec(z) = 1 / cos(z)
1625         cot(z) = 1 / tan(z)
1626
1627         asin(z) = -i * log(i*z + sqrt(1-z*z))
1628         acos(z) = -i * log(z + i*sqrt(1-z*z))
1629         atan(z) = i/2 * log((i+z) / (i-z))
1630
1631         acsc(z) = asin(1 / z)
1632         asec(z) = acos(1 / z)
1633         acot(z) = atan(1 / z) = -i/2 * log((i+z) / (z-i))
1634
1635         sinh(z) = 1/2 (exp(z) - exp(-z))
1636         cosh(z) = 1/2 (exp(z) + exp(-z))
1637         tanh(z) = sinh(z) / cosh(z) = (exp(z) - exp(-z)) / (exp(z) + exp(-z))
1638
1639         csch(z) = 1 / sinh(z)
1640         sech(z) = 1 / cosh(z)
1641         coth(z) = 1 / tanh(z)
1642
1643         asinh(z) = log(z + sqrt(z*z+1))
1644         acosh(z) = log(z + sqrt(z*z-1))
1645         atanh(z) = 1/2 * log((1+z) / (1-z))
1646
1647         acsch(z) = asinh(1 / z)
1648         asech(z) = acosh(1 / z)
1649         acoth(z) = atanh(1 / z) = 1/2 * log((1+z) / (z-1))
1650
1651 I<arg>, I<abs>, I<log>, I<csc>, I<cot>, I<acsc>, I<acot>, I<csch>,
1652 I<coth>, I<acosech>, I<acotanh>, have aliases I<rho>, I<theta>, I<ln>,
1653 I<cosec>, I<cotan>, I<acosec>, I<acotan>, I<cosech>, I<cotanh>,
1654 I<acosech>, I<acotanh>, respectively.  C<Re>, C<Im>, C<arg>, C<abs>,
1655 C<rho>, and C<theta> can be used also also mutators.  The C<cbrt>
1656 returns only one of the solutions: if you want all three, use the
1657 C<root> function.
1658
1659 The I<root> function is available to compute all the I<n>
1660 roots of some complex, where I<n> is a strictly positive integer.
1661 There are exactly I<n> such roots, returned as a list. Getting the
1662 number mathematicians call C<j> such that:
1663
1664         1 + j + j*j = 0;
1665
1666 is a simple matter of writing:
1667
1668         $j = ((root(1, 3))[1];
1669
1670 The I<k>th root for C<z = [r,t]> is given by:
1671
1672         (root(z, n))[k] = r**(1/n) * exp(i * (t + 2*k*pi)/n)
1673
1674 The I<spaceship> comparison operator, E<lt>=E<gt>, is also defined. In
1675 order to ensure its restriction to real numbers is conform to what you
1676 would expect, the comparison is run on the real part of the complex
1677 number first, and imaginary parts are compared only when the real
1678 parts match.
1679
1680 =head1 CREATION
1681
1682 To create a complex number, use either:
1683
1684         $z = Math::Complex->make(3, 4);
1685         $z = cplx(3, 4);
1686
1687 if you know the cartesian form of the number, or
1688
1689         $z = 3 + 4*i;
1690
1691 if you like. To create a number using the polar form, use either:
1692
1693         $z = Math::Complex->emake(5, pi/3);
1694         $x = cplxe(5, pi/3);
1695
1696 instead. The first argument is the modulus, the second is the angle
1697 (in radians, the full circle is 2*pi).  (Mnemonic: C<e> is used as a
1698 notation for complex numbers in the polar form).
1699
1700 It is possible to write:
1701
1702         $x = cplxe(-3, pi/4);
1703
1704 but that will be silently converted into C<[3,-3pi/4]>, since the
1705 modulus must be non-negative (it represents the distance to the origin
1706 in the complex plane).
1707
1708 It is also possible to have a complex number as either argument of
1709 either the C<make> or C<emake>: the appropriate component of
1710 the argument will be used.
1711
1712         $z1 = cplx(-2,  1);
1713         $z2 = cplx($z1, 4);
1714
1715 =head1 STRINGIFICATION
1716
1717 When printed, a complex number is usually shown under its cartesian
1718 style I<a+bi>, but there are legitimate cases where the polar style
1719 I<[r,t]> is more appropriate.
1720
1721 By calling the class method C<Math::Complex::display_format> and
1722 supplying either C<"polar"> or C<"cartesian"> as an argument, you
1723 override the default display style, which is C<"cartesian">. Not
1724 supplying any argument returns the current settings.
1725
1726 This default can be overridden on a per-number basis by calling the
1727 C<display_format> method instead. As before, not supplying any argument
1728 returns the current display style for this number. Otherwise whatever you
1729 specify will be the new display style for I<this> particular number.
1730
1731 For instance:
1732
1733         use Math::Complex;
1734
1735         Math::Complex::display_format('polar');
1736         $j = (root(1, 3))[1];
1737         print "j = $j\n";               # Prints "j = [1,2pi/3]"
1738         $j->display_format('cartesian');
1739         print "j = $j\n";               # Prints "j = -0.5+0.866025403784439i"
1740
1741 The polar style attempts to emphasize arguments like I<k*pi/n>
1742 (where I<n> is a positive integer and I<k> an integer within [-9, +9]),
1743 this is called I<polar pretty-printing>.
1744
1745 =head2 CHANGED IN PERL 5.6
1746
1747 The C<display_format> class method and the corresponding
1748 C<display_format> object method can now be called using
1749 a parameter hash instead of just a one parameter.
1750
1751 The old display format style, which can have values C<"cartesian"> or
1752 C<"polar">, can be changed using the C<"style"> parameter.
1753
1754         $j->display_format(style => "polar");
1755
1756 The one parameter calling convention also still works.
1757
1758         $j->display_format("polar");
1759
1760 There are two new display parameters.
1761
1762 The first one is C<"format">, which is a sprintf()-style format string
1763 to be used for both numeric parts of the complex number(s).  The is
1764 somewhat system-dependent but most often it corresponds to C<"%.15g">.
1765 You can revert to the default by setting the C<format> to C<undef>.
1766
1767         # the $j from the above example
1768
1769         $j->display_format('format' => '%.5f');
1770         print "j = $j\n";               # Prints "j = -0.50000+0.86603i"
1771         $j->display_format('format' => undef);
1772         print "j = $j\n";               # Prints "j = -0.5+0.86603i"
1773
1774 Notice that this affects also the return values of the
1775 C<display_format> methods: in list context the whole parameter hash
1776 will be returned, as opposed to only the style parameter value.
1777 This is a potential incompatibility with earlier versions if you
1778 have been calling the C<display_format> method in list context.
1779
1780 The second new display parameter is C<"polar_pretty_print">, which can
1781 be set to true or false, the default being true.  See the previous
1782 section for what this means.
1783
1784 =head1 USAGE
1785
1786 Thanks to overloading, the handling of arithmetics with complex numbers
1787 is simple and almost transparent.
1788
1789 Here are some examples:
1790
1791         use Math::Complex;
1792
1793         $j = cplxe(1, 2*pi/3);  # $j ** 3 == 1
1794         print "j = $j, j**3 = ", $j ** 3, "\n";
1795         print "1 + j + j**2 = ", 1 + $j + $j**2, "\n";
1796
1797         $z = -16 + 0*i;                 # Force it to be a complex
1798         print "sqrt($z) = ", sqrt($z), "\n";
1799
1800         $k = exp(i * 2*pi/3);
1801         print "$j - $k = ", $j - $k, "\n";
1802
1803         $z->Re(3);                      # Re, Im, arg, abs,
1804         $j->arg(2);                     # (the last two aka rho, theta)
1805                                         # can be used also as mutators.
1806
1807 =head1 ERRORS DUE TO DIVISION BY ZERO OR LOGARITHM OF ZERO
1808
1809 The division (/) and the following functions
1810
1811         log     ln      log10   logn
1812         tan     sec     csc     cot
1813         atan    asec    acsc    acot
1814         tanh    sech    csch    coth
1815         atanh   asech   acsch   acoth
1816
1817 cannot be computed for all arguments because that would mean dividing
1818 by zero or taking logarithm of zero. These situations cause fatal
1819 runtime errors looking like this
1820
1821         cot(0): Division by zero.
1822         (Because in the definition of cot(0), the divisor sin(0) is 0)
1823         Died at ...
1824
1825 or
1826
1827         atanh(-1): Logarithm of zero.
1828         Died at...
1829
1830 For the C<csc>, C<cot>, C<asec>, C<acsc>, C<acot>, C<csch>, C<coth>,
1831 C<asech>, C<acsch>, the argument cannot be C<0> (zero).  For the the
1832 logarithmic functions and the C<atanh>, C<acoth>, the argument cannot
1833 be C<1> (one).  For the C<atanh>, C<acoth>, the argument cannot be
1834 C<-1> (minus one).  For the C<atan>, C<acot>, the argument cannot be
1835 C<i> (the imaginary unit).  For the C<atan>, C<acoth>, the argument
1836 cannot be C<-i> (the negative imaginary unit).  For the C<tan>,
1837 C<sec>, C<tanh>, the argument cannot be I<pi/2 + k * pi>, where I<k>
1838 is any integer.
1839
1840 Note that because we are operating on approximations of real numbers,
1841 these errors can happen when merely `too close' to the singularities
1842 listed above.
1843
1844 =head1 ERRORS DUE TO INDIGESTIBLE ARGUMENTS
1845
1846 The C<make> and C<emake> accept both real and complex arguments.
1847 When they cannot recognize the arguments they will die with error
1848 messages like the following
1849
1850     Math::Complex::make: Cannot take real part of ...
1851     Math::Complex::make: Cannot take real part of ...
1852     Math::Complex::emake: Cannot take rho of ...
1853     Math::Complex::emake: Cannot take theta of ...
1854
1855 =head1 BUGS
1856
1857 Saying C<use Math::Complex;> exports many mathematical routines in the
1858 caller environment and even overrides some (C<sqrt>, C<log>).
1859 This is construed as a feature by the Authors, actually... ;-)
1860
1861 All routines expect to be given real or complex numbers. Don't attempt to
1862 use BigFloat, since Perl has currently no rule to disambiguate a '+'
1863 operation (for instance) between two overloaded entities.
1864
1865 In Cray UNICOS there is some strange numerical instability that results
1866 in root(), cos(), sin(), cosh(), sinh(), losing accuracy fast.  Beware.
1867 The bug may be in UNICOS math libs, in UNICOS C compiler, in Math::Complex.
1868 Whatever it is, it does not manifest itself anywhere else where Perl runs.
1869
1870 =head1 AUTHORS
1871
1872 Raphael Manfredi <F<Raphael_Manfredi@pobox.com>> and
1873 Jarkko Hietaniemi <F<jhi@iki.fi>>.
1874
1875 Extensive patches by Daniel S. Lewart <F<d-lewart@uiuc.edu>>.
1876
1877 =cut
1878
1879 1;
1880
1881 # eof