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