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
10 use vars qw($VERSION @ISA @EXPORT @EXPORT_OK %EXPORT_TAGS $Inf);
15 unless ($^O eq 'unicosmk') {
17 # We do want an arithmetic overflow, Inf INF inf Infinity:.
18 undef $Inf unless eval <<'EOE' and $Inf =~ /^inf(?:inity)?$/i;
19 local $SIG{FPE} = sub {die};
23 if (!defined $Inf) { # Try a different method
24 undef $Inf unless eval <<'EOE' and $Inf =~ /^inf(?:inity)?$/i;
25 local $SIG{FPE} = sub {die};
27 $Inf = $t + "1e99999999999999999999999999999999";
31 $Inf = "Inf" if !defined $Inf || !($Inf > 0); # Desperation.
39 # Regular expression for floating point numbers.
40 # These days we could use Scalar::Util::lln(), I guess.
41 my $gre = qr'\s*([\+\-]?(?:(?:(?:\d+(?:_\d+)*(?:\.\d*(?:_\d+)*)?|\.\d+(?:_\d+)*)(?:[eE][\+\-]?\d+(?:_\d+)*)?))|inf)'i;
50 csc cosec sec cot cotan
52 acsc acosec asec acot acotan
54 csch cosech sech coth cotanh
56 acsch acosech asech acoth acotanh
68 my @pi = qw(pi pi2 pi4 pip2 pip4);
84 '<=>' => \&_spaceship,
101 my %DISPLAY_FORMAT = ('style' => 'cartesian',
102 'polar_pretty_print' => 1);
103 my $eps = 1e-14; # Epsilon
106 # Object attributes (internal):
107 # cartesian [real, imaginary] -- cartesian form
108 # polar [rho, theta] -- polar form
109 # c_dirty cartesian form not up-to-date
110 # p_dirty polar form not up-to-date
111 # display display format (package's global when not set)
116 # Die on bad *make() arguments.
119 die "@{[(caller(1))[3]]}: Cannot take $_[0] of '$_[1]'.\n";
126 if ($arg =~ /^$gre$/) {
128 } elsif ($arg =~ /^(?:$gre)?$gre\s*i\s*$/) {
129 ($p, $q) = ($1 || 0, $2);
130 } elsif ($arg =~ /^\s*\(\s*$gre\s*(?:,\s*$gre\s*)?\)\s*$/) {
131 ($p, $q) = ($1, $2 || 0);
136 $p =~ s/^(-?)inf$/"${1}9**9**9"/e;
138 $q =~ s/^(-?)inf$/"${1}9**9**9"/e;
148 if ($arg =~ /^\s*\[\s*$gre\s*(?:,\s*$gre\s*)?\]\s*$/) {
149 ($p, $q) = ($1, $2 || 0);
150 } elsif ($arg =~ m!^\s*\[\s*$gre\s*(?:,\s*([-+]?\d*\s*)?pi(?:/\s*(\d+))?\s*)?\]\s*$!) {
151 ($p, $q) = ($1, ($2 eq '-' ? -1 : ($2 || 1)) * pi() / ($3 || 1));
152 } elsif ($arg =~ /^\s*\[\s*$gre\s*\]\s*$/) {
154 } elsif ($arg =~ /^\s*$gre\s*$/) {
161 $p =~ s/^(-?)inf$/"${1}9**9**9"/e;
162 $q =~ s/^(-?)inf$/"${1}9**9**9"/e;
171 # Create a new complex number (cartesian form)
174 my $self = bless {}, shift;
179 return (ref $self)->emake($_[0])
180 if ($_[0] =~ /^\s*\[/);
181 ($re, $im) = _make($_[0]);
186 _cannot_make("real part", $re) unless $re =~ /^$gre$/;
189 _cannot_make("imaginary part", $im) unless $im =~ /^$gre$/;
190 $self->_set_cartesian([$re, $im ]);
191 $self->display_format('cartesian');
199 # Create a new complex number (exponential form)
202 my $self = bless {}, shift;
205 ($rho, $theta) = (0, 0);
207 return (ref $self)->make($_[0])
208 if ($_[0] =~ /^\s*\(/ || $_[0] =~ /i\s*$/);
209 ($rho, $theta) = _emake($_[0]);
213 if (defined $rho && defined $theta) {
216 $theta = ($theta <= 0) ? $theta + pi() : $theta - pi();
220 _cannot_make("rho", $rho) unless $rho =~ /^$gre$/;
223 _cannot_make("theta", $theta) unless $theta =~ /^$gre$/;
224 $self->_set_polar([$rho, $theta]);
225 $self->display_format('polar');
230 sub new { &make } # For backward compatibility only.
235 # Creates a complex number from a (re, im) tuple.
236 # This avoids the burden of writing Math::Complex->make(re, im).
239 return __PACKAGE__->make(@_);
245 # Creates a complex number from a (rho, theta) tuple.
246 # This avoids the burden of writing Math::Complex->emake(rho, theta).
249 return __PACKAGE__->emake(@_);
255 # The number defined as pi = 180 degrees
257 sub pi () { 4 * CORE::atan2(1, 1) }
264 sub pi2 () { 2 * pi }
269 # The full circle twice.
271 sub pi4 () { 4 * pi }
278 sub pip2 () { pi / 2 }
285 sub pip4 () { pi / 4 }
292 sub _uplog10 () { 1 / CORE::log(10) }
297 # The number defined as i*i = -1;
302 $i->{'cartesian'} = [0, 1];
303 $i->{'polar'} = [1, pip2];
314 sub _ip2 () { i / 2 }
317 # Attribute access/set routines
320 sub _cartesian {$_[0]->{c_dirty} ?
321 $_[0]->_update_cartesian : $_[0]->{'cartesian'}}
322 sub _polar {$_[0]->{p_dirty} ?
323 $_[0]->_update_polar : $_[0]->{'polar'}}
325 sub _set_cartesian { $_[0]->{p_dirty}++; $_[0]->{c_dirty} = 0;
326 $_[0]->{'cartesian'} = $_[1] }
327 sub _set_polar { $_[0]->{c_dirty}++; $_[0]->{p_dirty} = 0;
328 $_[0]->{'polar'} = $_[1] }
331 # ->_update_cartesian
333 # Recompute and return the cartesian form, given accurate polar form.
335 sub _update_cartesian {
337 my ($r, $t) = @{$self->{'polar'}};
338 $self->{c_dirty} = 0;
339 return $self->{'cartesian'} = [$r * CORE::cos($t), $r * CORE::sin($t)];
346 # Recompute and return the polar form, given accurate cartesian form.
350 my ($x, $y) = @{$self->{'cartesian'}};
351 $self->{p_dirty} = 0;
352 return $self->{'polar'} = [0, 0] if $x == 0 && $y == 0;
353 return $self->{'polar'} = [CORE::sqrt($x*$x + $y*$y),
354 CORE::atan2($y, $x)];
363 my ($z1, $z2, $regular) = @_;
364 my ($re1, $im1) = @{$z1->_cartesian};
365 $z2 = cplx($z2) unless ref $z2;
366 my ($re2, $im2) = ref $z2 ? @{$z2->_cartesian} : ($z2, 0);
367 unless (defined $regular) {
368 $z1->_set_cartesian([$re1 + $re2, $im1 + $im2]);
371 return (ref $z1)->make($re1 + $re2, $im1 + $im2);
380 my ($z1, $z2, $inverted) = @_;
381 my ($re1, $im1) = @{$z1->_cartesian};
382 $z2 = cplx($z2) unless ref $z2;
383 my ($re2, $im2) = @{$z2->_cartesian};
384 unless (defined $inverted) {
385 $z1->_set_cartesian([$re1 - $re2, $im1 - $im2]);
389 (ref $z1)->make($re2 - $re1, $im2 - $im1) :
390 (ref $z1)->make($re1 - $re2, $im1 - $im2);
400 my ($z1, $z2, $regular) = @_;
401 if ($z1->{p_dirty} == 0 and ref $z2 and $z2->{p_dirty} == 0) {
402 # if both polar better use polar to avoid rounding errors
403 my ($r1, $t1) = @{$z1->_polar};
404 my ($r2, $t2) = @{$z2->_polar};
406 if ($t > pi()) { $t -= pi2 }
407 elsif ($t <= -pi()) { $t += pi2 }
408 unless (defined $regular) {
409 $z1->_set_polar([$r1 * $r2, $t]);
412 return (ref $z1)->emake($r1 * $r2, $t);
414 my ($x1, $y1) = @{$z1->_cartesian};
416 my ($x2, $y2) = @{$z2->_cartesian};
417 return (ref $z1)->make($x1*$x2-$y1*$y2, $x1*$y2+$y1*$x2);
419 return (ref $z1)->make($x1*$z2, $y1*$z2);
427 # Die on division by zero.
430 my $mess = "$_[0]: Division by zero.\n";
433 $mess .= "(Because in the definition of $_[0], the divisor ";
434 $mess .= "$_[1] " unless ("$_[1]" eq '0');
440 $mess .= "Died at $up[1] line $up[2].\n";
451 my ($z1, $z2, $inverted) = @_;
452 if ($z1->{p_dirty} == 0 and ref $z2 and $z2->{p_dirty} == 0) {
453 # if both polar better use polar to avoid rounding errors
454 my ($r1, $t1) = @{$z1->_polar};
455 my ($r2, $t2) = @{$z2->_polar};
458 _divbyzero "$z2/0" if ($r1 == 0);
460 if ($t > pi()) { $t -= pi2 }
461 elsif ($t <= -pi()) { $t += pi2 }
462 return (ref $z1)->emake($r2 / $r1, $t);
464 _divbyzero "$z1/0" if ($r2 == 0);
466 if ($t > pi()) { $t -= pi2 }
467 elsif ($t <= -pi()) { $t += pi2 }
468 return (ref $z1)->emake($r1 / $r2, $t);
473 ($x2, $y2) = @{$z1->_cartesian};
474 $d = $x2*$x2 + $y2*$y2;
475 _divbyzero "$z2/0" if $d == 0;
476 return (ref $z1)->make(($x2*$z2)/$d, -($y2*$z2)/$d);
478 my ($x1, $y1) = @{$z1->_cartesian};
480 ($x2, $y2) = @{$z2->_cartesian};
481 $d = $x2*$x2 + $y2*$y2;
482 _divbyzero "$z1/0" if $d == 0;
483 my $u = ($x1*$x2 + $y1*$y2)/$d;
484 my $v = ($y1*$x2 - $x1*$y2)/$d;
485 return (ref $z1)->make($u, $v);
487 _divbyzero "$z1/0" if $z2 == 0;
488 return (ref $z1)->make($x1/$z2, $y1/$z2);
497 # Computes z1**z2 = exp(z2 * log z1)).
500 my ($z1, $z2, $inverted) = @_;
502 return 1 if $z1 == 0 || $z2 == 1;
503 return 0 if $z2 == 0 && Re($z1) > 0;
505 return 1 if $z2 == 0 || $z1 == 1;
506 return 0 if $z1 == 0 && Re($z2) > 0;
508 my $w = $inverted ? &exp($z1 * &log($z2))
509 : &exp($z2 * &log($z1));
510 # If both arguments cartesian, return cartesian, else polar.
511 return $z1->{c_dirty} == 0 &&
512 (not ref $z2 or $z2->{c_dirty} == 0) ?
513 cplx(@{$w->_cartesian}) : $w;
519 # Computes z1 <=> z2.
520 # Sorts on the real part first, then on the imaginary part. Thus 2-4i < 3+8i.
523 my ($z1, $z2, $inverted) = @_;
524 my ($re1, $im1) = ref $z1 ? @{$z1->_cartesian} : ($z1, 0);
525 my ($re2, $im2) = ref $z2 ? @{$z2->_cartesian} : ($z2, 0);
526 my $sgn = $inverted ? -1 : 1;
527 return $sgn * ($re1 <=> $re2) if $re1 != $re2;
528 return $sgn * ($im1 <=> $im2);
536 # (Required in addition to _spaceship() because of NaNs.)
538 my ($z1, $z2, $inverted) = @_;
539 my ($re1, $im1) = ref $z1 ? @{$z1->_cartesian} : ($z1, 0);
540 my ($re2, $im2) = ref $z2 ? @{$z2->_cartesian} : ($z2, 0);
541 return $re1 == $re2 && $im1 == $im2 ? 1 : 0;
552 my ($r, $t) = @{$z->_polar};
553 $t = ($t <= 0) ? $t + pi : $t - pi;
554 return (ref $z)->emake($r, $t);
556 my ($re, $im) = @{$z->_cartesian};
557 return (ref $z)->make(-$re, -$im);
563 # Compute complex's _conjugate.
568 my ($r, $t) = @{$z->_polar};
569 return (ref $z)->emake($r, -$t);
571 my ($re, $im) = @{$z->_cartesian};
572 return (ref $z)->make($re, -$im);
578 # Compute or set complex's norm (rho).
586 return CORE::abs($z);
590 $z->{'polar'} = [ $rho, ${$z->_polar}[1] ];
595 return ${$z->_polar}[0];
602 if ($$theta > pi()) { $$theta -= pi2 }
603 elsif ($$theta <= -pi()) { $$theta += pi2 }
609 # Compute or set complex's argument (theta).
612 my ($z, $theta) = @_;
613 return $z unless ref $z;
614 if (defined $theta) {
616 $z->{'polar'} = [ ${$z->_polar}[0], $theta ];
620 $theta = ${$z->_polar}[1];
631 # It is quite tempting to use wantarray here so that in list context
632 # sqrt() would return the two solutions. This, however, would
635 # print "sqrt(z) = ", sqrt($z), "\n";
637 # The two values would be printed side by side without no intervening
638 # whitespace, quite confusing.
639 # Therefore if you want the two solutions use the root().
643 my ($re, $im) = ref $z ? @{$z->_cartesian} : ($z, 0);
644 return $re < 0 ? cplx(0, CORE::sqrt(-$re)) : CORE::sqrt($re)
646 my ($r, $t) = @{$z->_polar};
647 return (ref $z)->emake(CORE::sqrt($r), $t/2);
653 # Compute cbrt(z) (cubic root).
655 # Why are we not returning three values? The same answer as for sqrt().
660 -CORE::exp(CORE::log(-$z)/3) :
661 ($z > 0 ? CORE::exp(CORE::log($z)/3): 0)
663 my ($r, $t) = @{$z->_polar};
665 return (ref $z)->emake(CORE::exp(CORE::log($r)/3), $t/3);
674 my $mess = "Root '$_[0]' illegal, root rank must be positive integer.\n";
678 $mess .= "Died at $up[1] line $up[2].\n";
686 # Computes all nth root for z, returning an array whose size is n.
687 # `n' must be a positive integer.
689 # The roots are given by (for k = 0..n-1):
691 # z^(1/n) = r^(1/n) (cos ((t+2 k pi)/n) + i sin ((t+2 k pi)/n))
694 my ($z, $n, $k) = @_;
695 _rootbad($n) if ($n < 1 or int($n) != $n);
696 my ($r, $t) = ref $z ?
697 @{$z->_polar} : (CORE::abs($z), $z >= 0 ? 0 : pi);
698 my $theta_inc = pi2 / $n;
699 my $rho = $r ** (1/$n);
700 my $cartesian = ref $z && $z->{c_dirty} == 0;
703 for (my $i = 0, my $theta = $t / $n;
705 $i++, $theta += $theta_inc) {
706 my $w = cplxe($rho, $theta);
707 # Yes, $cartesian is loop invariant.
708 push @root, $cartesian ? cplx(@{$w->_cartesian}) : $w;
712 my $w = cplxe($rho, $t / $n + $k * $theta_inc);
713 return $cartesian ? cplx(@{$w->_cartesian}) : $w;
720 # Return or set Re(z).
724 return $z unless ref $z;
726 $z->{'cartesian'} = [ $Re, ${$z->_cartesian}[1] ];
730 return ${$z->_cartesian}[0];
737 # Return or set Im(z).
741 return 0 unless ref $z;
743 $z->{'cartesian'} = [ ${$z->_cartesian}[0], $Im ];
747 return ${$z->_cartesian}[1];
754 # Return or set rho(w).
757 Math::Complex::abs(@_);
763 # Return or set theta(w).
766 Math::Complex::arg(@_);
776 my ($x, $y) = @{$z->_cartesian};
777 return (ref $z)->emake(CORE::exp($x), $y);
783 # Die on logarithm of zero.
786 my $mess = "$_[0]: Logarithm of zero.\n";
789 $mess .= "(Because in the definition of $_[0], the argument ";
790 $mess .= "$_[1] " unless ($_[1] eq '0');
796 $mess .= "Died at $up[1] line $up[2].\n";
809 _logofzero("log") if $z == 0;
810 return $z > 0 ? CORE::log($z) : cplx(CORE::log(-$z), pi);
812 my ($r, $t) = @{$z->_polar};
813 _logofzero("log") if $r == 0;
814 if ($t > pi()) { $t -= pi2 }
815 elsif ($t <= -pi()) { $t += pi2 }
816 return (ref $z)->make(CORE::log($r), $t);
824 sub ln { Math::Complex::log(@_) }
833 return Math::Complex::log($_[0]) * _uplog10;
839 # Compute logn(z,n) = log(z) / log(n)
843 $z = cplx($z, 0) unless ref $z;
844 my $logn = $LOGN{$n};
845 $logn = $LOGN{$n} = CORE::log($n) unless defined $logn; # Cache log(n)
846 return &log($z) / $logn;
852 # Compute cos(z) = (exp(iz) + exp(-iz))/2.
856 return CORE::cos($z) unless ref $z;
857 my ($x, $y) = @{$z->_cartesian};
858 my $ey = CORE::exp($y);
859 my $sx = CORE::sin($x);
860 my $cx = CORE::cos($x);
861 my $ey_1 = $ey ? 1 / $ey : $Inf;
862 return (ref $z)->make($cx * ($ey + $ey_1)/2,
863 $sx * ($ey_1 - $ey)/2);
869 # Compute sin(z) = (exp(iz) - exp(-iz))/2.
873 return CORE::sin($z) unless ref $z;
874 my ($x, $y) = @{$z->_cartesian};
875 my $ey = CORE::exp($y);
876 my $sx = CORE::sin($x);
877 my $cx = CORE::cos($x);
878 my $ey_1 = $ey ? 1 / $ey : $Inf;
879 return (ref $z)->make($sx * ($ey + $ey_1)/2,
880 $cx * ($ey - $ey_1)/2);
886 # Compute tan(z) = sin(z) / cos(z).
891 _divbyzero "tan($z)", "cos($z)" if $cz == 0;
892 return &sin($z) / $cz;
898 # Computes the secant sec(z) = 1 / cos(z).
903 _divbyzero "sec($z)", "cos($z)" if ($cz == 0);
910 # Computes the cosecant csc(z) = 1 / sin(z).
915 _divbyzero "csc($z)", "sin($z)" if ($sz == 0);
924 sub cosec { Math::Complex::csc(@_) }
929 # Computes cot(z) = cos(z) / sin(z).
934 _divbyzero "cot($z)", "sin($z)" if ($sz == 0);
935 return &cos($z) / $sz;
943 sub cotan { Math::Complex::cot(@_) }
948 # Computes the arc cosine acos(z) = -i log(z + sqrt(z*z-1)).
952 return CORE::atan2(CORE::sqrt(1-$z*$z), $z)
953 if (! ref $z) && CORE::abs($z) <= 1;
954 $z = cplx($z, 0) unless ref $z;
955 my ($x, $y) = @{$z->_cartesian};
956 return 0 if $x == 1 && $y == 0;
957 my $t1 = CORE::sqrt(($x+1)*($x+1) + $y*$y);
958 my $t2 = CORE::sqrt(($x-1)*($x-1) + $y*$y);
959 my $alpha = ($t1 + $t2)/2;
960 my $beta = ($t1 - $t2)/2;
961 $alpha = 1 if $alpha < 1;
962 if ($beta > 1) { $beta = 1 }
963 elsif ($beta < -1) { $beta = -1 }
964 my $u = CORE::atan2(CORE::sqrt(1-$beta*$beta), $beta);
965 my $v = CORE::log($alpha + CORE::sqrt($alpha*$alpha-1));
966 $v = -$v if $y > 0 || ($y == 0 && $x < -1);
967 return (ref $z)->make($u, $v);
973 # Computes the arc sine asin(z) = -i log(iz + sqrt(1-z*z)).
977 return CORE::atan2($z, CORE::sqrt(1-$z*$z))
978 if (! ref $z) && CORE::abs($z) <= 1;
979 $z = cplx($z, 0) unless ref $z;
980 my ($x, $y) = @{$z->_cartesian};
981 return 0 if $x == 0 && $y == 0;
982 my $t1 = CORE::sqrt(($x+1)*($x+1) + $y*$y);
983 my $t2 = CORE::sqrt(($x-1)*($x-1) + $y*$y);
984 my $alpha = ($t1 + $t2)/2;
985 my $beta = ($t1 - $t2)/2;
986 $alpha = 1 if $alpha < 1;
987 if ($beta > 1) { $beta = 1 }
988 elsif ($beta < -1) { $beta = -1 }
989 my $u = CORE::atan2($beta, CORE::sqrt(1-$beta*$beta));
990 my $v = -CORE::log($alpha + CORE::sqrt($alpha*$alpha-1));
991 $v = -$v if $y > 0 || ($y == 0 && $x < -1);
992 return (ref $z)->make($u, $v);
998 # Computes the arc tangent atan(z) = i/2 log((i+z) / (i-z)).
1002 return CORE::atan2($z, 1) unless ref $z;
1003 my ($x, $y) = ref $z ? @{$z->_cartesian} : ($z, 0);
1004 return 0 if $x == 0 && $y == 0;
1005 _divbyzero "atan(i)" if ( $z == i);
1006 _logofzero "atan(-i)" if (-$z == i); # -i is a bad file test...
1007 my $log = &log((i + $z) / (i - $z));
1014 # Computes the arc secant asec(z) = acos(1 / z).
1018 _divbyzero "asec($z)", $z if ($z == 0);
1019 return acos(1 / $z);
1025 # Computes the arc cosecant acsc(z) = asin(1 / z).
1029 _divbyzero "acsc($z)", $z if ($z == 0);
1030 return asin(1 / $z);
1038 sub acosec { Math::Complex::acsc(@_) }
1043 # Computes the arc cotangent acot(z) = atan(1 / z)
1047 _divbyzero "acot(0)" if $z == 0;
1048 return ($z >= 0) ? CORE::atan2(1, $z) : CORE::atan2(-1, -$z)
1050 _divbyzero "acot(i)" if ($z - i == 0);
1051 _logofzero "acot(-i)" if ($z + i == 0);
1052 return atan(1 / $z);
1060 sub acotan { Math::Complex::acot(@_) }
1065 # Computes the hyperbolic cosine cosh(z) = (exp(z) + exp(-z))/2.
1071 $ex = CORE::exp($z);
1072 return $ex ? ($ex + 1/$ex)/2 : $Inf;
1074 my ($x, $y) = @{$z->_cartesian};
1075 $ex = CORE::exp($x);
1076 my $ex_1 = $ex ? 1 / $ex : $Inf;
1077 return (ref $z)->make(CORE::cos($y) * ($ex + $ex_1)/2,
1078 CORE::sin($y) * ($ex - $ex_1)/2);
1084 # Computes the hyperbolic sine sinh(z) = (exp(z) - exp(-z))/2.
1090 return 0 if $z == 0;
1091 $ex = CORE::exp($z);
1092 return $ex ? ($ex - 1/$ex)/2 : "-$Inf";
1094 my ($x, $y) = @{$z->_cartesian};
1095 my $cy = CORE::cos($y);
1096 my $sy = CORE::sin($y);
1097 $ex = CORE::exp($x);
1098 my $ex_1 = $ex ? 1 / $ex : $Inf;
1099 return (ref $z)->make(CORE::cos($y) * ($ex - $ex_1)/2,
1100 CORE::sin($y) * ($ex + $ex_1)/2);
1106 # Computes the hyperbolic tangent tanh(z) = sinh(z) / cosh(z).
1111 _divbyzero "tanh($z)", "cosh($z)" if ($cz == 0);
1112 return sinh($z) / $cz;
1118 # Computes the hyperbolic secant sech(z) = 1 / cosh(z).
1123 _divbyzero "sech($z)", "cosh($z)" if ($cz == 0);
1130 # Computes the hyperbolic cosecant csch(z) = 1 / sinh(z).
1135 _divbyzero "csch($z)", "sinh($z)" if ($sz == 0);
1144 sub cosech { Math::Complex::csch(@_) }
1149 # Computes the hyperbolic cotangent coth(z) = cosh(z) / sinh(z).
1154 _divbyzero "coth($z)", "sinh($z)" if $sz == 0;
1155 return cosh($z) / $sz;
1163 sub cotanh { Math::Complex::coth(@_) }
1168 # Computes the arc hyperbolic cosine acosh(z) = log(z + sqrt(z*z-1)).
1175 my ($re, $im) = @{$z->_cartesian};
1177 return CORE::log($re + CORE::sqrt($re*$re - 1))
1179 return cplx(0, CORE::atan2(CORE::sqrt(1 - $re*$re), $re))
1180 if CORE::abs($re) < 1;
1182 my $t = &sqrt($z * $z - 1) + $z;
1183 # Try Taylor if looking bad (this usually means that
1184 # $z was large negative, therefore the sqrt is really
1185 # close to abs(z), summing that with z...)
1186 $t = 1/(2 * $z) - 1/(8 * $z**3) + 1/(16 * $z**5) - 5/(128 * $z**7)
1189 $u->Im(-$u->Im) if $re < 0 && $im == 0;
1190 return $re < 0 ? -$u : $u;
1196 # Computes the arc hyperbolic sine asinh(z) = log(z + sqrt(z*z+1))
1201 my $t = $z + CORE::sqrt($z*$z + 1);
1202 return CORE::log($t) if $t;
1204 my $t = &sqrt($z * $z + 1) + $z;
1205 # Try Taylor if looking bad (this usually means that
1206 # $z was large negative, therefore the sqrt is really
1207 # close to abs(z), summing that with z...)
1208 $t = 1/(2 * $z) - 1/(8 * $z**3) + 1/(16 * $z**5) - 5/(128 * $z**7)
1216 # Computes the arc hyperbolic tangent atanh(z) = 1/2 log((1+z) / (1-z)).
1221 return CORE::log((1 + $z)/(1 - $z))/2 if CORE::abs($z) < 1;
1224 _divbyzero 'atanh(1)', "1 - $z" if (1 - $z == 0);
1225 _logofzero 'atanh(-1)' if (1 + $z == 0);
1226 return 0.5 * &log((1 + $z) / (1 - $z));
1232 # Computes the hyperbolic arc secant asech(z) = acosh(1 / z).
1236 _divbyzero 'asech(0)', "$z" if ($z == 0);
1237 return acosh(1 / $z);
1243 # Computes the hyperbolic arc cosecant acsch(z) = asinh(1 / z).
1247 _divbyzero 'acsch(0)', $z if ($z == 0);
1248 return asinh(1 / $z);
1254 # Alias for acosh().
1256 sub acosech { Math::Complex::acsch(@_) }
1261 # Computes the arc hyperbolic cotangent acoth(z) = 1/2 log((1+z) / (z-1)).
1265 _divbyzero 'acoth(0)' if ($z == 0);
1267 return CORE::log(($z + 1)/($z - 1))/2 if CORE::abs($z) > 1;
1270 _divbyzero 'acoth(1)', "$z - 1" if ($z - 1 == 0);
1271 _logofzero 'acoth(-1)', "1 + $z" if (1 + $z == 0);
1272 return &log((1 + $z) / ($z - 1)) / 2;
1280 sub acotanh { Math::Complex::acoth(@_) }
1285 # Compute atan(z1/z2), minding the right quadrant.
1288 my ($z1, $z2, $inverted) = @_;
1289 my ($re1, $im1, $re2, $im2);
1291 ($re1, $im1) = ref $z2 ? @{$z2->_cartesian} : ($z2, 0);
1292 ($re2, $im2) = ref $z1 ? @{$z1->_cartesian} : ($z1, 0);
1294 ($re1, $im1) = ref $z1 ? @{$z1->_cartesian} : ($z1, 0);
1295 ($re2, $im2) = ref $z2 ? @{$z2->_cartesian} : ($z2, 0);
1298 # In MATLAB the imaginary parts are ignored.
1299 # warn "atan2: Imaginary parts ignored";
1300 # http://documents.wolfram.com/mathematica/functions/ArcTan
1301 # NOTE: Mathematica ArcTan[x,y] while atan2(y,x)
1302 my $s = $z1 * $z1 + $z2 * $z2;
1303 _divbyzero("atan2") if $s == 0;
1305 my $r = $z2 + $z1 * $i;
1306 return -$i * &log($r / &sqrt( $s ));
1308 return CORE::atan2($re1, $re2);
1315 # Set (get if no argument) the display format for all complex numbers that
1316 # don't happen to have overridden it via ->display_format
1318 # When called as an object method, this actually sets the display format for
1319 # the current object.
1321 # Valid object formats are 'c' and 'p' for cartesian and polar. The first
1322 # letter is used actually, so the type can be fully spelled out for clarity.
1324 sub display_format {
1326 my %display_format = %DISPLAY_FORMAT;
1328 if (ref $self) { # Called as an object method
1329 if (exists $self->{display_format}) {
1330 my %obj = %{$self->{display_format}};
1331 @display_format{keys %obj} = values %obj;
1335 $display_format{style} = shift;
1338 @display_format{keys %new} = values %new;
1341 if (ref $self) { # Called as an object method
1342 $self->{display_format} = { %display_format };
1345 %{$self->{display_format}} :
1346 $self->{display_format}->{style};
1349 # Called as a class method
1350 %DISPLAY_FORMAT = %display_format;
1354 $DISPLAY_FORMAT{style};
1360 # Show nicely formatted complex number under its cartesian or polar form,
1361 # depending on the current display format:
1363 # . If a specific display format has been recorded for this object, use it.
1364 # . Otherwise, use the generic current default for all complex numbers,
1365 # which is a package global variable.
1370 my $style = $z->display_format;
1372 $style = $DISPLAY_FORMAT{style} unless defined $style;
1374 return $z->_stringify_polar if $style =~ /^p/i;
1375 return $z->_stringify_cartesian;
1379 # ->_stringify_cartesian
1381 # Stringify as a cartesian representation 'a+bi'.
1383 sub _stringify_cartesian {
1385 my ($x, $y) = @{$z->_cartesian};
1388 my %format = $z->display_format;
1389 my $format = $format{format};
1392 if ($x =~ /^NaN[QS]?$/i) {
1395 if ($x =~ /^-?$Inf$/oi) {
1398 $re = defined $format ? sprintf($format, $x) : $x;
1406 if ($y =~ /^(NaN[QS]?)$/i) {
1409 if ($y =~ /^-?$Inf$/oi) {
1414 sprintf($format, $y) :
1415 ($y == 1 ? "" : ($y == -1 ? "-" : $y));
1428 } elsif ($y > 0 || $im =~ /^NaN[QS]?i$/i) {
1429 $str .= "+" if defined $re;
1432 } elsif (!defined $re) {
1441 # ->_stringify_polar
1443 # Stringify as a polar representation '[r,t]'.
1445 sub _stringify_polar {
1447 my ($r, $t) = @{$z->_polar};
1450 my %format = $z->display_format;
1451 my $format = $format{format};
1453 if ($t =~ /^NaN[QS]?$/i || $t =~ /^-?$Inf$/oi) {
1455 } elsif ($t == pi) {
1457 } elsif ($r == 0 || $t == 0) {
1458 $theta = defined $format ? sprintf($format, $t) : $t;
1461 return "[$r,$theta]" if defined $theta;
1464 # Try to identify pi/n and friends.
1467 $t -= int(CORE::abs($t) / pi2) * pi2;
1469 if ($format{polar_pretty_print} && $t) {
1473 if ($b =~ /^-?\d+$/) {
1474 $b = $b < 0 ? "-" : "" if CORE::abs($b) == 1;
1475 $theta = "${b}pi/$a";
1481 if (defined $format) {
1482 $r = sprintf($format, $r);
1483 $theta = sprintf($format, $theta) unless defined $theta;
1485 $theta = $t unless defined $theta;
1488 return "[$r,$theta]";
1498 Math::Complex - complex numbers and associated mathematical functions
1504 $z = Math::Complex->make(5, 6);
1506 $j = cplxe(1, 2*pi/3);
1510 This package lets you create and manipulate complex numbers. By default,
1511 I<Perl> limits itself to real numbers, but an extra C<use> statement brings
1512 full complex support, along with a full set of mathematical functions
1513 typically associated with and/or extended to complex numbers.
1515 If you wonder what complex numbers are, they were invented to be able to solve
1516 the following equation:
1520 and by definition, the solution is noted I<i> (engineers use I<j> instead since
1521 I<i> usually denotes an intensity, but the name does not matter). The number
1522 I<i> is a pure I<imaginary> number.
1524 The arithmetics with pure imaginary numbers works just like you would expect
1525 it with real numbers... you just have to remember that
1531 5i + 7i = i * (5 + 7) = 12i
1532 4i - 3i = i * (4 - 3) = i
1537 Complex numbers are numbers that have both a real part and an imaginary
1538 part, and are usually noted:
1542 where C<a> is the I<real> part and C<b> is the I<imaginary> part. The
1543 arithmetic with complex numbers is straightforward. You have to
1544 keep track of the real and the imaginary parts, but otherwise the
1545 rules used for real numbers just apply:
1547 (4 + 3i) + (5 - 2i) = (4 + 5) + i(3 - 2) = 9 + i
1548 (2 + i) * (4 - i) = 2*4 + 4i -2i -i*i = 8 + 2i + 1 = 9 + 2i
1550 A graphical representation of complex numbers is possible in a plane
1551 (also called the I<complex plane>, but it's really a 2D plane).
1556 is the point whose coordinates are (a, b). Actually, it would
1557 be the vector originating from (0, 0) to (a, b). It follows that the addition
1558 of two complex numbers is a vectorial addition.
1560 Since there is a bijection between a point in the 2D plane and a complex
1561 number (i.e. the mapping is unique and reciprocal), a complex number
1562 can also be uniquely identified with polar coordinates:
1566 where C<rho> is the distance to the origin, and C<theta> the angle between
1567 the vector and the I<x> axis. There is a notation for this using the
1568 exponential form, which is:
1570 rho * exp(i * theta)
1572 where I<i> is the famous imaginary number introduced above. Conversion
1573 between this form and the cartesian form C<a + bi> is immediate:
1575 a = rho * cos(theta)
1576 b = rho * sin(theta)
1578 which is also expressed by this formula:
1580 z = rho * exp(i * theta) = rho * (cos theta + i * sin theta)
1582 In other words, it's the projection of the vector onto the I<x> and I<y>
1583 axes. Mathematicians call I<rho> the I<norm> or I<modulus> and I<theta>
1584 the I<argument> of the complex number. The I<norm> of C<z> is
1585 marked here as C<abs(z)>.
1587 The polar notation (also known as the trigonometric representation) is
1588 much more handy for performing multiplications and divisions of
1589 complex numbers, whilst the cartesian notation is better suited for
1590 additions and subtractions. Real numbers are on the I<x> axis, and
1591 therefore I<y> or I<theta> is zero or I<pi>.
1593 All the common operations that can be performed on a real number have
1594 been defined to work on complex numbers as well, and are merely
1595 I<extensions> of the operations defined on real numbers. This means
1596 they keep their natural meaning when there is no imaginary part, provided
1597 the number is within their definition set.
1599 For instance, the C<sqrt> routine which computes the square root of
1600 its argument is only defined for non-negative real numbers and yields a
1601 non-negative real number (it is an application from B<R+> to B<R+>).
1602 If we allow it to return a complex number, then it can be extended to
1603 negative real numbers to become an application from B<R> to B<C> (the
1604 set of complex numbers):
1606 sqrt(x) = x >= 0 ? sqrt(x) : sqrt(-x)*i
1608 It can also be extended to be an application from B<C> to B<C>,
1609 whilst its restriction to B<R> behaves as defined above by using
1610 the following definition:
1612 sqrt(z = [r,t]) = sqrt(r) * exp(i * t/2)
1614 Indeed, a negative real number can be noted C<[x,pi]> (the modulus
1615 I<x> is always non-negative, so C<[x,pi]> is really C<-x>, a negative
1616 number) and the above definition states that
1618 sqrt([x,pi]) = sqrt(x) * exp(i*pi/2) = [sqrt(x),pi/2] = sqrt(x)*i
1620 which is exactly what we had defined for negative real numbers above.
1621 The C<sqrt> returns only one of the solutions: if you want the both,
1622 use the C<root> function.
1624 All the common mathematical functions defined on real numbers that
1625 are extended to complex numbers share that same property of working
1626 I<as usual> when the imaginary part is zero (otherwise, it would not
1627 be called an extension, would it?).
1629 A I<new> operation possible on a complex number that is
1630 the identity for real numbers is called the I<conjugate>, and is noted
1631 with a horizontal bar above the number, or C<~z> here.
1638 z * ~z = (a + bi) * (a - bi) = a*a + b*b
1640 We saw that the norm of C<z> was noted C<abs(z)> and was defined as the
1641 distance to the origin, also known as:
1643 rho = abs(z) = sqrt(a*a + b*b)
1647 z * ~z = abs(z) ** 2
1649 If z is a pure real number (i.e. C<b == 0>), then the above yields:
1653 which is true (C<abs> has the regular meaning for real number, i.e. stands
1654 for the absolute value). This example explains why the norm of C<z> is
1655 noted C<abs(z)>: it extends the C<abs> function to complex numbers, yet
1656 is the regular C<abs> we know when the complex number actually has no
1657 imaginary part... This justifies I<a posteriori> our use of the C<abs>
1658 notation for the norm.
1662 Given the following notations:
1664 z1 = a + bi = r1 * exp(i * t1)
1665 z2 = c + di = r2 * exp(i * t2)
1666 z = <any complex or real number>
1668 the following (overloaded) operations are supported on complex numbers:
1670 z1 + z2 = (a + c) + i(b + d)
1671 z1 - z2 = (a - c) + i(b - d)
1672 z1 * z2 = (r1 * r2) * exp(i * (t1 + t2))
1673 z1 / z2 = (r1 / r2) * exp(i * (t1 - t2))
1674 z1 ** z2 = exp(z2 * log z1)
1676 abs(z) = r1 = sqrt(a*a + b*b)
1677 sqrt(z) = sqrt(r1) * exp(i * t/2)
1678 exp(z) = exp(a) * exp(i * b)
1679 log(z) = log(r1) + i*t
1680 sin(z) = 1/2i (exp(i * z1) - exp(-i * z))
1681 cos(z) = 1/2 (exp(i * z1) + exp(-i * z))
1682 atan2(y, x) = atan(y / x) # Minding the right quadrant, note the order.
1684 The definition used for complex arguments of atan2() is
1686 -i log((x + iy)/sqrt(x*x+y*y))
1688 Note that atan2(0, 0) is not well-defined.
1690 The following extra operations are supported on both real and complex
1698 cbrt(z) = z ** (1/3)
1699 log10(z) = log(z) / log(10)
1700 logn(z, n) = log(z) / log(n)
1702 tan(z) = sin(z) / cos(z)
1708 asin(z) = -i * log(i*z + sqrt(1-z*z))
1709 acos(z) = -i * log(z + i*sqrt(1-z*z))
1710 atan(z) = i/2 * log((i+z) / (i-z))
1712 acsc(z) = asin(1 / z)
1713 asec(z) = acos(1 / z)
1714 acot(z) = atan(1 / z) = -i/2 * log((i+z) / (z-i))
1716 sinh(z) = 1/2 (exp(z) - exp(-z))
1717 cosh(z) = 1/2 (exp(z) + exp(-z))
1718 tanh(z) = sinh(z) / cosh(z) = (exp(z) - exp(-z)) / (exp(z) + exp(-z))
1720 csch(z) = 1 / sinh(z)
1721 sech(z) = 1 / cosh(z)
1722 coth(z) = 1 / tanh(z)
1724 asinh(z) = log(z + sqrt(z*z+1))
1725 acosh(z) = log(z + sqrt(z*z-1))
1726 atanh(z) = 1/2 * log((1+z) / (1-z))
1728 acsch(z) = asinh(1 / z)
1729 asech(z) = acosh(1 / z)
1730 acoth(z) = atanh(1 / z) = 1/2 * log((1+z) / (z-1))
1732 I<arg>, I<abs>, I<log>, I<csc>, I<cot>, I<acsc>, I<acot>, I<csch>,
1733 I<coth>, I<acosech>, I<acotanh>, have aliases I<rho>, I<theta>, I<ln>,
1734 I<cosec>, I<cotan>, I<acosec>, I<acotan>, I<cosech>, I<cotanh>,
1735 I<acosech>, I<acotanh>, respectively. C<Re>, C<Im>, C<arg>, C<abs>,
1736 C<rho>, and C<theta> can be used also as mutators. The C<cbrt>
1737 returns only one of the solutions: if you want all three, use the
1740 The I<root> function is available to compute all the I<n>
1741 roots of some complex, where I<n> is a strictly positive integer.
1742 There are exactly I<n> such roots, returned as a list. Getting the
1743 number mathematicians call C<j> such that:
1747 is a simple matter of writing:
1749 $j = ((root(1, 3))[1];
1751 The I<k>th root for C<z = [r,t]> is given by:
1753 (root(z, n))[k] = r**(1/n) * exp(i * (t + 2*k*pi)/n)
1755 You can return the I<k>th root directly by C<root(z, n, k)>,
1756 indexing starting from I<zero> and ending at I<n - 1>.
1758 The I<spaceship> comparison operator, E<lt>=E<gt>, is also defined. In
1759 order to ensure its restriction to real numbers is conform to what you
1760 would expect, the comparison is run on the real part of the complex
1761 number first, and imaginary parts are compared only when the real
1766 To create a complex number, use either:
1768 $z = Math::Complex->make(3, 4);
1771 if you know the cartesian form of the number, or
1775 if you like. To create a number using the polar form, use either:
1777 $z = Math::Complex->emake(5, pi/3);
1778 $x = cplxe(5, pi/3);
1780 instead. The first argument is the modulus, the second is the angle
1781 (in radians, the full circle is 2*pi). (Mnemonic: C<e> is used as a
1782 notation for complex numbers in the polar form).
1784 It is possible to write:
1786 $x = cplxe(-3, pi/4);
1788 but that will be silently converted into C<[3,-3pi/4]>, since the
1789 modulus must be non-negative (it represents the distance to the origin
1790 in the complex plane).
1792 It is also possible to have a complex number as either argument of the
1793 C<make>, C<emake>, C<cplx>, and C<cplxe>: the appropriate component of
1794 the argument will be used.
1799 The C<new>, C<make>, C<emake>, C<cplx>, and C<cplxe> will also
1800 understand a single (string) argument of the forms
1808 in which case the appropriate cartesian and exponential components
1809 will be parsed from the string and used to create new complex numbers.
1810 The imaginary component and the theta, respectively, will default to zero.
1812 The C<new>, C<make>, C<emake>, C<cplx>, and C<cplxe> will also
1813 understand the case of no arguments: this means plain zero or (0, 0).
1817 When printed, a complex number is usually shown under its cartesian
1818 style I<a+bi>, but there are legitimate cases where the polar style
1819 I<[r,t]> is more appropriate. The process of converting the complex
1820 number into a string that can be displayed is known as I<stringification>.
1822 By calling the class method C<Math::Complex::display_format> and
1823 supplying either C<"polar"> or C<"cartesian"> as an argument, you
1824 override the default display style, which is C<"cartesian">. Not
1825 supplying any argument returns the current settings.
1827 This default can be overridden on a per-number basis by calling the
1828 C<display_format> method instead. As before, not supplying any argument
1829 returns the current display style for this number. Otherwise whatever you
1830 specify will be the new display style for I<this> particular number.
1836 Math::Complex::display_format('polar');
1837 $j = (root(1, 3))[1];
1838 print "j = $j\n"; # Prints "j = [1,2pi/3]"
1839 $j->display_format('cartesian');
1840 print "j = $j\n"; # Prints "j = -0.5+0.866025403784439i"
1842 The polar style attempts to emphasize arguments like I<k*pi/n>
1843 (where I<n> is a positive integer and I<k> an integer within [-9, +9]),
1844 this is called I<polar pretty-printing>.
1846 For the reverse of stringifying, see the C<make> and C<emake>.
1848 =head2 CHANGED IN PERL 5.6
1850 The C<display_format> class method and the corresponding
1851 C<display_format> object method can now be called using
1852 a parameter hash instead of just a one parameter.
1854 The old display format style, which can have values C<"cartesian"> or
1855 C<"polar">, can be changed using the C<"style"> parameter.
1857 $j->display_format(style => "polar");
1859 The one parameter calling convention also still works.
1861 $j->display_format("polar");
1863 There are two new display parameters.
1865 The first one is C<"format">, which is a sprintf()-style format string
1866 to be used for both numeric parts of the complex number(s). The is
1867 somewhat system-dependent but most often it corresponds to C<"%.15g">.
1868 You can revert to the default by setting the C<format> to C<undef>.
1870 # the $j from the above example
1872 $j->display_format('format' => '%.5f');
1873 print "j = $j\n"; # Prints "j = -0.50000+0.86603i"
1874 $j->display_format('format' => undef);
1875 print "j = $j\n"; # Prints "j = -0.5+0.86603i"
1877 Notice that this affects also the return values of the
1878 C<display_format> methods: in list context the whole parameter hash
1879 will be returned, as opposed to only the style parameter value.
1880 This is a potential incompatibility with earlier versions if you
1881 have been calling the C<display_format> method in list context.
1883 The second new display parameter is C<"polar_pretty_print">, which can
1884 be set to true or false, the default being true. See the previous
1885 section for what this means.
1889 Thanks to overloading, the handling of arithmetics with complex numbers
1890 is simple and almost transparent.
1892 Here are some examples:
1896 $j = cplxe(1, 2*pi/3); # $j ** 3 == 1
1897 print "j = $j, j**3 = ", $j ** 3, "\n";
1898 print "1 + j + j**2 = ", 1 + $j + $j**2, "\n";
1900 $z = -16 + 0*i; # Force it to be a complex
1901 print "sqrt($z) = ", sqrt($z), "\n";
1903 $k = exp(i * 2*pi/3);
1904 print "$j - $k = ", $j - $k, "\n";
1906 $z->Re(3); # Re, Im, arg, abs,
1907 $j->arg(2); # (the last two aka rho, theta)
1908 # can be used also as mutators.
1912 The constant C<pi> and some handy multiples of it (pi2, pi4,
1913 and pip2 (pi/2) and pip4 (pi/4)) are also available if separately
1916 use Math::Complex ':pi';
1917 $third_of_circle = pi2 / 3;
1919 =head1 ERRORS DUE TO DIVISION BY ZERO OR LOGARITHM OF ZERO
1921 The division (/) and the following functions
1927 atanh asech acsch acoth
1929 cannot be computed for all arguments because that would mean dividing
1930 by zero or taking logarithm of zero. These situations cause fatal
1931 runtime errors looking like this
1933 cot(0): Division by zero.
1934 (Because in the definition of cot(0), the divisor sin(0) is 0)
1939 atanh(-1): Logarithm of zero.
1942 For the C<csc>, C<cot>, C<asec>, C<acsc>, C<acot>, C<csch>, C<coth>,
1943 C<asech>, C<acsch>, the argument cannot be C<0> (zero). For the
1944 logarithmic functions and the C<atanh>, C<acoth>, the argument cannot
1945 be C<1> (one). For the C<atanh>, C<acoth>, the argument cannot be
1946 C<-1> (minus one). For the C<atan>, C<acot>, the argument cannot be
1947 C<i> (the imaginary unit). For the C<atan>, C<acoth>, the argument
1948 cannot be C<-i> (the negative imaginary unit). For the C<tan>,
1949 C<sec>, C<tanh>, the argument cannot be I<pi/2 + k * pi>, where I<k>
1950 is any integer. atan2(0, 0) is undefined, and if the complex arguments
1951 are used for atan2(), a division by zero will happen if z1**2+z2**2 == 0.
1953 Note that because we are operating on approximations of real numbers,
1954 these errors can happen when merely `too close' to the singularities
1957 =head1 ERRORS DUE TO INDIGESTIBLE ARGUMENTS
1959 The C<make> and C<emake> accept both real and complex arguments.
1960 When they cannot recognize the arguments they will die with error
1961 messages like the following
1963 Math::Complex::make: Cannot take real part of ...
1964 Math::Complex::make: Cannot take real part of ...
1965 Math::Complex::emake: Cannot take rho of ...
1966 Math::Complex::emake: Cannot take theta of ...
1970 Saying C<use Math::Complex;> exports many mathematical routines in the
1971 caller environment and even overrides some (C<sqrt>, C<log>, C<atan2>).
1972 This is construed as a feature by the Authors, actually... ;-)
1974 All routines expect to be given real or complex numbers. Don't attempt to
1975 use BigFloat, since Perl has currently no rule to disambiguate a '+'
1976 operation (for instance) between two overloaded entities.
1978 In Cray UNICOS there is some strange numerical instability that results
1979 in root(), cos(), sin(), cosh(), sinh(), losing accuracy fast. Beware.
1980 The bug may be in UNICOS math libs, in UNICOS C compiler, in Math::Complex.
1981 Whatever it is, it does not manifest itself anywhere else where Perl runs.
1985 Daniel S. Lewart <F<lewart!at!uiuc.edu>>
1986 Jarkko Hietaniemi <F<jhi!at!iki.fi>>
1987 Raphael Manfredi <F<Raphael_Manfredi!at!pobox.com>>