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 =pod
1409 =head1 NAME
1410
1411 Math::Complex - complex numbers and associated mathematical functions
1412
1413 =head1 SYNOPSIS
1414
1415         use Math::Complex;
1416
1417         $z = Math::Complex->make(5, 6);
1418         $t = 4 - 3*i + $z;
1419         $j = cplxe(1, 2*pi/3);
1420
1421 =head1 DESCRIPTION
1422
1423 This package lets you create and manipulate complex numbers. By default,
1424 I<Perl> limits itself to real numbers, but an extra C<use> statement brings
1425 full complex support, along with a full set of mathematical functions
1426 typically associated with and/or extended to complex numbers.
1427
1428 If you wonder what complex numbers are, they were invented to be able to solve
1429 the following equation:
1430
1431         x*x = -1
1432
1433 and by definition, the solution is noted I<i> (engineers use I<j> instead since
1434 I<i> usually denotes an intensity, but the name does not matter). The number
1435 I<i> is a pure I<imaginary> number.
1436
1437 The arithmetics with pure imaginary numbers works just like you would expect
1438 it with real numbers... you just have to remember that
1439
1440         i*i = -1
1441
1442 so you have:
1443
1444         5i + 7i = i * (5 + 7) = 12i
1445         4i - 3i = i * (4 - 3) = i
1446         4i * 2i = -8
1447         6i / 2i = 3
1448         1 / i = -i
1449
1450 Complex numbers are numbers that have both a real part and an imaginary
1451 part, and are usually noted:
1452
1453         a + bi
1454
1455 where C<a> is the I<real> part and C<b> is the I<imaginary> part. The
1456 arithmetic with complex numbers is straightforward. You have to
1457 keep track of the real and the imaginary parts, but otherwise the
1458 rules used for real numbers just apply:
1459
1460         (4 + 3i) + (5 - 2i) = (4 + 5) + i(3 - 2) = 9 + i
1461         (2 + i) * (4 - i) = 2*4 + 4i -2i -i*i = 8 + 2i + 1 = 9 + 2i
1462
1463 A graphical representation of complex numbers is possible in a plane
1464 (also called the I<complex plane>, but it's really a 2D plane).
1465 The number
1466
1467         z = a + bi
1468
1469 is the point whose coordinates are (a, b). Actually, it would
1470 be the vector originating from (0, 0) to (a, b). It follows that the addition
1471 of two complex numbers is a vectorial addition.
1472
1473 Since there is a bijection between a point in the 2D plane and a complex
1474 number (i.e. the mapping is unique and reciprocal), a complex number
1475 can also be uniquely identified with polar coordinates:
1476
1477         [rho, theta]
1478
1479 where C<rho> is the distance to the origin, and C<theta> the angle between
1480 the vector and the I<x> axis. There is a notation for this using the
1481 exponential form, which is:
1482
1483         rho * exp(i * theta)
1484
1485 where I<i> is the famous imaginary number introduced above. Conversion
1486 between this form and the cartesian form C<a + bi> is immediate:
1487
1488         a = rho * cos(theta)
1489         b = rho * sin(theta)
1490
1491 which is also expressed by this formula:
1492
1493         z = rho * exp(i * theta) = rho * (cos theta + i * sin theta)
1494
1495 In other words, it's the projection of the vector onto the I<x> and I<y>
1496 axes. Mathematicians call I<rho> the I<norm> or I<modulus> and I<theta>
1497 the I<argument> of the complex number. The I<norm> of C<z> will be
1498 noted C<abs(z)>.
1499
1500 The polar notation (also known as the trigonometric
1501 representation) is much more handy for performing multiplications and
1502 divisions of complex numbers, whilst the cartesian notation is better
1503 suited for additions and subtractions. Real numbers are on the I<x>
1504 axis, and therefore I<theta> is zero or I<pi>.
1505
1506 All the common operations that can be performed on a real number have
1507 been defined to work on complex numbers as well, and are merely
1508 I<extensions> of the operations defined on real numbers. This means
1509 they keep their natural meaning when there is no imaginary part, provided
1510 the number is within their definition set.
1511
1512 For instance, the C<sqrt> routine which computes the square root of
1513 its argument is only defined for non-negative real numbers and yields a
1514 non-negative real number (it is an application from B<R+> to B<R+>).
1515 If we allow it to return a complex number, then it can be extended to
1516 negative real numbers to become an application from B<R> to B<C> (the
1517 set of complex numbers):
1518
1519         sqrt(x) = x >= 0 ? sqrt(x) : sqrt(-x)*i
1520
1521 It can also be extended to be an application from B<C> to B<C>,
1522 whilst its restriction to B<R> behaves as defined above by using
1523 the following definition:
1524
1525         sqrt(z = [r,t]) = sqrt(r) * exp(i * t/2)
1526
1527 Indeed, a negative real number can be noted C<[x,pi]> (the modulus
1528 I<x> is always non-negative, so C<[x,pi]> is really C<-x>, a negative
1529 number) and the above definition states that
1530
1531         sqrt([x,pi]) = sqrt(x) * exp(i*pi/2) = [sqrt(x),pi/2] = sqrt(x)*i
1532
1533 which is exactly what we had defined for negative real numbers above.
1534 The C<sqrt> returns only one of the solutions: if you want the both,
1535 use the C<root> function.
1536
1537 All the common mathematical functions defined on real numbers that
1538 are extended to complex numbers share that same property of working
1539 I<as usual> when the imaginary part is zero (otherwise, it would not
1540 be called an extension, would it?).
1541
1542 A I<new> operation possible on a complex number that is
1543 the identity for real numbers is called the I<conjugate>, and is noted
1544 with an horizontal bar above the number, or C<~z> here.
1545
1546          z = a + bi
1547         ~z = a - bi
1548
1549 Simple... Now look:
1550
1551         z * ~z = (a + bi) * (a - bi) = a*a + b*b
1552
1553 We saw that the norm of C<z> was noted C<abs(z)> and was defined as the
1554 distance to the origin, also known as:
1555
1556         rho = abs(z) = sqrt(a*a + b*b)
1557
1558 so
1559
1560         z * ~z = abs(z) ** 2
1561
1562 If z is a pure real number (i.e. C<b == 0>), then the above yields:
1563
1564         a * a = abs(a) ** 2
1565
1566 which is true (C<abs> has the regular meaning for real number, i.e. stands
1567 for the absolute value). This example explains why the norm of C<z> is
1568 noted C<abs(z)>: it extends the C<abs> function to complex numbers, yet
1569 is the regular C<abs> we know when the complex number actually has no
1570 imaginary part... This justifies I<a posteriori> our use of the C<abs>
1571 notation for the norm.
1572
1573 =head1 OPERATIONS
1574
1575 Given the following notations:
1576
1577         z1 = a + bi = r1 * exp(i * t1)
1578         z2 = c + di = r2 * exp(i * t2)
1579         z = <any complex or real number>
1580
1581 the following (overloaded) operations are supported on complex numbers:
1582
1583         z1 + z2 = (a + c) + i(b + d)
1584         z1 - z2 = (a - c) + i(b - d)
1585         z1 * z2 = (r1 * r2) * exp(i * (t1 + t2))
1586         z1 / z2 = (r1 / r2) * exp(i * (t1 - t2))
1587         z1 ** z2 = exp(z2 * log z1)
1588         ~z = a - bi
1589         abs(z) = r1 = sqrt(a*a + b*b)
1590         sqrt(z) = sqrt(r1) * exp(i * t/2)
1591         exp(z) = exp(a) * exp(i * b)
1592         log(z) = log(r1) + i*t
1593         sin(z) = 1/2i (exp(i * z1) - exp(-i * z))
1594         cos(z) = 1/2 (exp(i * z1) + exp(-i * z))
1595         atan2(z1, z2) = atan(z1/z2)
1596
1597 The following extra operations are supported on both real and complex
1598 numbers:
1599
1600         Re(z) = a
1601         Im(z) = b
1602         arg(z) = t
1603         abs(z) = r
1604
1605         cbrt(z) = z ** (1/3)
1606         log10(z) = log(z) / log(10)
1607         logn(z, n) = log(z) / log(n)
1608
1609         tan(z) = sin(z) / cos(z)
1610
1611         csc(z) = 1 / sin(z)
1612         sec(z) = 1 / cos(z)
1613         cot(z) = 1 / tan(z)
1614
1615         asin(z) = -i * log(i*z + sqrt(1-z*z))
1616         acos(z) = -i * log(z + i*sqrt(1-z*z))
1617         atan(z) = i/2 * log((i+z) / (i-z))
1618
1619         acsc(z) = asin(1 / z)
1620         asec(z) = acos(1 / z)
1621         acot(z) = atan(1 / z) = -i/2 * log((i+z) / (z-i))
1622
1623         sinh(z) = 1/2 (exp(z) - exp(-z))
1624         cosh(z) = 1/2 (exp(z) + exp(-z))
1625         tanh(z) = sinh(z) / cosh(z) = (exp(z) - exp(-z)) / (exp(z) + exp(-z))
1626
1627         csch(z) = 1 / sinh(z)
1628         sech(z) = 1 / cosh(z)
1629         coth(z) = 1 / tanh(z)
1630
1631         asinh(z) = log(z + sqrt(z*z+1))
1632         acosh(z) = log(z + sqrt(z*z-1))
1633         atanh(z) = 1/2 * log((1+z) / (1-z))
1634
1635         acsch(z) = asinh(1 / z)
1636         asech(z) = acosh(1 / z)
1637         acoth(z) = atanh(1 / z) = 1/2 * log((1+z) / (z-1))
1638
1639 I<arg>, I<abs>, I<log>, I<csc>, I<cot>, I<acsc>, I<acot>, I<csch>,
1640 I<coth>, I<acosech>, I<acotanh>, have aliases I<rho>, I<theta>, I<ln>,
1641 I<cosec>, I<cotan>, I<acosec>, I<acotan>, I<cosech>, I<cotanh>,
1642 I<acosech>, I<acotanh>, respectively.  C<Re>, C<Im>, C<arg>, C<abs>,
1643 C<rho>, and C<theta> can be used also also mutators.  The C<cbrt>
1644 returns only one of the solutions: if you want all three, use the
1645 C<root> function.
1646
1647 The I<root> function is available to compute all the I<n>
1648 roots of some complex, where I<n> is a strictly positive integer.
1649 There are exactly I<n> such roots, returned as a list. Getting the
1650 number mathematicians call C<j> such that:
1651
1652         1 + j + j*j = 0;
1653
1654 is a simple matter of writing:
1655
1656         $j = ((root(1, 3))[1];
1657
1658 The I<k>th root for C<z = [r,t]> is given by:
1659
1660         (root(z, n))[k] = r**(1/n) * exp(i * (t + 2*k*pi)/n)
1661
1662 The I<spaceship> comparison operator, E<lt>=E<gt>, is also defined. In
1663 order to ensure its restriction to real numbers is conform to what you
1664 would expect, the comparison is run on the real part of the complex
1665 number first, and imaginary parts are compared only when the real
1666 parts match.
1667
1668 =head1 CREATION
1669
1670 To create a complex number, use either:
1671
1672         $z = Math::Complex->make(3, 4);
1673         $z = cplx(3, 4);
1674
1675 if you know the cartesian form of the number, or
1676
1677         $z = 3 + 4*i;
1678
1679 if you like. To create a number using the polar form, use either:
1680
1681         $z = Math::Complex->emake(5, pi/3);
1682         $x = cplxe(5, pi/3);
1683
1684 instead. The first argument is the modulus, the second is the angle
1685 (in radians, the full circle is 2*pi).  (Mnemonic: C<e> is used as a
1686 notation for complex numbers in the polar form).
1687
1688 It is possible to write:
1689
1690         $x = cplxe(-3, pi/4);
1691
1692 but that will be silently converted into C<[3,-3pi/4]>, since the
1693 modulus must be non-negative (it represents the distance to the origin
1694 in the complex plane).
1695
1696 It is also possible to have a complex number as either argument of
1697 either the C<make> or C<emake>: the appropriate component of
1698 the argument will be used.
1699
1700         $z1 = cplx(-2,  1);
1701         $z2 = cplx($z1, 4);
1702
1703 =head1 STRINGIFICATION
1704
1705 When printed, a complex number is usually shown under its cartesian
1706 style I<a+bi>, but there are legitimate cases where the polar style
1707 I<[r,t]> is more appropriate.
1708
1709 By calling the class method C<Math::Complex::display_format> and
1710 supplying either C<"polar"> or C<"cartesian"> as an argument, you
1711 override the default display style, which is C<"cartesian">. Not
1712 supplying any argument returns the current settings.
1713
1714 This default can be overridden on a per-number basis by calling the
1715 C<display_format> method instead. As before, not supplying any argument
1716 returns the current display style for this number. Otherwise whatever you
1717 specify will be the new display style for I<this> particular number.
1718
1719 For instance:
1720
1721         use Math::Complex;
1722
1723         Math::Complex::display_format('polar');
1724         $j = (root(1, 3))[1];
1725         print "j = $j\n";               # Prints "j = [1,2pi/3]"
1726         $j->display_format('cartesian');
1727         print "j = $j\n";               # Prints "j = -0.5+0.866025403784439i"
1728
1729 The polar style attempts to emphasize arguments like I<k*pi/n>
1730 (where I<n> is a positive integer and I<k> an integer within [-9,+9]),
1731 this is called I<polar pretty-printing>.
1732
1733 =head2 CHANGED IN PERL 5.6
1734
1735 The C<display_format> class method and the corresponding
1736 C<display_format> object method can now be called using
1737 a parameter hash instead of just a one parameter.
1738
1739 The old display format style, which can have values C<"cartesian"> or
1740 C<"polar">, can be changed using the C<"style"> parameter.  (The one
1741 parameter calling convention also still works.)
1742
1743 There are two new display parameters.
1744
1745 The first one is C<"format">, which is a sprintf()-style format
1746 string to be used for both parts of the complex number(s).  The
1747 default is C<undef>, which corresponds usually (this is somewhat
1748 system-dependent) to C<"%.15g">.  You can revert to the default by
1749 setting the format string to C<undef>.
1750
1751         # the $j from the above example
1752
1753         $j->display_format('format' => '%.5f');
1754         print "j = $j\n";               # Prints "j = -0.50000+0.86603i"
1755         $j->display_format('format' => '%.6f');
1756         print "j = $j\n";               # Prints "j = -0.5+0.86603i"
1757
1758 Notice that this affects also the return values of the
1759 C<display_format> methods: in list context the whole parameter hash
1760 will be returned, as opposed to only the style parameter value.  If
1761 you want to know the whole truth for a complex number, you must call
1762 both the class method and the object method:
1763
1764 The second new display parameter is C<"polar_pretty_print">, which can
1765 be set to true or false, the default being true.  See the previous
1766 section for what this means.
1767
1768 =head1 USAGE
1769
1770 Thanks to overloading, the handling of arithmetics with complex numbers
1771 is simple and almost transparent.
1772
1773 Here are some examples:
1774
1775         use Math::Complex;
1776
1777         $j = cplxe(1, 2*pi/3);  # $j ** 3 == 1
1778         print "j = $j, j**3 = ", $j ** 3, "\n";
1779         print "1 + j + j**2 = ", 1 + $j + $j**2, "\n";
1780
1781         $z = -16 + 0*i;                 # Force it to be a complex
1782         print "sqrt($z) = ", sqrt($z), "\n";
1783
1784         $k = exp(i * 2*pi/3);
1785         print "$j - $k = ", $j - $k, "\n";
1786
1787         $z->Re(3);                      # Re, Im, arg, abs,
1788         $j->arg(2);                     # (the last two aka rho, theta)
1789                                         # can be used also as mutators.
1790
1791 =head1 ERRORS DUE TO DIVISION BY ZERO OR LOGARITHM OF ZERO
1792
1793 The division (/) and the following functions
1794
1795         log     ln      log10   logn
1796         tan     sec     csc     cot
1797         atan    asec    acsc    acot
1798         tanh    sech    csch    coth
1799         atanh   asech   acsch   acoth
1800
1801 cannot be computed for all arguments because that would mean dividing
1802 by zero or taking logarithm of zero. These situations cause fatal
1803 runtime errors looking like this
1804
1805         cot(0): Division by zero.
1806         (Because in the definition of cot(0), the divisor sin(0) is 0)
1807         Died at ...
1808
1809 or
1810
1811         atanh(-1): Logarithm of zero.
1812         Died at...
1813
1814 For the C<csc>, C<cot>, C<asec>, C<acsc>, C<acot>, C<csch>, C<coth>,
1815 C<asech>, C<acsch>, the argument cannot be C<0> (zero).  For the the
1816 logarithmic functions and the C<atanh>, C<acoth>, the argument cannot
1817 be C<1> (one).  For the C<atanh>, C<acoth>, the argument cannot be
1818 C<-1> (minus one).  For the C<atan>, C<acot>, the argument cannot be
1819 C<i> (the imaginary unit).  For the C<atan>, C<acoth>, the argument
1820 cannot be C<-i> (the negative imaginary unit).  For the C<tan>,
1821 C<sec>, C<tanh>, the argument cannot be I<pi/2 + k * pi>, where I<k>
1822 is any integer.
1823
1824 Note that because we are operating on approximations of real numbers,
1825 these errors can happen when merely `too close' to the singularities
1826 listed above.  For example C<tan(2*atan2(1,1)+1e-15)> will die of
1827 division by zero.
1828
1829 =head1 ERRORS DUE TO INDIGESTIBLE ARGUMENTS
1830
1831 The C<make> and C<emake> accept both real and complex arguments.
1832 When they cannot recognize the arguments they will die with error
1833 messages like the following
1834
1835     Math::Complex::make: Cannot take real part of ...
1836     Math::Complex::make: Cannot take real part of ...
1837     Math::Complex::emake: Cannot take rho of ...
1838     Math::Complex::emake: Cannot take theta of ...
1839
1840 =head1 BUGS
1841
1842 Saying C<use Math::Complex;> exports many mathematical routines in the
1843 caller environment and even overrides some (C<sqrt>, C<log>).
1844 This is construed as a feature by the Authors, actually... ;-)
1845
1846 All routines expect to be given real or complex numbers. Don't attempt to
1847 use BigFloat, since Perl has currently no rule to disambiguate a '+'
1848 operation (for instance) between two overloaded entities.
1849
1850 In Cray UNICOS there is some strange numerical instability that results
1851 in root(), cos(), sin(), cosh(), sinh(), losing accuracy fast.  Beware.
1852 The bug may be in UNICOS math libs, in UNICOS C compiler, in Math::Complex.
1853 Whatever it is, it does not manifest itself anywhere else where Perl runs.
1854
1855 =head1 AUTHORS
1856
1857 Raphael Manfredi <F<Raphael_Manfredi@pobox.com>> and
1858 Jarkko Hietaniemi <F<jhi@iki.fi>>.
1859
1860 Extensive patches by Daniel S. Lewart <F<d-lewart@uiuc.edu>>.
1861
1862 =cut
1863
1864 1;
1865
1866 # eof