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