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