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