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