lib/Math/Complex.pm

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