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