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