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