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