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