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