lib/Math/Complex.pm

   1 # $RCSFile$
   2 #
   3 # Complex numbers and associated mathematical functions
   4 # -- Raphael Manfredi, September 1996
   5 # -- Jarkko Hietaniemi, March 1997
   6
   7 require Exporter;
   8 package Math::Complex;
   9
  10 use strict;
  11
  12 use vars qw($VERSION @ISA
  13             @EXPORT %EXPORT_TAGS
  14             $package $display
  15             $pi $i $ilog10 $logn %logn);
  16
  17 @ISA = qw(Exporter);
  18
  19 $VERSION = 1.01;
  20
  21 my @trig = qw(
  22               pi
  23               tan
  24               csc cosec sec cot cotan
  25               asin acos atan
  26               acsc acosec asec acot acotan
  27               sinh cosh tanh
  28               csch cosech sech coth cotanh
  29               asinh acosh atanh
  30               acsch acosech asech acoth acotanh
  31              );
  32
  33 @EXPORT = (qw(
  34              i Re Im arg
  35              sqrt exp log ln
  36              log10 logn cbrt root
  37              cplx cplxe
  38              ),
  39            @trig);
  40
  41 %EXPORT_TAGS = (
  42     'trig' => [@trig],
  43 );
  44
  45 use overload
  46         '+'     => \&plus,
  47         '-'     => \&minus,
  48         '*'     => \&multiply,
  49         '/'     => \&divide,
  50         '**'    => \&power,
  51         '<=>'   => \&spaceship,
  52         'neg'   => \&negate,
  53         '~'     => \&conjugate,
  54         'abs'   => \&abs,
  55         'sqrt'  => \&sqrt,
  56         'exp'   => \&exp,
  57         'log'   => \&log,
  58         'sin'   => \&sin,
  59         'cos'   => \&cos,
  60         'tan'   => \&tan,
  61         'atan2' => \&atan2,
  62         qw("" stringify);
  63
  64 #
  65 # Package globals
  66 #
  67
  68 $package = 'Math::Complex';             # Package name
  69 $display = 'cartesian';                 # Default display format
  70
  71 #
  72 # Object attributes (internal):
  73 #       cartesian       [real, imaginary] -- cartesian form
  74 #       polar           [rho, theta] -- polar form
  75 #       c_dirty         cartesian form not up-to-date
  76 #       p_dirty         polar form not up-to-date
  77 #       display         display format (package's global when not set)
  78 #
  79
  80 #
  81 # ->make
  82 #
  83 # Create a new complex number (cartesian form)
  84 #
  85 sub make {
  86         my $self = bless {}, shift;
  87         my ($re, $im) = @_;
  88         $self->{'cartesian'} = [$re, $im];
  89         $self->{c_dirty} = 0;
  90         $self->{p_dirty} = 1;
  91         return $self;
  92 }
  93
  94 #
  95 # ->emake
  96 #
  97 # Create a new complex number (exponential form)
  98 #
  99 sub emake {
 100         my $self = bless {}, shift;
 101         my ($rho, $theta) = @_;
 102         $theta += pi() if $rho < 0;
 103         $self->{'polar'} = [abs($rho), $theta];
 104         $self->{p_dirty} = 0;
 105         $self->{c_dirty} = 1;
 106         return $self;
 107 }
 108
 109 sub new { &make }               # For backward compatibility only.
 110
 111 #
 112 # cplx
 113 #
 114 # Creates a complex number from a (re, im) tuple.
 115 # This avoids the burden of writing Math::Complex->make(re, im).
 116 #
 117 sub cplx {
 118         my ($re, $im) = @_;
 119         return $package->make($re, defined $im ? $im : 0);
 120 }
 121
 122 #
 123 # cplxe
 124 #
 125 # Creates a complex number from a (rho, theta) tuple.
 126 # This avoids the burden of writing Math::Complex->emake(rho, theta).
 127 #
 128 sub cplxe {
 129         my ($rho, $theta) = @_;
 130         return $package->emake($rho, defined $theta ? $theta : 0);
 131 }
 132
 133 #
 134 # pi
 135 #
 136 # The number defined as 2 * pi = 360 degrees
 137 #
 138 sub pi () {
 139         $pi = 4 * atan2(1, 1) unless $pi;
 140         return $pi;
 141 }
 142
 143 #
 144 # i
 145 #
 146 # The number defined as i*i = -1;
 147 #
 148 sub i () {
 149         $i = bless {} unless $i;                # There can be only one i
 150         $i->{'cartesian'} = [0, 1];
 151         $i->{'polar'} = [1, pi/2];
 152         $i->{c_dirty} = 0;
 153         $i->{p_dirty} = 0;
 154         return $i;
 155 }
 156
 157 #
 158 # Attribute access/set routines
 159 #
 160
 161 sub cartesian {$_[0]->{c_dirty} ?
 162                    $_[0]->update_cartesian : $_[0]->{'cartesian'}}
 163 sub polar     {$_[0]->{p_dirty} ?
 164                    $_[0]->update_polar : $_[0]->{'polar'}}
 165
 166 sub set_cartesian { $_[0]->{p_dirty}++; $_[0]->{'cartesian'} = $_[1] }
 167 sub set_polar     { $_[0]->{c_dirty}++; $_[0]->{'polar'} = $_[1] }
 168
 169 #
 170 # ->update_cartesian
 171 #
 172 # Recompute and return the cartesian form, given accurate polar form.
 173 #
 174 sub update_cartesian {
 175         my $self = shift;
 176         my ($r, $t) = @{$self->{'polar'}};
 177         $self->{c_dirty} = 0;
 178         return $self->{'cartesian'} = [$r * cos $t, $r * sin $t];
 179 }
 180
 181 #
 182 #
 183 # ->update_polar
 184 #
 185 # Recompute and return the polar form, given accurate cartesian form.
 186 #
 187 sub update_polar {
 188         my $self = shift;
 189         my ($x, $y) = @{$self->{'cartesian'}};
 190         $self->{p_dirty} = 0;
 191         return $self->{'polar'} = [0, 0] if $x == 0 && $y == 0;
 192         return $self->{'polar'} = [sqrt($x*$x + $y*$y), atan2($y, $x)];
 193 }
 194
 195 #
 196 # (plus)
 197 #
 198 # Computes z1+z2.
 199 #
 200 sub plus {
 201         my ($z1, $z2, $regular) = @_;
 202         $z2 = cplx($z2, 0) unless ref $z2;
 203         my ($re1, $im1) = @{$z1->cartesian};
 204         my ($re2, $im2) = @{$z2->cartesian};
 205         unless (defined $regular) {
 206                 $z1->set_cartesian([$re1 + $re2, $im1 + $im2]);
 207                 return $z1;
 208         }
 209         return (ref $z1)->make($re1 + $re2, $im1 + $im2);
 210 }
 211
 212 #
 213 # (minus)
 214 #
 215 # Computes z1-z2.
 216 #
 217 sub minus {
 218         my ($z1, $z2, $inverted) = @_;
 219         $z2 = cplx($z2, 0) unless ref $z2;
 220         my ($re1, $im1) = @{$z1->cartesian};
 221         my ($re2, $im2) = @{$z2->cartesian};
 222         unless (defined $inverted) {
 223                 $z1->set_cartesian([$re1 - $re2, $im1 - $im2]);
 224                 return $z1;
 225         }
 226         return $inverted ?
 227                 (ref $z1)->make($re2 - $re1, $im2 - $im1) :
 228                 (ref $z1)->make($re1 - $re2, $im1 - $im2);
 229 }
 230
 231 #
 232 # (multiply)
 233 #
 234 # Computes z1*z2.
 235 #
 236 sub multiply {
 237         my ($z1, $z2, $regular) = @_;
 238         my ($r1, $t1) = @{$z1->polar};
 239         my ($r2, $t2) = ref $z2 ?
 240             @{$z2->polar} : (abs($z2), $z2 >= 0 ? 0 : pi);
 241         unless (defined $regular) {
 242                 $z1->set_polar([$r1 * $r2, $t1 + $t2]);
 243                 return $z1;
 244         }
 245         return (ref $z1)->emake($r1 * $r2, $t1 + $t2);
 246 }
 247
 248 #
 249 # divbyzero
 250 #
 251 # Die on division by zero.
 252 #
 253 sub divbyzero {
 254     warn "$_[0]: Division by zero.\n";
 255     warn "(Because in the definition of $_[0], $_[1] is 0)\n"
 256         if (defined $_[1]);
 257     my @up = caller(1);
 258     my $dmess = "Died at $up[1] line $up[2].\n";
 259     die $dmess;
 260 }
 261
 262 #
 263 # (divide)
 264 #
 265 # Computes z1/z2.
 266 #
 267 sub divide {
 268         my ($z1, $z2, $inverted) = @_;
 269         my ($r1, $t1) = @{$z1->polar};
 270         my ($r2, $t2) = ref $z2 ?
 271             @{$z2->polar} : (abs($z2), $z2 >= 0 ? 0 : pi);
 272         unless (defined $inverted) {
 273                 divbyzero "$z1/0" if ($r2 == 0);
 274                 $z1->set_polar([$r1 / $r2, $t1 - $t2]);
 275                 return $z1;
 276         }
 277         if ($inverted) {
 278                 divbyzero "$z2/0" if ($r1 == 0);
 279                 return (ref $z1)->emake($r2 / $r1, $t2 - $t1);
 280         } else {
 281                 divbyzero "$z1/0" if ($r2 == 0);
 282                 return (ref $z1)->emake($r1 / $r2, $t1 - $t2);
 283         }
 284 }
 285
 286 #
 287 # (power)
 288 #
 289 # Computes z1**z2 = exp(z2 * log z1)).
 290 #
 291 sub power {
 292         my ($z1, $z2, $inverted) = @_;
 293         return exp($z1 * log $z2) if defined $inverted && $inverted;
 294         return exp($z2 * log $z1);
 295 }
 296
 297 #
 298 # (spaceship)
 299 #
 300 # Computes z1 <=> z2.
 301 # Sorts on the real part first, then on the imaginary part. Thus 2-4i > 3+8i.
 302 #
 303 sub spaceship {
 304         my ($z1, $z2, $inverted) = @_;
 305         $z2 = cplx($z2, 0) unless ref $z2;
 306         my ($re1, $im1) = @{$z1->cartesian};
 307         my ($re2, $im2) = @{$z2->cartesian};
 308         my $sgn = $inverted ? -1 : 1;
 309         return $sgn * ($re1 <=> $re2) if $re1 != $re2;
 310         return $sgn * ($im1 <=> $im2);
 311 }
 312
 313 #
 314 # (negate)
 315 #
 316 # Computes -z.
 317 #
 318 sub negate {
 319         my ($z) = @_;
 320         if ($z->{c_dirty}) {
 321                 my ($r, $t) = @{$z->polar};
 322                 return (ref $z)->emake($r, pi + $t);
 323         }
 324         my ($re, $im) = @{$z->cartesian};
 325         return (ref $z)->make(-$re, -$im);
 326 }
 327
 328 #
 329 # (conjugate)
 330 #
 331 # Compute complex's conjugate.
 332 #
 333 sub conjugate {
 334         my ($z) = @_;
 335         if ($z->{c_dirty}) {
 336                 my ($r, $t) = @{$z->polar};
 337                 return (ref $z)->emake($r, -$t);
 338         }
 339         my ($re, $im) = @{$z->cartesian};
 340         return (ref $z)->make($re, -$im);
 341 }
 342
 343 #
 344 # (abs)
 345 #
 346 # Compute complex's norm (rho).
 347 #
 348 sub abs {
 349         my ($z) = @_;
 350         return abs($z) unless ref $z;
 351         my ($r, $t) = @{$z->polar};
 352         return abs($r);
 353 }
 354
 355 #
 356 # arg
 357 #
 358 # Compute complex's argument (theta).
 359 #
 360 sub arg {
 361         my ($z) = @_;
 362         return ($z < 0 ? pi : 0) unless ref $z;
 363         my ($r, $t) = @{$z->polar};
 364         return $t;
 365 }
 366
 367 #
 368 # (sqrt)
 369 #
 370 # Compute sqrt(z).
 371 #
 372 sub sqrt {
 373         my ($z) = @_;
 374         $z = cplx($z, 0) unless ref $z;
 375         my ($r, $t) = @{$z->polar};
 376         return (ref $z)->emake(sqrt($r), $t/2);
 377 }
 378
 379 #
 380 # cbrt
 381 #
 382 # Compute cbrt(z) (cubic root).
 383 #
 384 sub cbrt {
 385         my ($z) = @_;
 386         return cplx($z, 0) ** (1/3) unless ref $z;
 387         my ($r, $t) = @{$z->polar};
 388         return (ref $z)->emake($r**(1/3), $t/3);
 389 }
 390
 391 #
 392 # root
 393 #
 394 # Computes all nth root for z, returning an array whose size is n.
 395 # `n' must be a positive integer.
 396 #
 397 # The roots are given by (for k = 0..n-1):
 398 #
 399 # z^(1/n) = r^(1/n) (cos ((t+2 k pi)/n) + i sin ((t+2 k pi)/n))
 400 #
 401 sub root {
 402         my ($z, $n) = @_;
 403         $n = int($n + 0.5);
 404         return undef unless $n > 0;
 405         my ($r, $t) = ref $z ? @{$z->polar} : (abs($z), $z >= 0 ? 0 : pi);
 406         my @root;
 407         my $k;
 408         my $theta_inc = 2 * pi / $n;
 409         my $rho = $r ** (1/$n);
 410         my $theta;
 411         my $complex = ref($z) || $package;
 412         for ($k = 0, $theta = $t / $n; $k < $n; $k++, $theta += $theta_inc) {
 413                 push(@root, $complex->emake($rho, $theta));
 414         }
 415         return @root;
 416 }
 417
 418 #
 419 # Re
 420 #
 421 # Return Re(z).
 422 #
 423 sub Re {
 424         my ($z) = @_;
 425         return $z unless ref $z;
 426         my ($re, $im) = @{$z->cartesian};
 427         return $re;
 428 }
 429
 430 #
 431 # Im
 432 #
 433 # Return Im(z).
 434 #
 435 sub Im {
 436         my ($z) = @_;
 437         return 0 unless ref $z;
 438         my ($re, $im) = @{$z->cartesian};
 439         return $im;
 440 }
 441
 442 #
 443 # (exp)
 444 #
 445 # Computes exp(z).
 446 #
 447 sub exp {
 448         my ($z) = @_;
 449         $z = cplx($z, 0) unless ref $z;
 450         my ($x, $y) = @{$z->cartesian};
 451         return (ref $z)->emake(exp($x), $y);
 452 }
 453
 454 #
 455 # (log)
 456 #
 457 # Compute log(z).
 458 #
 459 sub log {
 460         my ($z) = @_;
 461         $z = cplx($z, 0) unless ref $z;
 462         my ($r, $t) = @{$z->polar};
 463         my ($x, $y) = @{$z->cartesian};
 464         $t -= 2 * pi if ($t >  pi() and $x < 0);
 465         $t += 2 * pi if ($t < -pi() and $x < 0);
 466         return (ref $z)->make(log($r), $t);
 467 }
 468
 469 #
 470 # ln
 471 #
 472 # Alias for log().
 473 #
 474 sub ln { Math::Complex::log(@_) }
 475
 476 #
 477 # log10
 478 #
 479 # Compute log10(z).
 480 #
 481 sub log10 {
 482         my ($z) = @_;
 483         my $ilog10 = 1 / log(10) unless defined $ilog10;
 484         return log(cplx($z, 0)) * $ilog10 unless ref $z;
 485         my ($r, $t) = @{$z->polar};
 486         return (ref $z)->make(log($r) * $ilog10, $t * $ilog10);
 487 }
 488
 489 #
 490 # logn
 491 #
 492 # Compute logn(z,n) = log(z) / log(n)
 493 #
 494 sub logn {
 495         my ($z, $n) = @_;
 496         $z = cplx($z, 0) unless ref $z;
 497         my $logn = $logn{$n};
 498         $logn = $logn{$n} = log($n) unless defined $logn;       # Cache log(n)
 499         return log($z) / $logn;
 500 }
 501
 502 #
 503 # (cos)
 504 #
 505 # Compute cos(z) = (exp(iz) + exp(-iz))/2.
 506 #
 507 sub cos {
 508         my ($z) = @_;
 509         my ($x, $y) = @{$z->cartesian};
 510         my $ey = exp($y);
 511         my $ey_1 = 1 / $ey;
 512         return (ref $z)->make(cos($x) * ($ey + $ey_1)/2,
 513                               sin($x) * ($ey_1 - $ey)/2);
 514 }
 515
 516 #
 517 # (sin)
 518 #
 519 # Compute sin(z) = (exp(iz) - exp(-iz))/2.
 520 #
 521 sub sin {
 522         my ($z) = @_;
 523         my ($x, $y) = @{$z->cartesian};
 524         my $ey = exp($y);
 525         my $ey_1 = 1 / $ey;
 526         return (ref $z)->make(sin($x) * ($ey + $ey_1)/2,
 527                               cos($x) * ($ey - $ey_1)/2);
 528 }
 529
 530 #
 531 # tan
 532 #
 533 # Compute tan(z) = sin(z) / cos(z).
 534 #
 535 sub tan {
 536         my ($z) = @_;
 537         my $cz = cos($z);
 538         divbyzero "tan($z)", "cos($z)" if ($cz == 0);
 539         return sin($z) / $cz;
 540 }
 541
 542 #
 543 # sec
 544 #
 545 # Computes the secant sec(z) = 1 / cos(z).
 546 #
 547 sub sec {
 548         my ($z) = @_;
 549         my $cz = cos($z);
 550         divbyzero "sec($z)", "cos($z)" if ($cz == 0);
 551         return 1 / $cz;
 552 }
 553
 554 #
 555 # csc
 556 #
 557 # Computes the cosecant csc(z) = 1 / sin(z).
 558 #
 559 sub csc {
 560         my ($z) = @_;
 561         my $sz = sin($z);
 562         divbyzero "csc($z)", "sin($z)" if ($sz == 0);
 563         return 1 / $sz;
 564 }
 565
 566 #
 567 # cosec
 568 #
 569 # Alias for csc().
 570 #
 571 sub cosec { Math::Complex::csc(@_) }
 572
 573 #
 574 # cot
 575 #
 576 # Computes cot(z) = 1 / tan(z).
 577 #
 578 sub cot {
 579         my ($z) = @_;
 580         my $sz = sin($z);
 581         divbyzero "cot($z)", "sin($z)" if ($sz == 0);
 582         return cos($z) / $sz;
 583 }
 584
 585 #
 586 # cotan
 587 #
 588 # Alias for cot().
 589 #
 590 sub cotan { Math::Complex::cot(@_) }
 591
 592 #
 593 # acos
 594 #
 595 # Computes the arc cosine acos(z) = -i log(z + sqrt(z*z-1)).
 596 #
 597 sub acos {
 598         my ($z) = @_;
 599         $z = cplx($z, 0) unless ref $z;
 600         return ~i * log($z + (Re($z) * Im($z) > 0 ? 1 : -1) * sqrt($z*$z - 1));
 601 }
 602
 603 #
 604 # asin
 605 #
 606 # Computes the arc sine asin(z) = -i log(iz + sqrt(1-z*z)).
 607 #
 608 sub asin {
 609         my ($z) = @_;
 610         $z = cplx($z, 0) unless ref $z;
 611         return ~i * log(i * $z + sqrt(1 - $z*$z));
 612 }
 613
 614 #
 615 # atan
 616 #
 617 # Computes the arc tangent atan(z) = i/2 log((i+z) / (i-z)).
 618 #
 619 sub atan {
 620         my ($z) = @_;
 621         divbyzero "atan($z)", "i - $z" if ($z == i);
 622         return i/2*log((i + $z) / (i - $z));
 623 }
 624
 625 #
 626 # asec
 627 #
 628 # Computes the arc secant asec(z) = acos(1 / z).
 629 #
 630 sub asec {
 631         my ($z) = @_;
 632         return acos(1 / $z);
 633 }
 634
 635 #
 636 # acosec
 637 #
 638 # Computes the arc cosecant sec(z) = asin(1 / z).
 639 #
 640 sub acosec {
 641         my ($z) = @_;
 642         return asin(1 / $z);
 643 }
 644
 645 #
 646 # acsc
 647 #
 648 # Alias for acosec().
 649 #
 650 sub acsc { Math::Complex::acosec(@_) }
 651
 652 #
 653 # acot
 654 #
 655 # Computes the arc cotangent acot(z) = -i/2 log((i+z) / (z-i))
 656 #
 657 sub acot {
 658         my ($z) = @_;
 659         divbyzero "acot($z)", "$z - i" if ($z == i);
 660         return i/-2 * log((i + $z) / ($z - i));
 661 }
 662
 663 #
 664 # acotan
 665 #
 666 # Alias for acot().
 667 #
 668 sub acotan { Math::Complex::acot(@_) }
 669
 670 #
 671 # cosh
 672 #
 673 # Computes the hyperbolic cosine cosh(z) = (exp(z) + exp(-z))/2.
 674 #
 675 sub cosh {
 676         my ($z) = @_;
 677         $z = cplx($z, 0) unless ref $z;
 678         my ($x, $y) = @{$z->cartesian};
 679         my $ex = exp($x);
 680         my $ex_1 = 1 / $ex;
 681         return ($ex + $ex_1)/2 unless ref $z;
 682         return (ref $z)->make(cos($y) * ($ex + $ex_1)/2,
 683                               sin($y) * ($ex - $ex_1)/2);
 684 }
 685
 686 #
 687 # sinh
 688 #
 689 # Computes the hyperbolic sine sinh(z) = (exp(z) - exp(-z))/2.
 690 #
 691 sub sinh {
 692         my ($z) = @_;
 693         $z = cplx($z, 0) unless ref $z;
 694         my ($x, $y) = @{$z->cartesian};
 695         my $ex = exp($x);
 696         my $ex_1 = 1 / $ex;
 697         return ($ex - $ex_1)/2 unless ref $z;
 698         return (ref $z)->make(cos($y) * ($ex - $ex_1)/2,
 699                               sin($y) * ($ex + $ex_1)/2);
 700 }
 701
 702 #
 703 # tanh
 704 #
 705 # Computes the hyperbolic tangent tanh(z) = sinh(z) / cosh(z).
 706 #
 707 sub tanh {
 708         my ($z) = @_;
 709         my $cz = cosh($z);
 710         divbyzero "tanh($z)", "cosh($z)" if ($cz == 0);
 711         return sinh($z) / $cz;
 712 }
 713
 714 #
 715 # sech
 716 #
 717 # Computes the hyperbolic secant sech(z) = 1 / cosh(z).
 718 #
 719 sub sech {
 720         my ($z) = @_;
 721         my $cz = cosh($z);
 722         divbyzero "sech($z)", "cosh($z)" if ($cz == 0);
 723         return 1 / $cz;
 724 }
 725
 726 #
 727 # csch
 728 #
 729 # Computes the hyperbolic cosecant csch(z) = 1 / sinh(z).
 730 #
 731 sub csch {
 732         my ($z) = @_;
 733         my $sz = sinh($z);
 734         divbyzero "csch($z)", "sinh($z)" if ($sz == 0);
 735         return 1 / $sz;
 736 }
 737
 738 #
 739 # cosech
 740 #
 741 # Alias for csch().
 742 #
 743 sub cosech { Math::Complex::csch(@_) }
 744
 745 #
 746 # coth
 747 #
 748 # Computes the hyperbolic cotangent coth(z) = cosh(z) / sinh(z).
 749 #
 750 sub coth {
 751         my ($z) = @_;
 752         my $sz = sinh($z);
 753         divbyzero "coth($z)", "sinh($z)" if ($sz == 0);
 754         return cosh($z) / $sz;
 755 }
 756
 757 #
 758 # cotanh
 759 #
 760 # Alias for coth().
 761 #
 762 sub cotanh { Math::Complex::coth(@_) }
 763
 764 #
 765 # acosh
 766 #
 767 # Computes the arc hyperbolic cosine acosh(z) = log(z + sqrt(z*z-1)).
 768 #
 769 sub acosh {
 770         my ($z) = @_;
 771         $z = cplx($z, 0) unless ref $z; # asinh(-2)
 772         return log($z + sqrt($z*$z - 1));
 773 }
 774
 775 #
 776 # asinh
 777 #
 778 # Computes the arc hyperbolic sine asinh(z) = log(z + sqrt(z*z-1))
 779 #
 780 sub asinh {
 781         my ($z) = @_;
 782         $z = cplx($z, 0) unless ref $z; # asinh(-2)
 783         return log($z + sqrt($z*$z + 1));
 784 }
 785
 786 #
 787 # atanh
 788 #
 789 # Computes the arc hyperbolic tangent atanh(z) = 1/2 log((1+z) / (1-z)).
 790 #
 791 sub atanh {
 792         my ($z) = @_;
 793         $z = cplx($z, 0) unless ref $z; # atanh(-2)
 794         divbyzero 'atanh(1)', "1 - $z" if ($z == 1);
 795         my $cz = (1 + $z) / (1 - $z);
 796         return log($cz) / 2;
 797 }
 798
 799 #
 800 # asech
 801 #
 802 # Computes the hyperbolic arc secant asech(z) = acosh(1 / z).
 803 #
 804 sub asech {
 805         my ($z) = @_;
 806         divbyzero 'asech(0)', $z if ($z == 0);
 807         return acosh(1 / $z);
 808 }
 809
 810 #
 811 # acsch
 812 #
 813 # Computes the hyperbolic arc cosecant acsch(z) = asinh(1 / z).
 814 #
 815 sub acsch {
 816         my ($z) = @_;
 817         divbyzero 'acsch(0)', $z if ($z == 0);
 818         return asinh(1 / $z);
 819 }
 820
 821 #
 822 # acosech
 823 #
 824 # Alias for acosh().
 825 #
 826 sub acosech { Math::Complex::acsch(@_) }
 827
 828 #
 829 # acoth
 830 #
 831 # Computes the arc hyperbolic cotangent acoth(z) = 1/2 log((1+z) / (z-1)).
 832 #
 833 sub acoth {
 834         my ($z) = @_;
 835         $z = cplx($z, 0) unless ref $z; # acoth(-2)
 836         divbyzero 'acoth(1)', "$z - 1" if ($z == 1);
 837         my $cz = (1 + $z) / ($z - 1);
 838         return log($cz) / 2;
 839 }
 840
 841 #
 842 # acotanh
 843 #
 844 # Alias for acot().
 845 #
 846 sub acotanh { Math::Complex::acoth(@_) }
 847
 848 #
 849 # (atan2)
 850 #
 851 # Compute atan(z1/z2).
 852 #
 853 sub atan2 {
 854         my ($z1, $z2, $inverted) = @_;
 855         my ($re1, $im1) = @{$z1->cartesian};
 856         my ($re2, $im2) = @{$z2->cartesian};
 857         my $tan;
 858         if (defined $inverted && $inverted) {   # atan(z2/z1)
 859                 return pi * ($re2 > 0 ? 1 : -1) if $re1 == 0 && $im1 == 0;
 860                 $tan = $z2 / $z1;
 861         } else {
 862                 return pi * ($re1 > 0 ? 1 : -1) if $re2 == 0 && $im2 == 0;
 863                 $tan = $z1 / $z2;
 864         }
 865         return atan($tan);
 866 }
 867
 868 #
 869 # display_format
 870 # ->display_format
 871 #
 872 # Set (fetch if no argument) display format for all complex numbers that
 873 # don't happen to have overrriden it via ->display_format
 874 #
 875 # When called as a method, this actually sets the display format for
 876 # the current object.
 877 #
 878 # Valid object formats are 'c' and 'p' for cartesian and polar. The first
 879 # letter is used actually, so the type can be fully spelled out for clarity.
 880 #
 881 sub display_format {
 882         my $self = shift;
 883         my $format = undef;
 884
 885         if (ref $self) {                        # Called as a method
 886                 $format = shift;
 887         } else {                                # Regular procedure call
 888                 $format = $self;
 889                 undef $self;
 890         }
 891
 892         if (defined $self) {
 893                 return defined $self->{display} ? $self->{display} : $display
 894                         unless defined $format;
 895                 return $self->{display} = $format;
 896         }
 897
 898         return $display unless defined $format;
 899         return $display = $format;
 900 }
 901
 902 #
 903 # (stringify)
 904 #
 905 # Show nicely formatted complex number under its cartesian or polar form,
 906 # depending on the current display format:
 907 #
 908 # . If a specific display format has been recorded for this object, use it.
 909 # . Otherwise, use the generic current default for all complex numbers,
 910 #   which is a package global variable.
 911 #
 912 sub stringify {
 913         my ($z) = shift;
 914         my $format;
 915
 916         $format = $display;
 917         $format = $z->{display} if defined $z->{display};
 918
 919         return $z->stringify_polar if $format =~ /^p/i;
 920         return $z->stringify_cartesian;
 921 }
 922
 923 #
 924 # ->stringify_cartesian
 925 #
 926 # Stringify as a cartesian representation 'a+bi'.
 927 #
 928 sub stringify_cartesian {
 929         my $z  = shift;
 930         my ($x, $y) = @{$z->cartesian};
 931         my ($re, $im);
 932
 933         $x = int($x + ($x < 0 ? -1 : 1) * 1e-14)
 934                 if int(abs($x)) != int(abs($x) + 1e-14);
 935         $y = int($y + ($y < 0 ? -1 : 1) * 1e-14)
 936                 if int(abs($y)) != int(abs($y) + 1e-14);
 937
 938         $re = "$x" if abs($x) >= 1e-14;
 939         if ($y == 1)                            { $im = 'i' }
 940         elsif ($y == -1)                        { $im = '-i' }
 941         elsif (abs($y) >= 1e-14)        { $im = $y . "i" }
 942
 943         my $str = '';
 944         $str = $re if defined $re;
 945         $str .= "+$im" if defined $im;
 946         $str =~ s/\+-/-/;
 947         $str =~ s/^\+//;
 948         $str = '0' unless $str;
 949
 950         return $str;
 951 }
 952
 953 #
 954 # ->stringify_polar
 955 #
 956 # Stringify as a polar representation '[r,t]'.
 957 #
 958 sub stringify_polar {
 959         my $z  = shift;
 960         my ($r, $t) = @{$z->polar};
 961         my $theta;
 962         my $eps = 1e-14;
 963
 964         return '[0,0]' if $r <= $eps;
 965
 966         my $tpi = 2 * pi;
 967         my $nt = $t / $tpi;
 968         $nt = ($nt - int($nt)) * $tpi;
 969         $nt += $tpi if $nt < 0;                 # Range [0, 2pi]
 970
 971         if (abs($nt) <= $eps)           { $theta = 0 }
 972         elsif (abs(pi-$nt) <= $eps)     { $theta = 'pi' }
 973
 974         if (defined $theta) {
 975                 $r = int($r + ($r < 0 ? -1 : 1) * $eps)
 976                         if int(abs($r)) != int(abs($r) + $eps);
 977                 $theta = int($theta + ($theta < 0 ? -1 : 1) * $eps)
 978                         if ($theta ne 'pi' and
 979                             int(abs($theta)) != int(abs($theta) + $eps));
 980                 return "\[$r,$theta\]";
 981         }
 982
 983         #
 984         # Okay, number is not a real. Try to identify pi/n and friends...
 985         #
 986
 987         $nt -= $tpi if $nt > pi;
 988         my ($n, $k, $kpi);
 989
 990         for ($k = 1, $kpi = pi; $k < 10; $k++, $kpi += pi) {
 991                 $n = int($kpi / $nt + ($nt > 0 ? 1 : -1) * 0.5);
 992                 if (abs($kpi/$n - $nt) <= $eps) {
 993                         $theta = ($nt < 0 ? '-':'').
 994                                  ($k == 1 ? 'pi':"${k}pi").'/'.abs($n);
 995                         last;
 996                 }
 997         }
 998
 999         $theta = $nt unless defined $theta;
1000
1001         $r = int($r + ($r < 0 ? -1 : 1) * $eps)
1002                 if int(abs($r)) != int(abs($r) + $eps);
1003         $theta = int($theta + ($theta < 0 ? -1 : 1) * $eps)
1004                 if ($theta !~ m(^-?\d*pi/\d+$) and
1005                     int(abs($theta)) != int(abs($theta) + $eps));
1006
1007         return "\[$r,$theta\]";
1008 }
1009
1010 1;
1011 __END__
1012
1013 =head1 NAME
1014
1015 Math::Complex - complex numbers and associated mathematical functions
1016
1017 =head1 SYNOPSIS
1018
1019         use Math::Complex;
1020
1021         $z = Math::Complex->make(5, 6);
1022         $t = 4 - 3*i + $z;
1023         $j = cplxe(1, 2*pi/3);
1024
1025 =head1 DESCRIPTION
1026
1027 This package lets you create and manipulate complex numbers. By default,
1028 I<Perl> limits itself to real numbers, but an extra C<use> statement brings
1029 full complex support, along with a full set of mathematical functions
1030 typically associated with and/or extended to complex numbers.
1031
1032 If you wonder what complex numbers are, they were invented to be able to solve
1033 the following equation:
1034
1035         x*x = -1
1036
1037 and by definition, the solution is noted I<i> (engineers use I<j> instead since
1038 I<i> usually denotes an intensity, but the name does not matter). The number
1039 I<i> is a pure I<imaginary> number.
1040
1041 The arithmetics with pure imaginary numbers works just like you would expect
1042 it with real numbers... you just have to remember that
1043
1044         i*i = -1
1045
1046 so you have:
1047
1048         5i + 7i = i * (5 + 7) = 12i
1049         4i - 3i = i * (4 - 3) = i
1050         4i * 2i = -8
1051         6i / 2i = 3
1052         1 / i = -i
1053
1054 Complex numbers are numbers that have both a real part and an imaginary
1055 part, and are usually noted:
1056
1057         a + bi
1058
1059 where C<a> is the I<real> part and C<b> is the I<imaginary> part. The
1060 arithmetic with complex numbers is straightforward. You have to
1061 keep track of the real and the imaginary parts, but otherwise the
1062 rules used for real numbers just apply:
1063
1064         (4 + 3i) + (5 - 2i) = (4 + 5) + i(3 - 2) = 9 + i
1065         (2 + i) * (4 - i) = 2*4 + 4i -2i -i*i = 8 + 2i + 1 = 9 + 2i
1066
1067 A graphical representation of complex numbers is possible in a plane
1068 (also called the I<complex plane>, but it's really a 2D plane).
1069 The number
1070
1071         z = a + bi
1072
1073 is the point whose coordinates are (a, b). Actually, it would
1074 be the vector originating from (0, 0) to (a, b). It follows that the addition
1075 of two complex numbers is a vectorial addition.
1076
1077 Since there is a bijection between a point in the 2D plane and a complex
1078 number (i.e. the mapping is unique and reciprocal), a complex number
1079 can also be uniquely identified with polar coordinates:
1080
1081         [rho, theta]
1082
1083 where C<rho> is the distance to the origin, and C<theta> the angle between
1084 the vector and the I<x> axis. There is a notation for this using the
1085 exponential form, which is:
1086
1087         rho * exp(i * theta)
1088
1089 where I<i> is the famous imaginary number introduced above. Conversion
1090 between this form and the cartesian form C<a + bi> is immediate:
1091
1092         a = rho * cos(theta)
1093         b = rho * sin(theta)
1094
1095 which is also expressed by this formula:
1096
1097         z = rho * exp(i * theta) = rho * (cos theta + i * sin theta)
1098
1099 In other words, it's the projection of the vector onto the I<x> and I<y>
1100 axes. Mathematicians call I<rho> the I<norm> or I<modulus> and I<theta>
1101 the I<argument> of the complex number. The I<norm> of C<z> will be
1102 noted C<abs(z)>.
1103
1104 The polar notation (also known as the trigonometric
1105 representation) is much more handy for performing multiplications and
1106 divisions of complex numbers, whilst the cartesian notation is better
1107 suited for additions and substractions. Real numbers are on the I<x>
1108 axis, and therefore I<theta> is zero.
1109
1110 All the common operations that can be performed on a real number have
1111 been defined to work on complex numbers as well, and are merely
1112 I<extensions> of the operations defined on real numbers. This means
1113 they keep their natural meaning when there is no imaginary part, provided
1114 the number is within their definition set.
1115
1116 For instance, the C<sqrt> routine which computes the square root of
1117 its argument is only defined for positive real numbers and yields a
1118 positive real number (it is an application from B<R+> to B<R+>).
1119 If we allow it to return a complex number, then it can be extended to
1120 negative real numbers to become an application from B<R> to B<C> (the
1121 set of complex numbers):
1122
1123         sqrt(x) = x >= 0 ? sqrt(x) : sqrt(-x)*i
1124
1125 It can also be extended to be an application from B<C> to B<C>,
1126 whilst its restriction to B<R> behaves as defined above by using
1127 the following definition:
1128
1129         sqrt(z = [r,t]) = sqrt(r) * exp(i * t/2)
1130
1131 Indeed, a negative real number can be noted C<[x,pi]>
1132 (the modulus I<x> is always positive, so C<[x,pi]> is really C<-x>, a
1133 negative number)
1134 and the above definition states that
1135
1136         sqrt([x,pi]) = sqrt(x) * exp(i*pi/2) = [sqrt(x),pi/2] = sqrt(x)*i
1137
1138 which is exactly what we had defined for negative real numbers above.
1139
1140 All the common mathematical functions defined on real numbers that
1141 are extended to complex numbers share that same property of working
1142 I<as usual> when the imaginary part is zero (otherwise, it would not
1143 be called an extension, would it?).
1144
1145 A I<new> operation possible on a complex number that is
1146 the identity for real numbers is called the I<conjugate>, and is noted
1147 with an horizontal bar above the number, or C<~z> here.
1148
1149          z = a + bi
1150         ~z = a - bi
1151
1152 Simple... Now look:
1153
1154         z * ~z = (a + bi) * (a - bi) = a*a + b*b
1155
1156 We saw that the norm of C<z> was noted C<abs(z)> and was defined as the
1157 distance to the origin, also known as:
1158
1159         rho = abs(z) = sqrt(a*a + b*b)
1160
1161 so
1162
1163         z * ~z = abs(z) ** 2
1164
1165 If z is a pure real number (i.e. C<b == 0>), then the above yields:
1166
1167         a * a = abs(a) ** 2
1168
1169 which is true (C<abs> has the regular meaning for real number, i.e. stands
1170 for the absolute value). This example explains why the norm of C<z> is
1171 noted C<abs(z)>: it extends the C<abs> function to complex numbers, yet
1172 is the regular C<abs> we know when the complex number actually has no
1173 imaginary part... This justifies I<a posteriori> our use of the C<abs>
1174 notation for the norm.
1175
1176 =head1 OPERATIONS
1177
1178 Given the following notations:
1179
1180         z1 = a + bi = r1 * exp(i * t1)
1181         z2 = c + di = r2 * exp(i * t2)
1182         z = <any complex or real number>
1183
1184 the following (overloaded) operations are supported on complex numbers:
1185
1186         z1 + z2 = (a + c) + i(b + d)
1187         z1 - z2 = (a - c) + i(b - d)
1188         z1 * z2 = (r1 * r2) * exp(i * (t1 + t2))
1189         z1 / z2 = (r1 / r2) * exp(i * (t1 - t2))
1190         z1 ** z2 = exp(z2 * log z1)
1191         ~z1 = a - bi
1192         abs(z1) = r1 = sqrt(a*a + b*b)
1193         sqrt(z1) = sqrt(r1) * exp(i * t1/2)
1194         exp(z1) = exp(a) * exp(i * b)
1195         log(z1) = log(r1) + i*t1
1196         sin(z1) = 1/2i (exp(i * z1) - exp(-i * z1))
1197         cos(z1) = 1/2 (exp(i * z1) + exp(-i * z1))
1198         abs(z1) = r1
1199         atan2(z1, z2) = atan(z1/z2)
1200
1201 The following extra operations are supported on both real and complex
1202 numbers:
1203
1204         Re(z) = a
1205         Im(z) = b
1206         arg(z) = t
1207
1208         cbrt(z) = z ** (1/3)
1209         log10(z) = log(z) / log(10)
1210         logn(z, n) = log(z) / log(n)
1211
1212         tan(z) = sin(z) / cos(z)
1213
1214         csc(z) = 1 / sin(z)
1215         sec(z) = 1 / cos(z)
1216         cot(z) = 1 / tan(z)
1217
1218         asin(z) = -i * log(i*z + sqrt(1-z*z))
1219         acos(z) = -i * log(z + sqrt(z*z-1))
1220         atan(z) = i/2 * log((i+z) / (i-z))
1221
1222         acsc(z) = asin(1 / z)
1223         asec(z) = acos(1 / z)
1224         acot(z) = -i/2 * log((i+z) / (z-i))
1225
1226         sinh(z) = 1/2 (exp(z) - exp(-z))
1227         cosh(z) = 1/2 (exp(z) + exp(-z))
1228         tanh(z) = sinh(z) / cosh(z) = (exp(z) - exp(-z)) / (exp(z) + exp(-z))
1229
1230         csch(z) = 1 / sinh(z)
1231         sech(z) = 1 / cosh(z)
1232         coth(z) = 1 / tanh(z)
1233
1234         asinh(z) = log(z + sqrt(z*z+1))
1235         acosh(z) = log(z + sqrt(z*z-1))
1236         atanh(z) = 1/2 * log((1+z) / (1-z))
1237
1238         acsch(z) = asinh(1 / z)
1239         asech(z) = acosh(1 / z)
1240         acoth(z) = atanh(1 / z) = 1/2 * log((1+z) / (z-1))
1241
1242 I<log>, I<csc>, I<cot>, I<acsc>, I<acot>, I<csch>, I<coth>,
1243 I<acosech>, I<acotanh>, have aliases I<ln>, I<cosec>, I<cotan>,
1244 I<acosec>, I<acotan>, I<cosech>, I<cotanh>, I<acosech>, I<acotanh>,
1245 respectively.
1246
1247 The I<root> function is available to compute all the I<n>
1248 roots of some complex, where I<n> is a strictly positive integer.
1249 There are exactly I<n> such roots, returned as a list. Getting the
1250 number mathematicians call C<j> such that:
1251
1252         1 + j + j*j = 0;
1253
1254 is a simple matter of writing:
1255
1256         $j = ((root(1, 3))[1];
1257
1258 The I<k>th root for C<z = [r,t]> is given by:
1259
1260         (root(z, n))[k] = r**(1/n) * exp(i * (t + 2*k*pi)/n)
1261
1262 The I<spaceship> comparison operator is also defined. In order to
1263 ensure its restriction to real numbers is conform to what you would
1264 expect, the comparison is run on the real part of the complex number
1265 first, and imaginary parts are compared only when the real parts
1266 match.
1267
1268 =head1 CREATION
1269
1270 To create a complex number, use either:
1271
1272         $z = Math::Complex->make(3, 4);
1273         $z = cplx(3, 4);
1274
1275 if you know the cartesian form of the number, or
1276
1277         $z = 3 + 4*i;
1278
1279 if you like. To create a number using the trigonometric form, use either:
1280
1281         $z = Math::Complex->emake(5, pi/3);
1282         $x = cplxe(5, pi/3);
1283
1284 instead. The first argument is the modulus, the second is the angle
1285 (in radians, the full circle is 2*pi).  (Mnmemonic: C<e> is used as a
1286 notation for complex numbers in the trigonometric form).
1287
1288 It is possible to write:
1289
1290         $x = cplxe(-3, pi/4);
1291
1292 but that will be silently converted into C<[3,-3pi/4]>, since the modulus
1293 must be positive (it represents the distance to the origin in the complex
1294 plane).
1295
1296 =head1 STRINGIFICATION
1297
1298 When printed, a complex number is usually shown under its cartesian
1299 form I<a+bi>, but there are legitimate cases where the polar format
1300 I<[r,t]> is more appropriate.
1301
1302 By calling the routine C<Math::Complex::display_format> and supplying either
1303 C<"polar"> or C<"cartesian">, you override the default display format,
1304 which is C<"cartesian">. Not supplying any argument returns the current
1305 setting.
1306
1307 This default can be overridden on a per-number basis by calling the
1308 C<display_format> method instead. As before, not supplying any argument
1309 returns the current display format for this number. Otherwise whatever you
1310 specify will be the new display format for I<this> particular number.
1311
1312 For instance:
1313
1314         use Math::Complex;
1315
1316         Math::Complex::display_format('polar');
1317         $j = ((root(1, 3))[1];
1318         print "j = $j\n";               # Prints "j = [1,2pi/3]
1319         $j->display_format('cartesian');
1320         print "j = $j\n";               # Prints "j = -0.5+0.866025403784439i"
1321
1322 The polar format attempts to emphasize arguments like I<k*pi/n>
1323 (where I<n> is a positive integer and I<k> an integer within [-9,+9]).
1324
1325 =head1 USAGE
1326
1327 Thanks to overloading, the handling of arithmetics with complex numbers
1328 is simple and almost transparent.
1329
1330 Here are some examples:
1331
1332         use Math::Complex;
1333
1334         $j = cplxe(1, 2*pi/3);  # $j ** 3 == 1
1335         print "j = $j, j**3 = ", $j ** 3, "\n";
1336         print "1 + j + j**2 = ", 1 + $j + $j**2, "\n";
1337
1338         $z = -16 + 0*i;                 # Force it to be a complex
1339         print "sqrt($z) = ", sqrt($z), "\n";
1340
1341         $k = exp(i * 2*pi/3);
1342         print "$j - $k = ", $j - $k, "\n";
1343
1344 =head1 CAVEATS
1345
1346 The division (/) and the following functions
1347
1348         tan
1349         sec
1350         csc
1351         cot
1352         atan
1353         acot
1354         tanh
1355         sech
1356         csch
1357         coth
1358         atanh
1359         asech
1360         acsch
1361         acoth
1362
1363 cannot be computed for all arguments because that would mean dividing
1364 by zero. These situations cause fatal runtime errors looking like this
1365
1366         cot(0): Division by zero.
1367         (Because in the definition of cot(0), sin(0) is 0)
1368         Died at ...
1369
1370 =head1 BUGS
1371
1372 Saying C<use Math::Complex;> exports many mathematical routines in the caller
1373 environment.  This is construed as a feature by the Author, actually... ;-)
1374
1375 The code is not optimized for speed, although we try to use the cartesian
1376 form for addition-like operators and the trigonometric form for all
1377 multiplication-like operators.
1378
1379 The arg() routine does not ensure the angle is within the range [-pi,+pi]
1380 (a side effect caused by multiplication and division using the trigonometric
1381 representation).
1382
1383 All routines expect to be given real or complex numbers. Don't attempt to
1384 use BigFloat, since Perl has currently no rule to disambiguate a '+'
1385 operation (for instance) between two overloaded entities.
1386
1387 =head1 AUTHORS
1388
1389         Raphael Manfredi <F<Raphael_Manfredi@grenoble.hp.com>>
1390         Jarkko Hietaniemi <F<jhi@iki.fi>>