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