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