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
12 use vars qw($VERSION @ISA @EXPORT @EXPORT_OK %EXPORT_TAGS $Inf $ExpInf);
21 4 => '1.70141183460469229e+38',
22 8 => '1.7976931348623157e+308',
23 # AFAICT the 10, 12, and 16-byte long doubles
24 # all have the same maximum.
25 10 => '1.1897314953572317650857593266280070162E+4932',
26 12 => '1.1897314953572317650857593266280070162E+4932',
27 16 => '1.1897314953572317650857593266280070162E+4932',
29 my $nvsize = $Config{nvsize} ||
30 ($Config{uselongdouble} && $Config{longdblsize}) ||
32 die "Math::Complex: Could not figure out nvsize\n"
33 unless defined $nvsize;
34 die "Math::Complex: Cannot not figure out max nv (nvsize = $nvsize)\n"
35 unless defined $DBL_MAX{$nvsize};
36 my $DBL_MAX = eval $DBL_MAX{$nvsize};
37 die "Math::Complex: Could not figure out max nv (nvsize = $nvsize)\n"
38 unless defined $DBL_MAX;
39 my $BIGGER_THAN_THIS = 1e30; # Must find something bigger than this.
40 if ($^O eq 'unicosmk') {
43 local $SIG{FPE} = { };
45 # We do want an arithmetic overflow, Inf INF inf Infinity.
47 'exp(99999)', # Enough even with 128-bit long doubles.
57 my $i = eval "$t+1.0";
58 if (defined $i && $i > $BIGGER_THAN_THIS) {
63 $Inf = $DBL_MAX unless defined $Inf; # Oh well, close enough.
64 die "Math::Complex: Could not get Infinity"
65 unless $Inf > $BIGGER_THAN_THIS;
68 # print "# On this machine, Inf = '$Inf'\n";
71 use Scalar::Util qw(set_prototype);
74 no warnings 'syntax'; # To avoid the (_) warnings.
77 # For certain functions that we override, in 5.10 or better
78 # we can set a smarter prototype that will handle the lexical $_
79 # (also a 5.10+ feature).
81 set_prototype \&abs, '_';
82 set_prototype \&cos, '_';
83 set_prototype \&exp, '_';
84 set_prototype \&log, '_';
85 set_prototype \&sin, '_';
86 set_prototype \&sqrt, '_';
93 # Regular expression for floating point numbers.
94 # These days we could use Scalar::Util::lln(), I guess.
95 my $gre = qr'\s*([\+\-]?(?:(?:(?:\d+(?:_\d+)*(?:\.\d*(?:_\d+)*)?|\.\d+(?:_\d+)*)(?:[eE][\+\-]?\d+(?:_\d+)*)?))|inf)'i;
104 csc cosec sec cot cotan
106 acsc acosec asec acot acotan
108 csch cosech sech coth cotanh
110 acsch acosech asech acoth acotanh
114 i Re Im rho theta arg
122 my @pi = qw(pi pi2 pi4 pip2 pip4 Inf);
138 '<=>' => \&_spaceship,
149 '""' => \&_stringify;
155 my %DISPLAY_FORMAT = ('style' => 'cartesian',
156 'polar_pretty_print' => 1);
157 my $eps = 1e-14; # Epsilon
160 # Object attributes (internal):
161 # cartesian [real, imaginary] -- cartesian form
162 # polar [rho, theta] -- polar form
163 # c_dirty cartesian form not up-to-date
164 # p_dirty polar form not up-to-date
165 # display display format (package's global when not set)
168 # Die on bad *make() arguments.
171 die "@{[(caller(1))[3]]}: Cannot take $_[0] of '$_[1]'.\n";
178 if ($arg =~ /^$gre$/) {
180 } elsif ($arg =~ /^(?:$gre)?$gre\s*i\s*$/) {
181 ($p, $q) = ($1 || 0, $2);
182 } elsif ($arg =~ /^\s*\(\s*$gre\s*(?:,\s*$gre\s*)?\)\s*$/) {
183 ($p, $q) = ($1, $2 || 0);
188 $p =~ s/^(-?)inf$/"${1}9**9**9"/e;
190 $q =~ s/^(-?)inf$/"${1}9**9**9"/e;
200 if ($arg =~ /^\s*\[\s*$gre\s*(?:,\s*$gre\s*)?\]\s*$/) {
201 ($p, $q) = ($1, $2 || 0);
202 } elsif ($arg =~ m!^\s*\[\s*$gre\s*(?:,\s*([-+]?\d*\s*)?pi(?:/\s*(\d+))?\s*)?\]\s*$!) {
203 ($p, $q) = ($1, ($2 eq '-' ? -1 : ($2 || 1)) * pi() / ($3 || 1));
204 } elsif ($arg =~ /^\s*\[\s*$gre\s*\]\s*$/) {
206 } elsif ($arg =~ /^\s*$gre\s*$/) {
213 $p =~ s/^(-?)inf$/"${1}9**9**9"/e;
214 $q =~ s/^(-?)inf$/"${1}9**9**9"/e;
223 # Create a new complex number (cartesian form)
226 my $self = bless {}, shift;
231 return (ref $self)->emake($_[0])
232 if ($_[0] =~ /^\s*\[/);
233 ($re, $im) = _make($_[0]);
238 _cannot_make("real part", $re) unless $re =~ /^$gre$/;
241 _cannot_make("imaginary part", $im) unless $im =~ /^$gre$/;
242 $self->_set_cartesian([$re, $im ]);
243 $self->display_format('cartesian');
251 # Create a new complex number (exponential form)
254 my $self = bless {}, shift;
257 ($rho, $theta) = (0, 0);
259 return (ref $self)->make($_[0])
260 if ($_[0] =~ /^\s*\(/ || $_[0] =~ /i\s*$/);
261 ($rho, $theta) = _emake($_[0]);
265 if (defined $rho && defined $theta) {
268 $theta = ($theta <= 0) ? $theta + pi() : $theta - pi();
272 _cannot_make("rho", $rho) unless $rho =~ /^$gre$/;
275 _cannot_make("theta", $theta) unless $theta =~ /^$gre$/;
276 $self->_set_polar([$rho, $theta]);
277 $self->display_format('polar');
282 sub new { &make } # For backward compatibility only.
287 # Creates a complex number from a (re, im) tuple.
288 # This avoids the burden of writing Math::Complex->make(re, im).
291 return __PACKAGE__->make(@_);
297 # Creates a complex number from a (rho, theta) tuple.
298 # This avoids the burden of writing Math::Complex->emake(rho, theta).
301 return __PACKAGE__->emake(@_);
307 # The number defined as pi = 180 degrees
309 sub pi () { 4 * CORE::atan2(1, 1) }
316 sub pi2 () { 2 * pi }
321 # The full circle twice.
323 sub pi4 () { 4 * pi }
330 sub pip2 () { pi / 2 }
337 sub pip4 () { pi / 4 }
344 sub _uplog10 () { 1 / CORE::log(10) }
349 # The number defined as i*i = -1;
354 $i->{'cartesian'} = [0, 1];
355 $i->{'polar'} = [1, pip2];
366 sub _ip2 () { i / 2 }
369 # Attribute access/set routines
372 sub _cartesian {$_[0]->{c_dirty} ?
373 $_[0]->_update_cartesian : $_[0]->{'cartesian'}}
374 sub _polar {$_[0]->{p_dirty} ?
375 $_[0]->_update_polar : $_[0]->{'polar'}}
377 sub _set_cartesian { $_[0]->{p_dirty}++; $_[0]->{c_dirty} = 0;
378 $_[0]->{'cartesian'} = $_[1] }
379 sub _set_polar { $_[0]->{c_dirty}++; $_[0]->{p_dirty} = 0;
380 $_[0]->{'polar'} = $_[1] }
383 # ->_update_cartesian
385 # Recompute and return the cartesian form, given accurate polar form.
387 sub _update_cartesian {
389 my ($r, $t) = @{$self->{'polar'}};
390 $self->{c_dirty} = 0;
391 return $self->{'cartesian'} = [$r * CORE::cos($t), $r * CORE::sin($t)];
398 # Recompute and return the polar form, given accurate cartesian form.
402 my ($x, $y) = @{$self->{'cartesian'}};
403 $self->{p_dirty} = 0;
404 return $self->{'polar'} = [0, 0] if $x == 0 && $y == 0;
405 return $self->{'polar'} = [CORE::sqrt($x*$x + $y*$y),
406 CORE::atan2($y, $x)];
415 my ($z1, $z2, $regular) = @_;
416 my ($re1, $im1) = @{$z1->_cartesian};
417 $z2 = cplx($z2) unless ref $z2;
418 my ($re2, $im2) = ref $z2 ? @{$z2->_cartesian} : ($z2, 0);
419 unless (defined $regular) {
420 $z1->_set_cartesian([$re1 + $re2, $im1 + $im2]);
423 return (ref $z1)->make($re1 + $re2, $im1 + $im2);
432 my ($z1, $z2, $inverted) = @_;
433 my ($re1, $im1) = @{$z1->_cartesian};
434 $z2 = cplx($z2) unless ref $z2;
435 my ($re2, $im2) = @{$z2->_cartesian};
436 unless (defined $inverted) {
437 $z1->_set_cartesian([$re1 - $re2, $im1 - $im2]);
441 (ref $z1)->make($re2 - $re1, $im2 - $im1) :
442 (ref $z1)->make($re1 - $re2, $im1 - $im2);
452 my ($z1, $z2, $regular) = @_;
453 if ($z1->{p_dirty} == 0 and ref $z2 and $z2->{p_dirty} == 0) {
454 # if both polar better use polar to avoid rounding errors
455 my ($r1, $t1) = @{$z1->_polar};
456 my ($r2, $t2) = @{$z2->_polar};
458 if ($t > pi()) { $t -= pi2 }
459 elsif ($t <= -pi()) { $t += pi2 }
460 unless (defined $regular) {
461 $z1->_set_polar([$r1 * $r2, $t]);
464 return (ref $z1)->emake($r1 * $r2, $t);
466 my ($x1, $y1) = @{$z1->_cartesian};
468 my ($x2, $y2) = @{$z2->_cartesian};
469 return (ref $z1)->make($x1*$x2-$y1*$y2, $x1*$y2+$y1*$x2);
471 return (ref $z1)->make($x1*$z2, $y1*$z2);
479 # Die on division by zero.
482 my $mess = "$_[0]: Division by zero.\n";
485 $mess .= "(Because in the definition of $_[0], the divisor ";
486 $mess .= "$_[1] " unless ("$_[1]" eq '0');
492 $mess .= "Died at $up[1] line $up[2].\n";
503 my ($z1, $z2, $inverted) = @_;
504 if ($z1->{p_dirty} == 0 and ref $z2 and $z2->{p_dirty} == 0) {
505 # if both polar better use polar to avoid rounding errors
506 my ($r1, $t1) = @{$z1->_polar};
507 my ($r2, $t2) = @{$z2->_polar};
510 _divbyzero "$z2/0" if ($r1 == 0);
512 if ($t > pi()) { $t -= pi2 }
513 elsif ($t <= -pi()) { $t += pi2 }
514 return (ref $z1)->emake($r2 / $r1, $t);
516 _divbyzero "$z1/0" if ($r2 == 0);
518 if ($t > pi()) { $t -= pi2 }
519 elsif ($t <= -pi()) { $t += pi2 }
520 return (ref $z1)->emake($r1 / $r2, $t);
525 ($x2, $y2) = @{$z1->_cartesian};
526 $d = $x2*$x2 + $y2*$y2;
527 _divbyzero "$z2/0" if $d == 0;
528 return (ref $z1)->make(($x2*$z2)/$d, -($y2*$z2)/$d);
530 my ($x1, $y1) = @{$z1->_cartesian};
532 ($x2, $y2) = @{$z2->_cartesian};
533 $d = $x2*$x2 + $y2*$y2;
534 _divbyzero "$z1/0" if $d == 0;
535 my $u = ($x1*$x2 + $y1*$y2)/$d;
536 my $v = ($y1*$x2 - $x1*$y2)/$d;
537 return (ref $z1)->make($u, $v);
539 _divbyzero "$z1/0" if $z2 == 0;
540 return (ref $z1)->make($x1/$z2, $y1/$z2);
549 # Computes z1**z2 = exp(z2 * log z1)).
552 my ($z1, $z2, $inverted) = @_;
554 return 1 if $z1 == 0 || $z2 == 1;
555 return 0 if $z2 == 0 && Re($z1) > 0;
557 return 1 if $z2 == 0 || $z1 == 1;
558 return 0 if $z1 == 0 && Re($z2) > 0;
560 my $w = $inverted ? &exp($z1 * &log($z2))
561 : &exp($z2 * &log($z1));
562 # If both arguments cartesian, return cartesian, else polar.
563 return $z1->{c_dirty} == 0 &&
564 (not ref $z2 or $z2->{c_dirty} == 0) ?
565 cplx(@{$w->_cartesian}) : $w;
571 # Computes z1 <=> z2.
572 # Sorts on the real part first, then on the imaginary part. Thus 2-4i < 3+8i.
575 my ($z1, $z2, $inverted) = @_;
576 my ($re1, $im1) = ref $z1 ? @{$z1->_cartesian} : ($z1, 0);
577 my ($re2, $im2) = ref $z2 ? @{$z2->_cartesian} : ($z2, 0);
578 my $sgn = $inverted ? -1 : 1;
579 return $sgn * ($re1 <=> $re2) if $re1 != $re2;
580 return $sgn * ($im1 <=> $im2);
588 # (Required in addition to _spaceship() because of NaNs.)
590 my ($z1, $z2, $inverted) = @_;
591 my ($re1, $im1) = ref $z1 ? @{$z1->_cartesian} : ($z1, 0);
592 my ($re2, $im2) = ref $z2 ? @{$z2->_cartesian} : ($z2, 0);
593 return $re1 == $re2 && $im1 == $im2 ? 1 : 0;
604 my ($r, $t) = @{$z->_polar};
605 $t = ($t <= 0) ? $t + pi : $t - pi;
606 return (ref $z)->emake($r, $t);
608 my ($re, $im) = @{$z->_cartesian};
609 return (ref $z)->make(-$re, -$im);
615 # Compute complex's _conjugate.
620 my ($r, $t) = @{$z->_polar};
621 return (ref $z)->emake($r, -$t);
623 my ($re, $im) = @{$z->_cartesian};
624 return (ref $z)->make($re, -$im);
630 # Compute or set complex's norm (rho).
633 my ($z, $rho) = @_ ? @_ : $_;
638 return CORE::abs($z);
642 $z->{'polar'} = [ $rho, ${$z->_polar}[1] ];
647 return ${$z->_polar}[0];
654 if ($$theta > pi()) { $$theta -= pi2 }
655 elsif ($$theta <= -pi()) { $$theta += pi2 }
661 # Compute or set complex's argument (theta).
664 my ($z, $theta) = @_;
665 return $z unless ref $z;
666 if (defined $theta) {
668 $z->{'polar'} = [ ${$z->_polar}[0], $theta ];
672 $theta = ${$z->_polar}[1];
683 # It is quite tempting to use wantarray here so that in list context
684 # sqrt() would return the two solutions. This, however, would
687 # print "sqrt(z) = ", sqrt($z), "\n";
689 # The two values would be printed side by side without no intervening
690 # whitespace, quite confusing.
691 # Therefore if you want the two solutions use the root().
694 my ($z) = @_ ? $_[0] : $_;
695 my ($re, $im) = ref $z ? @{$z->_cartesian} : ($z, 0);
696 return $re < 0 ? cplx(0, CORE::sqrt(-$re)) : CORE::sqrt($re)
698 my ($r, $t) = @{$z->_polar};
699 return (ref $z)->emake(CORE::sqrt($r), $t/2);
705 # Compute cbrt(z) (cubic root).
707 # Why are we not returning three values? The same answer as for sqrt().
712 -CORE::exp(CORE::log(-$z)/3) :
713 ($z > 0 ? CORE::exp(CORE::log($z)/3): 0)
715 my ($r, $t) = @{$z->_polar};
717 return (ref $z)->emake(CORE::exp(CORE::log($r)/3), $t/3);
726 my $mess = "Root '$_[0]' illegal, root rank must be positive integer.\n";
730 $mess .= "Died at $up[1] line $up[2].\n";
738 # Computes all nth root for z, returning an array whose size is n.
739 # `n' must be a positive integer.
741 # The roots are given by (for k = 0..n-1):
743 # z^(1/n) = r^(1/n) (cos ((t+2 k pi)/n) + i sin ((t+2 k pi)/n))
746 my ($z, $n, $k) = @_;
747 _rootbad($n) if ($n < 1 or int($n) != $n);
748 my ($r, $t) = ref $z ?
749 @{$z->_polar} : (CORE::abs($z), $z >= 0 ? 0 : pi);
750 my $theta_inc = pi2 / $n;
751 my $rho = $r ** (1/$n);
752 my $cartesian = ref $z && $z->{c_dirty} == 0;
755 for (my $i = 0, my $theta = $t / $n;
757 $i++, $theta += $theta_inc) {
758 my $w = cplxe($rho, $theta);
759 # Yes, $cartesian is loop invariant.
760 push @root, $cartesian ? cplx(@{$w->_cartesian}) : $w;
764 my $w = cplxe($rho, $t / $n + $k * $theta_inc);
765 return $cartesian ? cplx(@{$w->_cartesian}) : $w;
772 # Return or set Re(z).
776 return $z unless ref $z;
778 $z->{'cartesian'} = [ $Re, ${$z->_cartesian}[1] ];
782 return ${$z->_cartesian}[0];
789 # Return or set Im(z).
793 return 0 unless ref $z;
795 $z->{'cartesian'} = [ ${$z->_cartesian}[0], $Im ];
799 return ${$z->_cartesian}[1];
806 # Return or set rho(w).
809 Math::Complex::abs(@_);
815 # Return or set theta(w).
818 Math::Complex::arg(@_);
827 my ($z) = @_ ? @_ : $_;
828 return CORE::exp($z) unless ref $z;
829 my ($x, $y) = @{$z->_cartesian};
830 return (ref $z)->emake(CORE::exp($x), $y);
836 # Die on logarithm of zero.
839 my $mess = "$_[0]: Logarithm of zero.\n";
842 $mess .= "(Because in the definition of $_[0], the argument ";
843 $mess .= "$_[1] " unless ($_[1] eq '0');
849 $mess .= "Died at $up[1] line $up[2].\n";
860 my ($z) = @_ ? @_ : $_;
862 _logofzero("log") if $z == 0;
863 return $z > 0 ? CORE::log($z) : cplx(CORE::log(-$z), pi);
865 my ($r, $t) = @{$z->_polar};
866 _logofzero("log") if $r == 0;
867 if ($t > pi()) { $t -= pi2 }
868 elsif ($t <= -pi()) { $t += pi2 }
869 return (ref $z)->make(CORE::log($r), $t);
877 sub ln { Math::Complex::log(@_) }
886 return Math::Complex::log($_[0]) * _uplog10;
892 # Compute logn(z,n) = log(z) / log(n)
896 $z = cplx($z, 0) unless ref $z;
897 my $logn = $LOGN{$n};
898 $logn = $LOGN{$n} = CORE::log($n) unless defined $logn; # Cache log(n)
899 return &log($z) / $logn;
905 # Compute cos(z) = (exp(iz) + exp(-iz))/2.
908 my ($z) = @_ ? @_ : $_;
909 return CORE::cos($z) unless ref $z;
910 my ($x, $y) = @{$z->_cartesian};
911 my $ey = CORE::exp($y);
912 my $sx = CORE::sin($x);
913 my $cx = CORE::cos($x);
914 my $ey_1 = $ey ? 1 / $ey : Inf();
915 return (ref $z)->make($cx * ($ey + $ey_1)/2,
916 $sx * ($ey_1 - $ey)/2);
922 # Compute sin(z) = (exp(iz) - exp(-iz))/2.
925 my ($z) = @_ ? @_ : $_;
926 return CORE::sin($z) unless ref $z;
927 my ($x, $y) = @{$z->_cartesian};
928 my $ey = CORE::exp($y);
929 my $sx = CORE::sin($x);
930 my $cx = CORE::cos($x);
931 my $ey_1 = $ey ? 1 / $ey : Inf();
932 return (ref $z)->make($sx * ($ey + $ey_1)/2,
933 $cx * ($ey - $ey_1)/2);
939 # Compute tan(z) = sin(z) / cos(z).
944 _divbyzero "tan($z)", "cos($z)" if $cz == 0;
945 return &sin($z) / $cz;
951 # Computes the secant sec(z) = 1 / cos(z).
956 _divbyzero "sec($z)", "cos($z)" if ($cz == 0);
963 # Computes the cosecant csc(z) = 1 / sin(z).
968 _divbyzero "csc($z)", "sin($z)" if ($sz == 0);
977 sub cosec { Math::Complex::csc(@_) }
982 # Computes cot(z) = cos(z) / sin(z).
987 _divbyzero "cot($z)", "sin($z)" if ($sz == 0);
988 return &cos($z) / $sz;
996 sub cotan { Math::Complex::cot(@_) }
1001 # Computes the arc cosine acos(z) = -i log(z + sqrt(z*z-1)).
1005 return CORE::atan2(CORE::sqrt(1-$z*$z), $z)
1006 if (! ref $z) && CORE::abs($z) <= 1;
1007 $z = cplx($z, 0) unless ref $z;
1008 my ($x, $y) = @{$z->_cartesian};
1009 return 0 if $x == 1 && $y == 0;
1010 my $t1 = CORE::sqrt(($x+1)*($x+1) + $y*$y);
1011 my $t2 = CORE::sqrt(($x-1)*($x-1) + $y*$y);
1012 my $alpha = ($t1 + $t2)/2;
1013 my $beta = ($t1 - $t2)/2;
1014 $alpha = 1 if $alpha < 1;
1015 if ($beta > 1) { $beta = 1 }
1016 elsif ($beta < -1) { $beta = -1 }
1017 my $u = CORE::atan2(CORE::sqrt(1-$beta*$beta), $beta);
1018 my $v = CORE::log($alpha + CORE::sqrt($alpha*$alpha-1));
1019 $v = -$v if $y > 0 || ($y == 0 && $x < -1);
1020 return (ref $z)->make($u, $v);
1026 # Computes the arc sine asin(z) = -i log(iz + sqrt(1-z*z)).
1030 return CORE::atan2($z, CORE::sqrt(1-$z*$z))
1031 if (! ref $z) && CORE::abs($z) <= 1;
1032 $z = cplx($z, 0) unless ref $z;
1033 my ($x, $y) = @{$z->_cartesian};
1034 return 0 if $x == 0 && $y == 0;
1035 my $t1 = CORE::sqrt(($x+1)*($x+1) + $y*$y);
1036 my $t2 = CORE::sqrt(($x-1)*($x-1) + $y*$y);
1037 my $alpha = ($t1 + $t2)/2;
1038 my $beta = ($t1 - $t2)/2;
1039 $alpha = 1 if $alpha < 1;
1040 if ($beta > 1) { $beta = 1 }
1041 elsif ($beta < -1) { $beta = -1 }
1042 my $u = CORE::atan2($beta, CORE::sqrt(1-$beta*$beta));
1043 my $v = -CORE::log($alpha + CORE::sqrt($alpha*$alpha-1));
1044 $v = -$v if $y > 0 || ($y == 0 && $x < -1);
1045 return (ref $z)->make($u, $v);
1051 # Computes the arc tangent atan(z) = i/2 log((i+z) / (i-z)).
1055 return CORE::atan2($z, 1) unless ref $z;
1056 my ($x, $y) = ref $z ? @{$z->_cartesian} : ($z, 0);
1057 return 0 if $x == 0 && $y == 0;
1058 _divbyzero "atan(i)" if ( $z == i);
1059 _logofzero "atan(-i)" if (-$z == i); # -i is a bad file test...
1060 my $log = &log((i + $z) / (i - $z));
1067 # Computes the arc secant asec(z) = acos(1 / z).
1071 _divbyzero "asec($z)", $z if ($z == 0);
1072 return acos(1 / $z);
1078 # Computes the arc cosecant acsc(z) = asin(1 / z).
1082 _divbyzero "acsc($z)", $z if ($z == 0);
1083 return asin(1 / $z);
1091 sub acosec { Math::Complex::acsc(@_) }
1096 # Computes the arc cotangent acot(z) = atan(1 / z)
1100 _divbyzero "acot(0)" if $z == 0;
1101 return ($z >= 0) ? CORE::atan2(1, $z) : CORE::atan2(-1, -$z)
1103 _divbyzero "acot(i)" if ($z - i == 0);
1104 _logofzero "acot(-i)" if ($z + i == 0);
1105 return atan(1 / $z);
1113 sub acotan { Math::Complex::acot(@_) }
1118 # Computes the hyperbolic cosine cosh(z) = (exp(z) + exp(-z))/2.
1124 $ex = CORE::exp($z);
1125 return $ex ? ($ex == $ExpInf ? Inf() : ($ex + 1/$ex)/2) : Inf();
1127 my ($x, $y) = @{$z->_cartesian};
1128 $ex = CORE::exp($x);
1129 my $ex_1 = $ex ? 1 / $ex : Inf();
1130 return (ref $z)->make(CORE::cos($y) * ($ex + $ex_1)/2,
1131 CORE::sin($y) * ($ex - $ex_1)/2);
1137 # Computes the hyperbolic sine sinh(z) = (exp(z) - exp(-z))/2.
1143 return 0 if $z == 0;
1144 $ex = CORE::exp($z);
1145 return $ex ? ($ex == $ExpInf ? Inf() : ($ex - 1/$ex)/2) : -Inf();
1147 my ($x, $y) = @{$z->_cartesian};
1148 my $cy = CORE::cos($y);
1149 my $sy = CORE::sin($y);
1150 $ex = CORE::exp($x);
1151 my $ex_1 = $ex ? 1 / $ex : Inf();
1152 return (ref $z)->make(CORE::cos($y) * ($ex - $ex_1)/2,
1153 CORE::sin($y) * ($ex + $ex_1)/2);
1159 # Computes the hyperbolic tangent tanh(z) = sinh(z) / cosh(z).
1164 _divbyzero "tanh($z)", "cosh($z)" if ($cz == 0);
1166 return 1 if $cz == $sz;
1167 return -1 if $cz == -$sz;
1174 # Computes the hyperbolic secant sech(z) = 1 / cosh(z).
1179 _divbyzero "sech($z)", "cosh($z)" if ($cz == 0);
1186 # Computes the hyperbolic cosecant csch(z) = 1 / sinh(z).
1191 _divbyzero "csch($z)", "sinh($z)" if ($sz == 0);
1200 sub cosech { Math::Complex::csch(@_) }
1205 # Computes the hyperbolic cotangent coth(z) = cosh(z) / sinh(z).
1210 _divbyzero "coth($z)", "sinh($z)" if $sz == 0;
1212 return 1 if $cz == $sz;
1213 return -1 if $cz == -$sz;
1222 sub cotanh { Math::Complex::coth(@_) }
1227 # Computes the area/inverse hyperbolic cosine acosh(z) = log(z + sqrt(z*z-1)).
1234 my ($re, $im) = @{$z->_cartesian};
1236 return CORE::log($re + CORE::sqrt($re*$re - 1))
1238 return cplx(0, CORE::atan2(CORE::sqrt(1 - $re*$re), $re))
1239 if CORE::abs($re) < 1;
1241 my $t = &sqrt($z * $z - 1) + $z;
1242 # Try Taylor if looking bad (this usually means that
1243 # $z was large negative, therefore the sqrt is really
1244 # close to abs(z), summing that with z...)
1245 $t = 1/(2 * $z) - 1/(8 * $z**3) + 1/(16 * $z**5) - 5/(128 * $z**7)
1248 $u->Im(-$u->Im) if $re < 0 && $im == 0;
1249 return $re < 0 ? -$u : $u;
1255 # Computes the area/inverse hyperbolic sine asinh(z) = log(z + sqrt(z*z+1))
1260 my $t = $z + CORE::sqrt($z*$z + 1);
1261 return CORE::log($t) if $t;
1263 my $t = &sqrt($z * $z + 1) + $z;
1264 # Try Taylor if looking bad (this usually means that
1265 # $z was large negative, therefore the sqrt is really
1266 # close to abs(z), summing that with z...)
1267 $t = 1/(2 * $z) - 1/(8 * $z**3) + 1/(16 * $z**5) - 5/(128 * $z**7)
1275 # Computes the area/inverse hyperbolic tangent atanh(z) = 1/2 log((1+z) / (1-z)).
1280 return CORE::log((1 + $z)/(1 - $z))/2 if CORE::abs($z) < 1;
1283 _divbyzero 'atanh(1)', "1 - $z" if (1 - $z == 0);
1284 _logofzero 'atanh(-1)' if (1 + $z == 0);
1285 return 0.5 * &log((1 + $z) / (1 - $z));
1291 # Computes the area/inverse hyperbolic secant asech(z) = acosh(1 / z).
1295 _divbyzero 'asech(0)', "$z" if ($z == 0);
1296 return acosh(1 / $z);
1302 # Computes the area/inverse hyperbolic cosecant acsch(z) = asinh(1 / z).
1306 _divbyzero 'acsch(0)', $z if ($z == 0);
1307 return asinh(1 / $z);
1313 # Alias for acosh().
1315 sub acosech { Math::Complex::acsch(@_) }
1320 # Computes the area/inverse hyperbolic cotangent acoth(z) = 1/2 log((1+z) / (z-1)).
1324 _divbyzero 'acoth(0)' if ($z == 0);
1326 return CORE::log(($z + 1)/($z - 1))/2 if CORE::abs($z) > 1;
1329 _divbyzero 'acoth(1)', "$z - 1" if ($z - 1 == 0);
1330 _logofzero 'acoth(-1)', "1 + $z" if (1 + $z == 0);
1331 return &log((1 + $z) / ($z - 1)) / 2;
1339 sub acotanh { Math::Complex::acoth(@_) }
1344 # Compute atan(z1/z2), minding the right quadrant.
1347 my ($z1, $z2, $inverted) = @_;
1348 my ($re1, $im1, $re2, $im2);
1350 ($re1, $im1) = ref $z2 ? @{$z2->_cartesian} : ($z2, 0);
1351 ($re2, $im2) = ref $z1 ? @{$z1->_cartesian} : ($z1, 0);
1353 ($re1, $im1) = ref $z1 ? @{$z1->_cartesian} : ($z1, 0);
1354 ($re2, $im2) = ref $z2 ? @{$z2->_cartesian} : ($z2, 0);
1357 # In MATLAB the imaginary parts are ignored.
1358 # warn "atan2: Imaginary parts ignored";
1359 # http://documents.wolfram.com/mathematica/functions/ArcTan
1360 # NOTE: Mathematica ArcTan[x,y] while atan2(y,x)
1361 my $s = $z1 * $z1 + $z2 * $z2;
1362 _divbyzero("atan2") if $s == 0;
1364 my $r = $z2 + $z1 * $i;
1365 return -$i * &log($r / &sqrt( $s ));
1367 return CORE::atan2($re1, $re2);
1374 # Set (get if no argument) the display format for all complex numbers that
1375 # don't happen to have overridden it via ->display_format
1377 # When called as an object method, this actually sets the display format for
1378 # the current object.
1380 # Valid object formats are 'c' and 'p' for cartesian and polar. The first
1381 # letter is used actually, so the type can be fully spelled out for clarity.
1383 sub display_format {
1385 my %display_format = %DISPLAY_FORMAT;
1387 if (ref $self) { # Called as an object method
1388 if (exists $self->{display_format}) {
1389 my %obj = %{$self->{display_format}};
1390 @display_format{keys %obj} = values %obj;
1394 $display_format{style} = shift;
1397 @display_format{keys %new} = values %new;
1400 if (ref $self) { # Called as an object method
1401 $self->{display_format} = { %display_format };
1404 %{$self->{display_format}} :
1405 $self->{display_format}->{style};
1408 # Called as a class method
1409 %DISPLAY_FORMAT = %display_format;
1413 $DISPLAY_FORMAT{style};
1419 # Show nicely formatted complex number under its cartesian or polar form,
1420 # depending on the current display format:
1422 # . If a specific display format has been recorded for this object, use it.
1423 # . Otherwise, use the generic current default for all complex numbers,
1424 # which is a package global variable.
1429 my $style = $z->display_format;
1431 $style = $DISPLAY_FORMAT{style} unless defined $style;
1433 return $z->_stringify_polar if $style =~ /^p/i;
1434 return $z->_stringify_cartesian;
1438 # ->_stringify_cartesian
1440 # Stringify as a cartesian representation 'a+bi'.
1442 sub _stringify_cartesian {
1444 my ($x, $y) = @{$z->_cartesian};
1447 my %format = $z->display_format;
1448 my $format = $format{format};
1451 if ($x =~ /^NaN[QS]?$/i) {
1454 if ($x =~ /^-?\Q$Inf\E$/oi) {
1457 $re = defined $format ? sprintf($format, $x) : $x;
1465 if ($y =~ /^(NaN[QS]?)$/i) {
1468 if ($y =~ /^-?\Q$Inf\E$/oi) {
1473 sprintf($format, $y) :
1474 ($y == 1 ? "" : ($y == -1 ? "-" : $y));
1487 } elsif ($y > 0 || $im =~ /^NaN[QS]?i$/i) {
1488 $str .= "+" if defined $re;
1491 } elsif (!defined $re) {
1500 # ->_stringify_polar
1502 # Stringify as a polar representation '[r,t]'.
1504 sub _stringify_polar {
1506 my ($r, $t) = @{$z->_polar};
1509 my %format = $z->display_format;
1510 my $format = $format{format};
1512 if ($t =~ /^NaN[QS]?$/i || $t =~ /^-?\Q$Inf\E$/oi) {
1514 } elsif ($t == pi) {
1516 } elsif ($r == 0 || $t == 0) {
1517 $theta = defined $format ? sprintf($format, $t) : $t;
1520 return "[$r,$theta]" if defined $theta;
1523 # Try to identify pi/n and friends.
1526 $t -= int(CORE::abs($t) / pi2) * pi2;
1528 if ($format{polar_pretty_print} && $t) {
1532 if ($b =~ /^-?\d+$/) {
1533 $b = $b < 0 ? "-" : "" if CORE::abs($b) == 1;
1534 $theta = "${b}pi/$a";
1540 if (defined $format) {
1541 $r = sprintf($format, $r);
1542 $theta = sprintf($format, $theta) unless defined $theta;
1544 $theta = $t unless defined $theta;
1547 return "[$r,$theta]";
1561 Math::Complex - complex numbers and associated mathematical functions
1567 $z = Math::Complex->make(5, 6);
1569 $j = cplxe(1, 2*pi/3);
1573 This package lets you create and manipulate complex numbers. By default,
1574 I<Perl> limits itself to real numbers, but an extra C<use> statement brings
1575 full complex support, along with a full set of mathematical functions
1576 typically associated with and/or extended to complex numbers.
1578 If you wonder what complex numbers are, they were invented to be able to solve
1579 the following equation:
1583 and by definition, the solution is noted I<i> (engineers use I<j> instead since
1584 I<i> usually denotes an intensity, but the name does not matter). The number
1585 I<i> is a pure I<imaginary> number.
1587 The arithmetics with pure imaginary numbers works just like you would expect
1588 it with real numbers... you just have to remember that
1594 5i + 7i = i * (5 + 7) = 12i
1595 4i - 3i = i * (4 - 3) = i
1600 Complex numbers are numbers that have both a real part and an imaginary
1601 part, and are usually noted:
1605 where C<a> is the I<real> part and C<b> is the I<imaginary> part. The
1606 arithmetic with complex numbers is straightforward. You have to
1607 keep track of the real and the imaginary parts, but otherwise the
1608 rules used for real numbers just apply:
1610 (4 + 3i) + (5 - 2i) = (4 + 5) + i(3 - 2) = 9 + i
1611 (2 + i) * (4 - i) = 2*4 + 4i -2i -i*i = 8 + 2i + 1 = 9 + 2i
1613 A graphical representation of complex numbers is possible in a plane
1614 (also called the I<complex plane>, but it's really a 2D plane).
1619 is the point whose coordinates are (a, b). Actually, it would
1620 be the vector originating from (0, 0) to (a, b). It follows that the addition
1621 of two complex numbers is a vectorial addition.
1623 Since there is a bijection between a point in the 2D plane and a complex
1624 number (i.e. the mapping is unique and reciprocal), a complex number
1625 can also be uniquely identified with polar coordinates:
1629 where C<rho> is the distance to the origin, and C<theta> the angle between
1630 the vector and the I<x> axis. There is a notation for this using the
1631 exponential form, which is:
1633 rho * exp(i * theta)
1635 where I<i> is the famous imaginary number introduced above. Conversion
1636 between this form and the cartesian form C<a + bi> is immediate:
1638 a = rho * cos(theta)
1639 b = rho * sin(theta)
1641 which is also expressed by this formula:
1643 z = rho * exp(i * theta) = rho * (cos theta + i * sin theta)
1645 In other words, it's the projection of the vector onto the I<x> and I<y>
1646 axes. Mathematicians call I<rho> the I<norm> or I<modulus> and I<theta>
1647 the I<argument> of the complex number. The I<norm> of C<z> is
1648 marked here as C<abs(z)>.
1650 The polar notation (also known as the trigonometric representation) is
1651 much more handy for performing multiplications and divisions of
1652 complex numbers, whilst the cartesian notation is better suited for
1653 additions and subtractions. Real numbers are on the I<x> axis, and
1654 therefore I<y> or I<theta> is zero or I<pi>.
1656 All the common operations that can be performed on a real number have
1657 been defined to work on complex numbers as well, and are merely
1658 I<extensions> of the operations defined on real numbers. This means
1659 they keep their natural meaning when there is no imaginary part, provided
1660 the number is within their definition set.
1662 For instance, the C<sqrt> routine which computes the square root of
1663 its argument is only defined for non-negative real numbers and yields a
1664 non-negative real number (it is an application from B<R+> to B<R+>).
1665 If we allow it to return a complex number, then it can be extended to
1666 negative real numbers to become an application from B<R> to B<C> (the
1667 set of complex numbers):
1669 sqrt(x) = x >= 0 ? sqrt(x) : sqrt(-x)*i
1671 It can also be extended to be an application from B<C> to B<C>,
1672 whilst its restriction to B<R> behaves as defined above by using
1673 the following definition:
1675 sqrt(z = [r,t]) = sqrt(r) * exp(i * t/2)
1677 Indeed, a negative real number can be noted C<[x,pi]> (the modulus
1678 I<x> is always non-negative, so C<[x,pi]> is really C<-x>, a negative
1679 number) and the above definition states that
1681 sqrt([x,pi]) = sqrt(x) * exp(i*pi/2) = [sqrt(x),pi/2] = sqrt(x)*i
1683 which is exactly what we had defined for negative real numbers above.
1684 The C<sqrt> returns only one of the solutions: if you want the both,
1685 use the C<root> function.
1687 All the common mathematical functions defined on real numbers that
1688 are extended to complex numbers share that same property of working
1689 I<as usual> when the imaginary part is zero (otherwise, it would not
1690 be called an extension, would it?).
1692 A I<new> operation possible on a complex number that is
1693 the identity for real numbers is called the I<conjugate>, and is noted
1694 with a horizontal bar above the number, or C<~z> here.
1701 z * ~z = (a + bi) * (a - bi) = a*a + b*b
1703 We saw that the norm of C<z> was noted C<abs(z)> and was defined as the
1704 distance to the origin, also known as:
1706 rho = abs(z) = sqrt(a*a + b*b)
1710 z * ~z = abs(z) ** 2
1712 If z is a pure real number (i.e. C<b == 0>), then the above yields:
1716 which is true (C<abs> has the regular meaning for real number, i.e. stands
1717 for the absolute value). This example explains why the norm of C<z> is
1718 noted C<abs(z)>: it extends the C<abs> function to complex numbers, yet
1719 is the regular C<abs> we know when the complex number actually has no
1720 imaginary part... This justifies I<a posteriori> our use of the C<abs>
1721 notation for the norm.
1725 Given the following notations:
1727 z1 = a + bi = r1 * exp(i * t1)
1728 z2 = c + di = r2 * exp(i * t2)
1729 z = <any complex or real number>
1731 the following (overloaded) operations are supported on complex numbers:
1733 z1 + z2 = (a + c) + i(b + d)
1734 z1 - z2 = (a - c) + i(b - d)
1735 z1 * z2 = (r1 * r2) * exp(i * (t1 + t2))
1736 z1 / z2 = (r1 / r2) * exp(i * (t1 - t2))
1737 z1 ** z2 = exp(z2 * log z1)
1739 abs(z) = r1 = sqrt(a*a + b*b)
1740 sqrt(z) = sqrt(r1) * exp(i * t/2)
1741 exp(z) = exp(a) * exp(i * b)
1742 log(z) = log(r1) + i*t
1743 sin(z) = 1/2i (exp(i * z1) - exp(-i * z))
1744 cos(z) = 1/2 (exp(i * z1) + exp(-i * z))
1745 atan2(y, x) = atan(y / x) # Minding the right quadrant, note the order.
1747 The definition used for complex arguments of atan2() is
1749 -i log((x + iy)/sqrt(x*x+y*y))
1751 Note that atan2(0, 0) is not well-defined.
1753 The following extra operations are supported on both real and complex
1761 cbrt(z) = z ** (1/3)
1762 log10(z) = log(z) / log(10)
1763 logn(z, n) = log(z) / log(n)
1765 tan(z) = sin(z) / cos(z)
1771 asin(z) = -i * log(i*z + sqrt(1-z*z))
1772 acos(z) = -i * log(z + i*sqrt(1-z*z))
1773 atan(z) = i/2 * log((i+z) / (i-z))
1775 acsc(z) = asin(1 / z)
1776 asec(z) = acos(1 / z)
1777 acot(z) = atan(1 / z) = -i/2 * log((i+z) / (z-i))
1779 sinh(z) = 1/2 (exp(z) - exp(-z))
1780 cosh(z) = 1/2 (exp(z) + exp(-z))
1781 tanh(z) = sinh(z) / cosh(z) = (exp(z) - exp(-z)) / (exp(z) + exp(-z))
1783 csch(z) = 1 / sinh(z)
1784 sech(z) = 1 / cosh(z)
1785 coth(z) = 1 / tanh(z)
1787 asinh(z) = log(z + sqrt(z*z+1))
1788 acosh(z) = log(z + sqrt(z*z-1))
1789 atanh(z) = 1/2 * log((1+z) / (1-z))
1791 acsch(z) = asinh(1 / z)
1792 asech(z) = acosh(1 / z)
1793 acoth(z) = atanh(1 / z) = 1/2 * log((1+z) / (z-1))
1795 I<arg>, I<abs>, I<log>, I<csc>, I<cot>, I<acsc>, I<acot>, I<csch>,
1796 I<coth>, I<acosech>, I<acotanh>, have aliases I<rho>, I<theta>, I<ln>,
1797 I<cosec>, I<cotan>, I<acosec>, I<acotan>, I<cosech>, I<cotanh>,
1798 I<acosech>, I<acotanh>, respectively. C<Re>, C<Im>, C<arg>, C<abs>,
1799 C<rho>, and C<theta> can be used also as mutators. The C<cbrt>
1800 returns only one of the solutions: if you want all three, use the
1803 The I<root> function is available to compute all the I<n>
1804 roots of some complex, where I<n> is a strictly positive integer.
1805 There are exactly I<n> such roots, returned as a list. Getting the
1806 number mathematicians call C<j> such that:
1810 is a simple matter of writing:
1812 $j = ((root(1, 3))[1];
1814 The I<k>th root for C<z = [r,t]> is given by:
1816 (root(z, n))[k] = r**(1/n) * exp(i * (t + 2*k*pi)/n)
1818 You can return the I<k>th root directly by C<root(z, n, k)>,
1819 indexing starting from I<zero> and ending at I<n - 1>.
1821 The I<spaceship> numeric comparison operator, E<lt>=E<gt>, is also
1822 defined. In order to ensure its restriction to real numbers is conform
1823 to what you would expect, the comparison is run on the real part of
1824 the complex number first, and imaginary parts are compared only when
1825 the real parts match.
1829 To create a complex number, use either:
1831 $z = Math::Complex->make(3, 4);
1834 if you know the cartesian form of the number, or
1838 if you like. To create a number using the polar form, use either:
1840 $z = Math::Complex->emake(5, pi/3);
1841 $x = cplxe(5, pi/3);
1843 instead. The first argument is the modulus, the second is the angle
1844 (in radians, the full circle is 2*pi). (Mnemonic: C<e> is used as a
1845 notation for complex numbers in the polar form).
1847 It is possible to write:
1849 $x = cplxe(-3, pi/4);
1851 but that will be silently converted into C<[3,-3pi/4]>, since the
1852 modulus must be non-negative (it represents the distance to the origin
1853 in the complex plane).
1855 It is also possible to have a complex number as either argument of the
1856 C<make>, C<emake>, C<cplx>, and C<cplxe>: the appropriate component of
1857 the argument will be used.
1862 The C<new>, C<make>, C<emake>, C<cplx>, and C<cplxe> will also
1863 understand a single (string) argument of the forms
1871 in which case the appropriate cartesian and exponential components
1872 will be parsed from the string and used to create new complex numbers.
1873 The imaginary component and the theta, respectively, will default to zero.
1875 The C<new>, C<make>, C<emake>, C<cplx>, and C<cplxe> will also
1876 understand the case of no arguments: this means plain zero or (0, 0).
1880 When printed, a complex number is usually shown under its cartesian
1881 style I<a+bi>, but there are legitimate cases where the polar style
1882 I<[r,t]> is more appropriate. The process of converting the complex
1883 number into a string that can be displayed is known as I<stringification>.
1885 By calling the class method C<Math::Complex::display_format> and
1886 supplying either C<"polar"> or C<"cartesian"> as an argument, you
1887 override the default display style, which is C<"cartesian">. Not
1888 supplying any argument returns the current settings.
1890 This default can be overridden on a per-number basis by calling the
1891 C<display_format> method instead. As before, not supplying any argument
1892 returns the current display style for this number. Otherwise whatever you
1893 specify will be the new display style for I<this> particular number.
1899 Math::Complex::display_format('polar');
1900 $j = (root(1, 3))[1];
1901 print "j = $j\n"; # Prints "j = [1,2pi/3]"
1902 $j->display_format('cartesian');
1903 print "j = $j\n"; # Prints "j = -0.5+0.866025403784439i"
1905 The polar style attempts to emphasize arguments like I<k*pi/n>
1906 (where I<n> is a positive integer and I<k> an integer within [-9, +9]),
1907 this is called I<polar pretty-printing>.
1909 For the reverse of stringifying, see the C<make> and C<emake>.
1911 =head2 CHANGED IN PERL 5.6
1913 The C<display_format> class method and the corresponding
1914 C<display_format> object method can now be called using
1915 a parameter hash instead of just a one parameter.
1917 The old display format style, which can have values C<"cartesian"> or
1918 C<"polar">, can be changed using the C<"style"> parameter.
1920 $j->display_format(style => "polar");
1922 The one parameter calling convention also still works.
1924 $j->display_format("polar");
1926 There are two new display parameters.
1928 The first one is C<"format">, which is a sprintf()-style format string
1929 to be used for both numeric parts of the complex number(s). The is
1930 somewhat system-dependent but most often it corresponds to C<"%.15g">.
1931 You can revert to the default by setting the C<format> to C<undef>.
1933 # the $j from the above example
1935 $j->display_format('format' => '%.5f');
1936 print "j = $j\n"; # Prints "j = -0.50000+0.86603i"
1937 $j->display_format('format' => undef);
1938 print "j = $j\n"; # Prints "j = -0.5+0.86603i"
1940 Notice that this affects also the return values of the
1941 C<display_format> methods: in list context the whole parameter hash
1942 will be returned, as opposed to only the style parameter value.
1943 This is a potential incompatibility with earlier versions if you
1944 have been calling the C<display_format> method in list context.
1946 The second new display parameter is C<"polar_pretty_print">, which can
1947 be set to true or false, the default being true. See the previous
1948 section for what this means.
1952 Thanks to overloading, the handling of arithmetics with complex numbers
1953 is simple and almost transparent.
1955 Here are some examples:
1959 $j = cplxe(1, 2*pi/3); # $j ** 3 == 1
1960 print "j = $j, j**3 = ", $j ** 3, "\n";
1961 print "1 + j + j**2 = ", 1 + $j + $j**2, "\n";
1963 $z = -16 + 0*i; # Force it to be a complex
1964 print "sqrt($z) = ", sqrt($z), "\n";
1966 $k = exp(i * 2*pi/3);
1967 print "$j - $k = ", $j - $k, "\n";
1969 $z->Re(3); # Re, Im, arg, abs,
1970 $j->arg(2); # (the last two aka rho, theta)
1971 # can be used also as mutators.
1977 The constant C<pi> and some handy multiples of it (pi2, pi4,
1978 and pip2 (pi/2) and pip4 (pi/4)) are also available if separately
1981 use Math::Complex ':pi';
1982 $third_of_circle = pi2 / 3;
1986 The floating point infinity can be exported as a subroutine Inf():
1988 use Math::Complex qw(Inf sinh);
1989 my $AlsoInf = Inf() + 42;
1990 my $AnotherInf = sinh(1e42);
1991 print "$AlsoInf is $AnotherInf\n" if $AlsoInf == $AnotherInf;
1993 Note that the stringified form of infinity varies between platforms:
1994 it can be for example any of
2001 or it can be something else.
2003 Also note that in some platforms trying to use the infinity in
2004 arithmetic operations may result in Perl crashing because using
2005 an infinity causes SIGFPE or its moral equivalent to be sent.
2006 The way to ignore this is
2008 local $SIG{FPE} = sub { };
2010 =head1 ERRORS DUE TO DIVISION BY ZERO OR LOGARITHM OF ZERO
2012 The division (/) and the following functions
2018 atanh asech acsch acoth
2020 cannot be computed for all arguments because that would mean dividing
2021 by zero or taking logarithm of zero. These situations cause fatal
2022 runtime errors looking like this
2024 cot(0): Division by zero.
2025 (Because in the definition of cot(0), the divisor sin(0) is 0)
2030 atanh(-1): Logarithm of zero.
2033 For the C<csc>, C<cot>, C<asec>, C<acsc>, C<acot>, C<csch>, C<coth>,
2034 C<asech>, C<acsch>, the argument cannot be C<0> (zero). For the
2035 logarithmic functions and the C<atanh>, C<acoth>, the argument cannot
2036 be C<1> (one). For the C<atanh>, C<acoth>, the argument cannot be
2037 C<-1> (minus one). For the C<atan>, C<acot>, the argument cannot be
2038 C<i> (the imaginary unit). For the C<atan>, C<acoth>, the argument
2039 cannot be C<-i> (the negative imaginary unit). For the C<tan>,
2040 C<sec>, C<tanh>, the argument cannot be I<pi/2 + k * pi>, where I<k>
2041 is any integer. atan2(0, 0) is undefined, and if the complex arguments
2042 are used for atan2(), a division by zero will happen if z1**2+z2**2 == 0.
2044 Note that because we are operating on approximations of real numbers,
2045 these errors can happen when merely `too close' to the singularities
2048 =head1 ERRORS DUE TO INDIGESTIBLE ARGUMENTS
2050 The C<make> and C<emake> accept both real and complex arguments.
2051 When they cannot recognize the arguments they will die with error
2052 messages like the following
2054 Math::Complex::make: Cannot take real part of ...
2055 Math::Complex::make: Cannot take real part of ...
2056 Math::Complex::emake: Cannot take rho of ...
2057 Math::Complex::emake: Cannot take theta of ...
2061 Saying C<use Math::Complex;> exports many mathematical routines in the
2062 caller environment and even overrides some (C<sqrt>, C<log>, C<atan2>).
2063 This is construed as a feature by the Authors, actually... ;-)
2065 All routines expect to be given real or complex numbers. Don't attempt to
2066 use BigFloat, since Perl has currently no rule to disambiguate a '+'
2067 operation (for instance) between two overloaded entities.
2069 In Cray UNICOS there is some strange numerical instability that results
2070 in root(), cos(), sin(), cosh(), sinh(), losing accuracy fast. Beware.
2071 The bug may be in UNICOS math libs, in UNICOS C compiler, in Math::Complex.
2072 Whatever it is, it does not manifest itself anywhere else where Perl runs.
2080 Daniel S. Lewart <F<lewart!at!uiuc.edu>>
2081 Jarkko Hietaniemi <F<jhi!at!iki.fi>>
2082 Raphael Manfredi <F<Raphael_Manfredi!at!pobox.com>>
2086 This library is free software; you can redistribute it and/or modify
2087 it under the same terms as Perl itself.