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