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 our($VERSION, @ISA, @EXPORT, %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";
30 $! = $e; # Clear ERANGE.
32 $Inf = "Inf" if !defined $Inf || !($Inf > 0); # Desperation.
40 # Regular expression for floating point numbers.
41 my $gre = qr'\s*([\+\-]?(?:(?:(?:\d+(?:_\d+)*(?:\.\d*(?:_\d+)*)?|\.\d+(?:_\d+)*)(?:[eE][\+\-]?\d+(?:_\d+)*)?)))';
50 csc cosec sec cot cotan
52 acsc acosec asec acot acotan
54 csch cosech sech coth cotanh
56 acsch acosech asech acoth acotanh
95 my %DISPLAY_FORMAT = ('style' => 'cartesian',
96 'polar_pretty_print' => 1);
97 my $eps = 1e-14; # Epsilon
100 # Object attributes (internal):
101 # cartesian [real, imaginary] -- cartesian form
102 # polar [rho, theta] -- polar form
103 # c_dirty cartesian form not up-to-date
104 # p_dirty polar form not up-to-date
105 # display display format (package's global when not set)
108 # Die on bad *make() arguments.
111 die "@{[(caller(1))[3]]}: Cannot take $_[0] of $_[1].\n";
118 if ($arg =~ /^(?:$gre)?$gre\s*i\s*$/) {
119 ($p, $q) = ($1 || 0, $2);
121 } elsif ($arg =~ /^\s*\[\s*$gre\s*(?:,\s*$gre\s*)?\]\s*$/) {
122 ($p, $q) = ($1, $2 || 0);
131 return ($made, $p, $q);
137 # Create a new complex number (cartesian form)
140 my $self = bless {}, shift;
143 my ($remade, $p, $q) = _remake($re);
145 if ($remade eq 'cart') {
146 ($re, $im) = ($p, $q);
148 return (ref $self)->emake($p, $q);
154 if ( $rre eq ref $self ) {
157 _cannot_make("real part", $rre);
162 if ( $rim eq ref $self ) {
165 _cannot_make("imaginary part", $rim);
168 _cannot_make("real part", $re) unless $re =~ /^$gre$/;
170 _cannot_make("imaginary part", $im) unless $im =~ /^$gre$/;
171 $self->{'cartesian'} = [ $re, $im ];
172 $self->{c_dirty} = 0;
173 $self->{p_dirty} = 1;
174 $self->display_format('cartesian');
181 # Create a new complex number (exponential form)
184 my $self = bless {}, shift;
185 my ($rho, $theta) = @_;
187 my ($remade, $p, $q) = _remake($rho);
189 if ($remade eq 'exp') {
190 ($rho, $theta) = ($p, $q);
192 return (ref $self)->make($p, $q);
198 if ( $rrh eq ref $self ) {
201 _cannot_make("rho", $rrh);
204 my $rth = ref $theta;
206 if ( $rth eq ref $self ) {
207 $theta = theta($theta);
209 _cannot_make("theta", $rth);
214 $theta = ($theta <= 0) ? $theta + pi() : $theta - pi();
216 _cannot_make("rho", $rho) unless $rho =~ /^$gre$/;
218 _cannot_make("theta", $theta) unless $theta =~ /^$gre$/;
219 $self->{'polar'} = [$rho, $theta];
220 $self->{p_dirty} = 0;
221 $self->{c_dirty} = 1;
222 $self->display_format('polar');
226 sub new { &make } # For backward compatibility only.
231 # Creates a complex number from a (re, im) tuple.
232 # This avoids the burden of writing Math::Complex->make(re, im).
235 return __PACKAGE__->make(@_);
241 # Creates a complex number from a (rho, theta) tuple.
242 # This avoids the burden of writing Math::Complex->emake(rho, theta).
245 return __PACKAGE__->emake(@_);
251 # The number defined as pi = 180 degrees
253 sub pi () { 4 * CORE::atan2(1, 1) }
260 sub pit2 () { 2 * pi }
267 sub pip2 () { pi / 2 }
272 # One degree in radians, used in stringify_polar.
275 sub deg1 () { pi / 180 }
282 sub uplog10 () { 1 / CORE::log(10) }
287 # The number defined as i*i = -1;
292 $i->{'cartesian'} = [0, 1];
293 $i->{'polar'} = [1, pip2];
307 # Attribute access/set routines
310 sub cartesian {$_[0]->{c_dirty} ?
311 $_[0]->update_cartesian : $_[0]->{'cartesian'}}
312 sub polar {$_[0]->{p_dirty} ?
313 $_[0]->update_polar : $_[0]->{'polar'}}
315 sub set_cartesian { $_[0]->{p_dirty}++; $_[0]->{'cartesian'} = $_[1] }
316 sub set_polar { $_[0]->{c_dirty}++; $_[0]->{'polar'} = $_[1] }
321 # Recompute and return the cartesian form, given accurate polar form.
323 sub update_cartesian {
325 my ($r, $t) = @{$self->{'polar'}};
326 $self->{c_dirty} = 0;
327 return $self->{'cartesian'} = [$r * CORE::cos($t), $r * CORE::sin($t)];
334 # Recompute and return the polar form, given accurate cartesian form.
338 my ($x, $y) = @{$self->{'cartesian'}};
339 $self->{p_dirty} = 0;
340 return $self->{'polar'} = [0, 0] if $x == 0 && $y == 0;
341 return $self->{'polar'} = [CORE::sqrt($x*$x + $y*$y),
342 CORE::atan2($y, $x)];
351 my ($z1, $z2, $regular) = @_;
352 my ($re1, $im1) = @{$z1->cartesian};
353 $z2 = cplx($z2) unless ref $z2;
354 my ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
355 unless (defined $regular) {
356 $z1->set_cartesian([$re1 + $re2, $im1 + $im2]);
359 return (ref $z1)->make($re1 + $re2, $im1 + $im2);
368 my ($z1, $z2, $inverted) = @_;
369 my ($re1, $im1) = @{$z1->cartesian};
370 $z2 = cplx($z2) unless ref $z2;
371 my ($re2, $im2) = @{$z2->cartesian};
372 unless (defined $inverted) {
373 $z1->set_cartesian([$re1 - $re2, $im1 - $im2]);
377 (ref $z1)->make($re2 - $re1, $im2 - $im1) :
378 (ref $z1)->make($re1 - $re2, $im1 - $im2);
388 my ($z1, $z2, $regular) = @_;
389 if ($z1->{p_dirty} == 0 and ref $z2 and $z2->{p_dirty} == 0) {
390 # if both polar better use polar to avoid rounding errors
391 my ($r1, $t1) = @{$z1->polar};
392 my ($r2, $t2) = @{$z2->polar};
394 if ($t > pi()) { $t -= pit2 }
395 elsif ($t <= -pi()) { $t += pit2 }
396 unless (defined $regular) {
397 $z1->set_polar([$r1 * $r2, $t]);
400 return (ref $z1)->emake($r1 * $r2, $t);
402 my ($x1, $y1) = @{$z1->cartesian};
404 my ($x2, $y2) = @{$z2->cartesian};
405 return (ref $z1)->make($x1*$x2-$y1*$y2, $x1*$y2+$y1*$x2);
407 return (ref $z1)->make($x1*$z2, $y1*$z2);
415 # Die on division by zero.
418 my $mess = "$_[0]: Division by zero.\n";
421 $mess .= "(Because in the definition of $_[0], the divisor ";
422 $mess .= "$_[1] " unless ("$_[1]" eq '0');
428 $mess .= "Died at $up[1] line $up[2].\n";
439 my ($z1, $z2, $inverted) = @_;
440 if ($z1->{p_dirty} == 0 and ref $z2 and $z2->{p_dirty} == 0) {
441 # if both polar better use polar to avoid rounding errors
442 my ($r1, $t1) = @{$z1->polar};
443 my ($r2, $t2) = @{$z2->polar};
446 _divbyzero "$z2/0" if ($r1 == 0);
448 if ($t > pi()) { $t -= pit2 }
449 elsif ($t <= -pi()) { $t += pit2 }
450 return (ref $z1)->emake($r2 / $r1, $t);
452 _divbyzero "$z1/0" if ($r2 == 0);
454 if ($t > pi()) { $t -= pit2 }
455 elsif ($t <= -pi()) { $t += pit2 }
456 return (ref $z1)->emake($r1 / $r2, $t);
461 ($x2, $y2) = @{$z1->cartesian};
462 $d = $x2*$x2 + $y2*$y2;
463 _divbyzero "$z2/0" if $d == 0;
464 return (ref $z1)->make(($x2*$z2)/$d, -($y2*$z2)/$d);
466 my ($x1, $y1) = @{$z1->cartesian};
468 ($x2, $y2) = @{$z2->cartesian};
469 $d = $x2*$x2 + $y2*$y2;
470 _divbyzero "$z1/0" if $d == 0;
471 my $u = ($x1*$x2 + $y1*$y2)/$d;
472 my $v = ($y1*$x2 - $x1*$y2)/$d;
473 return (ref $z1)->make($u, $v);
475 _divbyzero "$z1/0" if $z2 == 0;
476 return (ref $z1)->make($x1/$z2, $y1/$z2);
485 # Computes z1**z2 = exp(z2 * log z1)).
488 my ($z1, $z2, $inverted) = @_;
490 return 1 if $z1 == 0 || $z2 == 1;
491 return 0 if $z2 == 0 && Re($z1) > 0;
493 return 1 if $z2 == 0 || $z1 == 1;
494 return 0 if $z1 == 0 && Re($z2) > 0;
496 my $w = $inverted ? &exp($z1 * &log($z2))
497 : &exp($z2 * &log($z1));
498 # If both arguments cartesian, return cartesian, else polar.
499 return $z1->{c_dirty} == 0 &&
500 (not ref $z2 or $z2->{c_dirty} == 0) ?
501 cplx(@{$w->cartesian}) : $w;
507 # Computes z1 <=> z2.
508 # Sorts on the real part first, then on the imaginary part. Thus 2-4i < 3+8i.
511 my ($z1, $z2, $inverted) = @_;
512 my ($re1, $im1) = ref $z1 ? @{$z1->cartesian} : ($z1, 0);
513 my ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
514 my $sgn = $inverted ? -1 : 1;
515 return $sgn * ($re1 <=> $re2) if $re1 != $re2;
516 return $sgn * ($im1 <=> $im2);
524 # (Required in addition to spaceship() because of NaNs.)
526 my ($z1, $z2, $inverted) = @_;
527 my ($re1, $im1) = ref $z1 ? @{$z1->cartesian} : ($z1, 0);
528 my ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
529 return $re1 == $re2 && $im1 == $im2 ? 1 : 0;
540 my ($r, $t) = @{$z->polar};
541 $t = ($t <= 0) ? $t + pi : $t - pi;
542 return (ref $z)->emake($r, $t);
544 my ($re, $im) = @{$z->cartesian};
545 return (ref $z)->make(-$re, -$im);
551 # Compute complex's conjugate.
556 my ($r, $t) = @{$z->polar};
557 return (ref $z)->emake($r, -$t);
559 my ($re, $im) = @{$z->cartesian};
560 return (ref $z)->make($re, -$im);
566 # Compute or set complex's norm (rho).
574 return CORE::abs($z);
578 $z->{'polar'} = [ $rho, ${$z->polar}[1] ];
583 return ${$z->polar}[0];
590 if ($$theta > pi()) { $$theta -= pit2 }
591 elsif ($$theta <= -pi()) { $$theta += pit2 }
597 # Compute or set complex's argument (theta).
600 my ($z, $theta) = @_;
601 return $z unless ref $z;
602 if (defined $theta) {
604 $z->{'polar'} = [ ${$z->polar}[0], $theta ];
608 $theta = ${$z->polar}[1];
619 # It is quite tempting to use wantarray here so that in list context
620 # sqrt() would return the two solutions. This, however, would
623 # print "sqrt(z) = ", sqrt($z), "\n";
625 # The two values would be printed side by side without no intervening
626 # whitespace, quite confusing.
627 # Therefore if you want the two solutions use the root().
631 my ($re, $im) = ref $z ? @{$z->cartesian} : ($z, 0);
632 return $re < 0 ? cplx(0, CORE::sqrt(-$re)) : CORE::sqrt($re)
634 my ($r, $t) = @{$z->polar};
635 return (ref $z)->emake(CORE::sqrt($r), $t/2);
641 # Compute cbrt(z) (cubic root).
643 # Why are we not returning three values? The same answer as for sqrt().
648 -CORE::exp(CORE::log(-$z)/3) :
649 ($z > 0 ? CORE::exp(CORE::log($z)/3): 0)
651 my ($r, $t) = @{$z->polar};
653 return (ref $z)->emake(CORE::exp(CORE::log($r)/3), $t/3);
662 my $mess = "Root $_[0] illegal, root rank must be positive integer.\n";
666 $mess .= "Died at $up[1] line $up[2].\n";
674 # Computes all nth root for z, returning an array whose size is n.
675 # `n' must be a positive integer.
677 # The roots are given by (for k = 0..n-1):
679 # z^(1/n) = r^(1/n) (cos ((t+2 k pi)/n) + i sin ((t+2 k pi)/n))
683 _rootbad($n) if ($n < 1 or int($n) != $n);
684 my ($r, $t) = ref $z ?
685 @{$z->polar} : (CORE::abs($z), $z >= 0 ? 0 : pi);
688 my $theta_inc = pit2 / $n;
689 my $rho = $r ** (1/$n);
691 my $cartesian = ref $z && $z->{c_dirty} == 0;
692 for ($k = 0, $theta = $t / $n; $k < $n; $k++, $theta += $theta_inc) {
693 my $w = cplxe($rho, $theta);
694 # Yes, $cartesian is loop invariant.
695 push @root, $cartesian ? cplx(@{$w->cartesian}) : $w;
703 # Return or set Re(z).
707 return $z unless ref $z;
709 $z->{'cartesian'} = [ $Re, ${$z->cartesian}[1] ];
713 return ${$z->cartesian}[0];
720 # Return or set Im(z).
724 return 0 unless ref $z;
726 $z->{'cartesian'} = [ ${$z->cartesian}[0], $Im ];
730 return ${$z->cartesian}[1];
737 # Return or set rho(w).
740 Math::Complex::abs(@_);
746 # Return or set theta(w).
749 Math::Complex::arg(@_);
759 my ($x, $y) = @{$z->cartesian};
760 return (ref $z)->emake(CORE::exp($x), $y);
766 # Die on logarithm of zero.
769 my $mess = "$_[0]: Logarithm of zero.\n";
772 $mess .= "(Because in the definition of $_[0], the argument ";
773 $mess .= "$_[1] " unless ($_[1] eq '0');
779 $mess .= "Died at $up[1] line $up[2].\n";
792 _logofzero("log") if $z == 0;
793 return $z > 0 ? CORE::log($z) : cplx(CORE::log(-$z), pi);
795 my ($r, $t) = @{$z->polar};
796 _logofzero("log") if $r == 0;
797 if ($t > pi()) { $t -= pit2 }
798 elsif ($t <= -pi()) { $t += pit2 }
799 return (ref $z)->make(CORE::log($r), $t);
807 sub ln { Math::Complex::log(@_) }
816 return Math::Complex::log($_[0]) * uplog10;
822 # Compute logn(z,n) = log(z) / log(n)
826 $z = cplx($z, 0) unless ref $z;
827 my $logn = $LOGN{$n};
828 $logn = $LOGN{$n} = CORE::log($n) unless defined $logn; # Cache log(n)
829 return &log($z) / $logn;
835 # Compute cos(z) = (exp(iz) + exp(-iz))/2.
839 return CORE::cos($z) unless ref $z;
840 my ($x, $y) = @{$z->cartesian};
841 my $ey = CORE::exp($y);
842 my $sx = CORE::sin($x);
843 my $cx = CORE::cos($x);
844 my $ey_1 = $ey ? 1 / $ey : $Inf;
845 return (ref $z)->make($cx * ($ey + $ey_1)/2,
846 $sx * ($ey_1 - $ey)/2);
852 # Compute sin(z) = (exp(iz) - exp(-iz))/2.
856 return CORE::sin($z) unless ref $z;
857 my ($x, $y) = @{$z->cartesian};
858 my $ey = CORE::exp($y);
859 my $sx = CORE::sin($x);
860 my $cx = CORE::cos($x);
861 my $ey_1 = $ey ? 1 / $ey : $Inf;
862 return (ref $z)->make($sx * ($ey + $ey_1)/2,
863 $cx * ($ey - $ey_1)/2);
869 # Compute tan(z) = sin(z) / cos(z).
874 _divbyzero "tan($z)", "cos($z)" if $cz == 0;
875 return &sin($z) / $cz;
881 # Computes the secant sec(z) = 1 / cos(z).
886 _divbyzero "sec($z)", "cos($z)" if ($cz == 0);
893 # Computes the cosecant csc(z) = 1 / sin(z).
898 _divbyzero "csc($z)", "sin($z)" if ($sz == 0);
907 sub cosec { Math::Complex::csc(@_) }
912 # Computes cot(z) = cos(z) / sin(z).
917 _divbyzero "cot($z)", "sin($z)" if ($sz == 0);
918 return &cos($z) / $sz;
926 sub cotan { Math::Complex::cot(@_) }
931 # Computes the arc cosine acos(z) = -i log(z + sqrt(z*z-1)).
935 return CORE::atan2(CORE::sqrt(1-$z*$z), $z)
936 if (! ref $z) && CORE::abs($z) <= 1;
937 $z = cplx($z, 0) unless ref $z;
938 my ($x, $y) = @{$z->cartesian};
939 return 0 if $x == 1 && $y == 0;
940 my $t1 = CORE::sqrt(($x+1)*($x+1) + $y*$y);
941 my $t2 = CORE::sqrt(($x-1)*($x-1) + $y*$y);
942 my $alpha = ($t1 + $t2)/2;
943 my $beta = ($t1 - $t2)/2;
944 $alpha = 1 if $alpha < 1;
945 if ($beta > 1) { $beta = 1 }
946 elsif ($beta < -1) { $beta = -1 }
947 my $u = CORE::atan2(CORE::sqrt(1-$beta*$beta), $beta);
948 my $v = CORE::log($alpha + CORE::sqrt($alpha*$alpha-1));
949 $v = -$v if $y > 0 || ($y == 0 && $x < -1);
950 return (ref $z)->make($u, $v);
956 # Computes the arc sine asin(z) = -i log(iz + sqrt(1-z*z)).
960 return CORE::atan2($z, CORE::sqrt(1-$z*$z))
961 if (! ref $z) && CORE::abs($z) <= 1;
962 $z = cplx($z, 0) unless ref $z;
963 my ($x, $y) = @{$z->cartesian};
964 return 0 if $x == 0 && $y == 0;
965 my $t1 = CORE::sqrt(($x+1)*($x+1) + $y*$y);
966 my $t2 = CORE::sqrt(($x-1)*($x-1) + $y*$y);
967 my $alpha = ($t1 + $t2)/2;
968 my $beta = ($t1 - $t2)/2;
969 $alpha = 1 if $alpha < 1;
970 if ($beta > 1) { $beta = 1 }
971 elsif ($beta < -1) { $beta = -1 }
972 my $u = CORE::atan2($beta, CORE::sqrt(1-$beta*$beta));
973 my $v = -CORE::log($alpha + CORE::sqrt($alpha*$alpha-1));
974 $v = -$v if $y > 0 || ($y == 0 && $x < -1);
975 return (ref $z)->make($u, $v);
981 # Computes the arc tangent atan(z) = i/2 log((i+z) / (i-z)).
985 return CORE::atan2($z, 1) unless ref $z;
986 my ($x, $y) = ref $z ? @{$z->cartesian} : ($z, 0);
987 return 0 if $x == 0 && $y == 0;
988 _divbyzero "atan(i)" if ( $z == i);
989 _logofzero "atan(-i)" if (-$z == i); # -i is a bad file test...
990 my $log = &log((i + $z) / (i - $z));
997 # Computes the arc secant asec(z) = acos(1 / z).
1001 _divbyzero "asec($z)", $z if ($z == 0);
1002 return acos(1 / $z);
1008 # Computes the arc cosecant acsc(z) = asin(1 / z).
1012 _divbyzero "acsc($z)", $z if ($z == 0);
1013 return asin(1 / $z);
1021 sub acosec { Math::Complex::acsc(@_) }
1026 # Computes the arc cotangent acot(z) = atan(1 / z)
1030 _divbyzero "acot(0)" if $z == 0;
1031 return ($z >= 0) ? CORE::atan2(1, $z) : CORE::atan2(-1, -$z)
1033 _divbyzero "acot(i)" if ($z - i == 0);
1034 _logofzero "acot(-i)" if ($z + i == 0);
1035 return atan(1 / $z);
1043 sub acotan { Math::Complex::acot(@_) }
1048 # Computes the hyperbolic cosine cosh(z) = (exp(z) + exp(-z))/2.
1054 $ex = CORE::exp($z);
1055 return $ex ? ($ex + 1/$ex)/2 : $Inf;
1057 my ($x, $y) = @{$z->cartesian};
1058 $ex = CORE::exp($x);
1059 my $ex_1 = $ex ? 1 / $ex : $Inf;
1060 return (ref $z)->make(CORE::cos($y) * ($ex + $ex_1)/2,
1061 CORE::sin($y) * ($ex - $ex_1)/2);
1067 # Computes the hyperbolic sine sinh(z) = (exp(z) - exp(-z))/2.
1073 return 0 if $z == 0;
1074 $ex = CORE::exp($z);
1075 return $ex ? ($ex - 1/$ex)/2 : "-$Inf";
1077 my ($x, $y) = @{$z->cartesian};
1078 my $cy = CORE::cos($y);
1079 my $sy = CORE::sin($y);
1080 $ex = CORE::exp($x);
1081 my $ex_1 = $ex ? 1 / $ex : $Inf;
1082 return (ref $z)->make(CORE::cos($y) * ($ex - $ex_1)/2,
1083 CORE::sin($y) * ($ex + $ex_1)/2);
1089 # Computes the hyperbolic tangent tanh(z) = sinh(z) / cosh(z).
1094 _divbyzero "tanh($z)", "cosh($z)" if ($cz == 0);
1095 return sinh($z) / $cz;
1101 # Computes the hyperbolic secant sech(z) = 1 / cosh(z).
1106 _divbyzero "sech($z)", "cosh($z)" if ($cz == 0);
1113 # Computes the hyperbolic cosecant csch(z) = 1 / sinh(z).
1118 _divbyzero "csch($z)", "sinh($z)" if ($sz == 0);
1127 sub cosech { Math::Complex::csch(@_) }
1132 # Computes the hyperbolic cotangent coth(z) = cosh(z) / sinh(z).
1137 _divbyzero "coth($z)", "sinh($z)" if $sz == 0;
1138 return cosh($z) / $sz;
1146 sub cotanh { Math::Complex::coth(@_) }
1151 # Computes the arc hyperbolic cosine acosh(z) = log(z + sqrt(z*z-1)).
1158 my ($re, $im) = @{$z->cartesian};
1160 return CORE::log($re + CORE::sqrt($re*$re - 1))
1162 return cplx(0, CORE::atan2(CORE::sqrt(1 - $re*$re), $re))
1163 if CORE::abs($re) < 1;
1165 my $t = &sqrt($z * $z - 1) + $z;
1166 # Try Taylor if looking bad (this usually means that
1167 # $z was large negative, therefore the sqrt is really
1168 # close to abs(z), summing that with z...)
1169 $t = 1/(2 * $z) - 1/(8 * $z**3) + 1/(16 * $z**5) - 5/(128 * $z**7)
1172 $u->Im(-$u->Im) if $re < 0 && $im == 0;
1173 return $re < 0 ? -$u : $u;
1179 # Computes the arc hyperbolic sine asinh(z) = log(z + sqrt(z*z+1))
1184 my $t = $z + CORE::sqrt($z*$z + 1);
1185 return CORE::log($t) if $t;
1187 my $t = &sqrt($z * $z + 1) + $z;
1188 # Try Taylor if looking bad (this usually means that
1189 # $z was large negative, therefore the sqrt is really
1190 # close to abs(z), summing that with z...)
1191 $t = 1/(2 * $z) - 1/(8 * $z**3) + 1/(16 * $z**5) - 5/(128 * $z**7)
1199 # Computes the arc hyperbolic tangent atanh(z) = 1/2 log((1+z) / (1-z)).
1204 return CORE::log((1 + $z)/(1 - $z))/2 if CORE::abs($z) < 1;
1207 _divbyzero 'atanh(1)', "1 - $z" if (1 - $z == 0);
1208 _logofzero 'atanh(-1)' if (1 + $z == 0);
1209 return 0.5 * &log((1 + $z) / (1 - $z));
1215 # Computes the hyperbolic arc secant asech(z) = acosh(1 / z).
1219 _divbyzero 'asech(0)', "$z" if ($z == 0);
1220 return acosh(1 / $z);
1226 # Computes the hyperbolic arc cosecant acsch(z) = asinh(1 / z).
1230 _divbyzero 'acsch(0)', $z if ($z == 0);
1231 return asinh(1 / $z);
1237 # Alias for acosh().
1239 sub acosech { Math::Complex::acsch(@_) }
1244 # Computes the arc hyperbolic cotangent acoth(z) = 1/2 log((1+z) / (z-1)).
1248 _divbyzero 'acoth(0)' if ($z == 0);
1250 return CORE::log(($z + 1)/($z - 1))/2 if CORE::abs($z) > 1;
1253 _divbyzero 'acoth(1)', "$z - 1" if ($z - 1 == 0);
1254 _logofzero 'acoth(-1)', "1 + $z" if (1 + $z == 0);
1255 return &log((1 + $z) / ($z - 1)) / 2;
1263 sub acotanh { Math::Complex::acoth(@_) }
1268 # Compute atan(z1/z2).
1271 my ($z1, $z2, $inverted) = @_;
1272 my ($re1, $im1, $re2, $im2);
1274 ($re1, $im1) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
1275 ($re2, $im2) = @{$z1->cartesian};
1277 ($re1, $im1) = @{$z1->cartesian};
1278 ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
1281 return CORE::atan2($re1, $re2) if $im1 == 0;
1282 return ($im1<=>0) * pip2 if $re2 == 0;
1284 my $w = atan($z1/$z2);
1285 my ($u, $v) = ref $w ? @{$w->cartesian} : ($w, 0);
1286 $u += pi if $re2 < 0;
1287 $u -= pit2 if $u > pi;
1288 return cplx($u, $v);
1295 # Set (get if no argument) the display format for all complex numbers that
1296 # don't happen to have overridden it via ->display_format
1298 # When called as an object method, this actually sets the display format for
1299 # the current object.
1301 # Valid object formats are 'c' and 'p' for cartesian and polar. The first
1302 # letter is used actually, so the type can be fully spelled out for clarity.
1304 sub display_format {
1306 my %display_format = %DISPLAY_FORMAT;
1308 if (ref $self) { # Called as an object method
1309 if (exists $self->{display_format}) {
1310 my %obj = %{$self->{display_format}};
1311 @display_format{keys %obj} = values %obj;
1315 $display_format{style} = shift;
1318 @display_format{keys %new} = values %new;
1321 if (ref $self) { # Called as an object method
1322 $self->{display_format} = { %display_format };
1325 %{$self->{display_format}} :
1326 $self->{display_format}->{style};
1329 # Called as a class method
1330 %DISPLAY_FORMAT = %display_format;
1334 $DISPLAY_FORMAT{style};
1340 # Show nicely formatted complex number under its cartesian or polar form,
1341 # depending on the current display format:
1343 # . If a specific display format has been recorded for this object, use it.
1344 # . Otherwise, use the generic current default for all complex numbers,
1345 # which is a package global variable.
1350 my $style = $z->display_format;
1352 $style = $DISPLAY_FORMAT{style} unless defined $style;
1354 return $z->stringify_polar if $style =~ /^p/i;
1355 return $z->stringify_cartesian;
1359 # ->stringify_cartesian
1361 # Stringify as a cartesian representation 'a+bi'.
1363 sub stringify_cartesian {
1365 my ($x, $y) = @{$z->cartesian};
1368 my %format = $z->display_format;
1369 my $format = $format{format};
1372 if ($x =~ /^NaN[QS]?$/i) {
1375 if ($x =~ /^-?$Inf$/oi) {
1378 $re = defined $format ? sprintf($format, $x) : $x;
1386 if ($y =~ /^(NaN[QS]?)$/i) {
1389 if ($y =~ /^-?$Inf$/oi) {
1394 sprintf($format, $y) :
1395 ($y == 1 ? "" : ($y == -1 ? "-" : $y));
1408 } elsif ($y > 0 || $im =~ /^NaN[QS]?i$/i) {
1409 $str .= "+" if defined $re;
1412 } elsif (!defined $re) {
1423 # Stringify as a polar representation '[r,t]'.
1425 sub stringify_polar {
1427 my ($r, $t) = @{$z->polar};
1430 my %format = $z->display_format;
1431 my $format = $format{format};
1433 if ($t =~ /^NaN[QS]?$/i || $t =~ /^-?$Inf$/oi) {
1435 } elsif ($t == pi) {
1437 } elsif ($r == 0 || $t == 0) {
1438 $theta = defined $format ? sprintf($format, $t) : $t;
1441 return "[$r,$theta]" if defined $theta;
1444 # Try to identify pi/n and friends.
1447 $t -= int(CORE::abs($t) / pit2) * pit2;
1449 if ($format{polar_pretty_print} && $t) {
1453 if ($b =~ /^-?\d+$/) {
1454 $b = $b < 0 ? "-" : "" if CORE::abs($b) == 1;
1455 $theta = "${b}pi/$a";
1461 if (defined $format) {
1462 $r = sprintf($format, $r);
1463 $theta = sprintf($format, $theta) unless defined $theta;
1465 $theta = $t unless defined $theta;
1468 return "[$r,$theta]";
1478 Math::Complex - complex numbers and associated mathematical functions
1484 $z = Math::Complex->make(5, 6);
1486 $j = cplxe(1, 2*pi/3);
1490 This package lets you create and manipulate complex numbers. By default,
1491 I<Perl> limits itself to real numbers, but an extra C<use> statement brings
1492 full complex support, along with a full set of mathematical functions
1493 typically associated with and/or extended to complex numbers.
1495 If you wonder what complex numbers are, they were invented to be able to solve
1496 the following equation:
1500 and by definition, the solution is noted I<i> (engineers use I<j> instead since
1501 I<i> usually denotes an intensity, but the name does not matter). The number
1502 I<i> is a pure I<imaginary> number.
1504 The arithmetics with pure imaginary numbers works just like you would expect
1505 it with real numbers... you just have to remember that
1511 5i + 7i = i * (5 + 7) = 12i
1512 4i - 3i = i * (4 - 3) = i
1517 Complex numbers are numbers that have both a real part and an imaginary
1518 part, and are usually noted:
1522 where C<a> is the I<real> part and C<b> is the I<imaginary> part. The
1523 arithmetic with complex numbers is straightforward. You have to
1524 keep track of the real and the imaginary parts, but otherwise the
1525 rules used for real numbers just apply:
1527 (4 + 3i) + (5 - 2i) = (4 + 5) + i(3 - 2) = 9 + i
1528 (2 + i) * (4 - i) = 2*4 + 4i -2i -i*i = 8 + 2i + 1 = 9 + 2i
1530 A graphical representation of complex numbers is possible in a plane
1531 (also called the I<complex plane>, but it's really a 2D plane).
1536 is the point whose coordinates are (a, b). Actually, it would
1537 be the vector originating from (0, 0) to (a, b). It follows that the addition
1538 of two complex numbers is a vectorial addition.
1540 Since there is a bijection between a point in the 2D plane and a complex
1541 number (i.e. the mapping is unique and reciprocal), a complex number
1542 can also be uniquely identified with polar coordinates:
1546 where C<rho> is the distance to the origin, and C<theta> the angle between
1547 the vector and the I<x> axis. There is a notation for this using the
1548 exponential form, which is:
1550 rho * exp(i * theta)
1552 where I<i> is the famous imaginary number introduced above. Conversion
1553 between this form and the cartesian form C<a + bi> is immediate:
1555 a = rho * cos(theta)
1556 b = rho * sin(theta)
1558 which is also expressed by this formula:
1560 z = rho * exp(i * theta) = rho * (cos theta + i * sin theta)
1562 In other words, it's the projection of the vector onto the I<x> and I<y>
1563 axes. Mathematicians call I<rho> the I<norm> or I<modulus> and I<theta>
1564 the I<argument> of the complex number. The I<norm> of C<z> will be
1567 The polar notation (also known as the trigonometric
1568 representation) is much more handy for performing multiplications and
1569 divisions of complex numbers, whilst the cartesian notation is better
1570 suited for additions and subtractions. Real numbers are on the I<x>
1571 axis, and therefore I<theta> is zero or I<pi>.
1573 All the common operations that can be performed on a real number have
1574 been defined to work on complex numbers as well, and are merely
1575 I<extensions> of the operations defined on real numbers. This means
1576 they keep their natural meaning when there is no imaginary part, provided
1577 the number is within their definition set.
1579 For instance, the C<sqrt> routine which computes the square root of
1580 its argument is only defined for non-negative real numbers and yields a
1581 non-negative real number (it is an application from B<R+> to B<R+>).
1582 If we allow it to return a complex number, then it can be extended to
1583 negative real numbers to become an application from B<R> to B<C> (the
1584 set of complex numbers):
1586 sqrt(x) = x >= 0 ? sqrt(x) : sqrt(-x)*i
1588 It can also be extended to be an application from B<C> to B<C>,
1589 whilst its restriction to B<R> behaves as defined above by using
1590 the following definition:
1592 sqrt(z = [r,t]) = sqrt(r) * exp(i * t/2)
1594 Indeed, a negative real number can be noted C<[x,pi]> (the modulus
1595 I<x> is always non-negative, so C<[x,pi]> is really C<-x>, a negative
1596 number) and the above definition states that
1598 sqrt([x,pi]) = sqrt(x) * exp(i*pi/2) = [sqrt(x),pi/2] = sqrt(x)*i
1600 which is exactly what we had defined for negative real numbers above.
1601 The C<sqrt> returns only one of the solutions: if you want the both,
1602 use the C<root> function.
1604 All the common mathematical functions defined on real numbers that
1605 are extended to complex numbers share that same property of working
1606 I<as usual> when the imaginary part is zero (otherwise, it would not
1607 be called an extension, would it?).
1609 A I<new> operation possible on a complex number that is
1610 the identity for real numbers is called the I<conjugate>, and is noted
1611 with a horizontal bar above the number, or C<~z> here.
1618 z * ~z = (a + bi) * (a - bi) = a*a + b*b
1620 We saw that the norm of C<z> was noted C<abs(z)> and was defined as the
1621 distance to the origin, also known as:
1623 rho = abs(z) = sqrt(a*a + b*b)
1627 z * ~z = abs(z) ** 2
1629 If z is a pure real number (i.e. C<b == 0>), then the above yields:
1633 which is true (C<abs> has the regular meaning for real number, i.e. stands
1634 for the absolute value). This example explains why the norm of C<z> is
1635 noted C<abs(z)>: it extends the C<abs> function to complex numbers, yet
1636 is the regular C<abs> we know when the complex number actually has no
1637 imaginary part... This justifies I<a posteriori> our use of the C<abs>
1638 notation for the norm.
1642 Given the following notations:
1644 z1 = a + bi = r1 * exp(i * t1)
1645 z2 = c + di = r2 * exp(i * t2)
1646 z = <any complex or real number>
1648 the following (overloaded) operations are supported on complex numbers:
1650 z1 + z2 = (a + c) + i(b + d)
1651 z1 - z2 = (a - c) + i(b - d)
1652 z1 * z2 = (r1 * r2) * exp(i * (t1 + t2))
1653 z1 / z2 = (r1 / r2) * exp(i * (t1 - t2))
1654 z1 ** z2 = exp(z2 * log z1)
1656 abs(z) = r1 = sqrt(a*a + b*b)
1657 sqrt(z) = sqrt(r1) * exp(i * t/2)
1658 exp(z) = exp(a) * exp(i * b)
1659 log(z) = log(r1) + i*t
1660 sin(z) = 1/2i (exp(i * z1) - exp(-i * z))
1661 cos(z) = 1/2 (exp(i * z1) + exp(-i * z))
1662 atan2(z1, z2) = atan(z1/z2)
1664 The following extra operations are supported on both real and complex
1672 cbrt(z) = z ** (1/3)
1673 log10(z) = log(z) / log(10)
1674 logn(z, n) = log(z) / log(n)
1676 tan(z) = sin(z) / cos(z)
1682 asin(z) = -i * log(i*z + sqrt(1-z*z))
1683 acos(z) = -i * log(z + i*sqrt(1-z*z))
1684 atan(z) = i/2 * log((i+z) / (i-z))
1686 acsc(z) = asin(1 / z)
1687 asec(z) = acos(1 / z)
1688 acot(z) = atan(1 / z) = -i/2 * log((i+z) / (z-i))
1690 sinh(z) = 1/2 (exp(z) - exp(-z))
1691 cosh(z) = 1/2 (exp(z) + exp(-z))
1692 tanh(z) = sinh(z) / cosh(z) = (exp(z) - exp(-z)) / (exp(z) + exp(-z))
1694 csch(z) = 1 / sinh(z)
1695 sech(z) = 1 / cosh(z)
1696 coth(z) = 1 / tanh(z)
1698 asinh(z) = log(z + sqrt(z*z+1))
1699 acosh(z) = log(z + sqrt(z*z-1))
1700 atanh(z) = 1/2 * log((1+z) / (1-z))
1702 acsch(z) = asinh(1 / z)
1703 asech(z) = acosh(1 / z)
1704 acoth(z) = atanh(1 / z) = 1/2 * log((1+z) / (z-1))
1706 I<arg>, I<abs>, I<log>, I<csc>, I<cot>, I<acsc>, I<acot>, I<csch>,
1707 I<coth>, I<acosech>, I<acotanh>, have aliases I<rho>, I<theta>, I<ln>,
1708 I<cosec>, I<cotan>, I<acosec>, I<acotan>, I<cosech>, I<cotanh>,
1709 I<acosech>, I<acotanh>, respectively. C<Re>, C<Im>, C<arg>, C<abs>,
1710 C<rho>, and C<theta> can be used also as mutators. The C<cbrt>
1711 returns only one of the solutions: if you want all three, use the
1714 The I<root> function is available to compute all the I<n>
1715 roots of some complex, where I<n> is a strictly positive integer.
1716 There are exactly I<n> such roots, returned as a list. Getting the
1717 number mathematicians call C<j> such that:
1721 is a simple matter of writing:
1723 $j = ((root(1, 3))[1];
1725 The I<k>th root for C<z = [r,t]> is given by:
1727 (root(z, n))[k] = r**(1/n) * exp(i * (t + 2*k*pi)/n)
1729 The I<spaceship> comparison operator, E<lt>=E<gt>, is also defined. In
1730 order to ensure its restriction to real numbers is conform to what you
1731 would expect, the comparison is run on the real part of the complex
1732 number first, and imaginary parts are compared only when the real
1737 To create a complex number, use either:
1739 $z = Math::Complex->make(3, 4);
1742 if you know the cartesian form of the number, or
1746 if you like. To create a number using the polar form, use either:
1748 $z = Math::Complex->emake(5, pi/3);
1749 $x = cplxe(5, pi/3);
1751 instead. The first argument is the modulus, the second is the angle
1752 (in radians, the full circle is 2*pi). (Mnemonic: C<e> is used as a
1753 notation for complex numbers in the polar form).
1755 It is possible to write:
1757 $x = cplxe(-3, pi/4);
1759 but that will be silently converted into C<[3,-3pi/4]>, since the
1760 modulus must be non-negative (it represents the distance to the origin
1761 in the complex plane).
1763 It is also possible to have a complex number as either argument of the
1764 C<make>, C<emake>, C<cplx>, and C<cplxe>: the appropriate component of
1765 the argument will be used.
1770 The C<new>, C<make>, C<emake>, C<cplx>, and C<cplxe> will also
1771 understand a single (string) argument of the forms
1778 in which case the appropriate cartesian and exponential components
1779 will be parsed from the string and used to create new complex numbers.
1780 The imaginary component and the theta, respectively, will default to zero.
1782 =head1 STRINGIFICATION
1784 When printed, a complex number is usually shown under its cartesian
1785 style I<a+bi>, but there are legitimate cases where the polar style
1786 I<[r,t]> is more appropriate.
1788 By calling the class method C<Math::Complex::display_format> and
1789 supplying either C<"polar"> or C<"cartesian"> as an argument, you
1790 override the default display style, which is C<"cartesian">. Not
1791 supplying any argument returns the current settings.
1793 This default can be overridden on a per-number basis by calling the
1794 C<display_format> method instead. As before, not supplying any argument
1795 returns the current display style for this number. Otherwise whatever you
1796 specify will be the new display style for I<this> particular number.
1802 Math::Complex::display_format('polar');
1803 $j = (root(1, 3))[1];
1804 print "j = $j\n"; # Prints "j = [1,2pi/3]"
1805 $j->display_format('cartesian');
1806 print "j = $j\n"; # Prints "j = -0.5+0.866025403784439i"
1808 The polar style attempts to emphasize arguments like I<k*pi/n>
1809 (where I<n> is a positive integer and I<k> an integer within [-9, +9]),
1810 this is called I<polar pretty-printing>.
1812 =head2 CHANGED IN PERL 5.6
1814 The C<display_format> class method and the corresponding
1815 C<display_format> object method can now be called using
1816 a parameter hash instead of just a one parameter.
1818 The old display format style, which can have values C<"cartesian"> or
1819 C<"polar">, can be changed using the C<"style"> parameter.
1821 $j->display_format(style => "polar");
1823 The one parameter calling convention also still works.
1825 $j->display_format("polar");
1827 There are two new display parameters.
1829 The first one is C<"format">, which is a sprintf()-style format string
1830 to be used for both numeric parts of the complex number(s). The is
1831 somewhat system-dependent but most often it corresponds to C<"%.15g">.
1832 You can revert to the default by setting the C<format> to C<undef>.
1834 # the $j from the above example
1836 $j->display_format('format' => '%.5f');
1837 print "j = $j\n"; # Prints "j = -0.50000+0.86603i"
1838 $j->display_format('format' => undef);
1839 print "j = $j\n"; # Prints "j = -0.5+0.86603i"
1841 Notice that this affects also the return values of the
1842 C<display_format> methods: in list context the whole parameter hash
1843 will be returned, as opposed to only the style parameter value.
1844 This is a potential incompatibility with earlier versions if you
1845 have been calling the C<display_format> method in list context.
1847 The second new display parameter is C<"polar_pretty_print">, which can
1848 be set to true or false, the default being true. See the previous
1849 section for what this means.
1853 Thanks to overloading, the handling of arithmetics with complex numbers
1854 is simple and almost transparent.
1856 Here are some examples:
1860 $j = cplxe(1, 2*pi/3); # $j ** 3 == 1
1861 print "j = $j, j**3 = ", $j ** 3, "\n";
1862 print "1 + j + j**2 = ", 1 + $j + $j**2, "\n";
1864 $z = -16 + 0*i; # Force it to be a complex
1865 print "sqrt($z) = ", sqrt($z), "\n";
1867 $k = exp(i * 2*pi/3);
1868 print "$j - $k = ", $j - $k, "\n";
1870 $z->Re(3); # Re, Im, arg, abs,
1871 $j->arg(2); # (the last two aka rho, theta)
1872 # can be used also as mutators.
1874 =head1 ERRORS DUE TO DIVISION BY ZERO OR LOGARITHM OF ZERO
1876 The division (/) and the following functions
1882 atanh asech acsch acoth
1884 cannot be computed for all arguments because that would mean dividing
1885 by zero or taking logarithm of zero. These situations cause fatal
1886 runtime errors looking like this
1888 cot(0): Division by zero.
1889 (Because in the definition of cot(0), the divisor sin(0) is 0)
1894 atanh(-1): Logarithm of zero.
1897 For the C<csc>, C<cot>, C<asec>, C<acsc>, C<acot>, C<csch>, C<coth>,
1898 C<asech>, C<acsch>, the argument cannot be C<0> (zero). For the
1899 logarithmic functions and the C<atanh>, C<acoth>, the argument cannot
1900 be C<1> (one). For the C<atanh>, C<acoth>, the argument cannot be
1901 C<-1> (minus one). For the C<atan>, C<acot>, the argument cannot be
1902 C<i> (the imaginary unit). For the C<atan>, C<acoth>, the argument
1903 cannot be C<-i> (the negative imaginary unit). For the C<tan>,
1904 C<sec>, C<tanh>, the argument cannot be I<pi/2 + k * pi>, where I<k>
1907 Note that because we are operating on approximations of real numbers,
1908 these errors can happen when merely `too close' to the singularities
1911 =head1 ERRORS DUE TO INDIGESTIBLE ARGUMENTS
1913 The C<make> and C<emake> accept both real and complex arguments.
1914 When they cannot recognize the arguments they will die with error
1915 messages like the following
1917 Math::Complex::make: Cannot take real part of ...
1918 Math::Complex::make: Cannot take real part of ...
1919 Math::Complex::emake: Cannot take rho of ...
1920 Math::Complex::emake: Cannot take theta of ...
1924 Saying C<use Math::Complex;> exports many mathematical routines in the
1925 caller environment and even overrides some (C<sqrt>, C<log>).
1926 This is construed as a feature by the Authors, actually... ;-)
1928 All routines expect to be given real or complex numbers. Don't attempt to
1929 use BigFloat, since Perl has currently no rule to disambiguate a '+'
1930 operation (for instance) between two overloaded entities.
1932 In Cray UNICOS there is some strange numerical instability that results
1933 in root(), cos(), sin(), cosh(), sinh(), losing accuracy fast. Beware.
1934 The bug may be in UNICOS math libs, in UNICOS C compiler, in Math::Complex.
1935 Whatever it is, it does not manifest itself anywhere else where Perl runs.
1939 Daniel S. Lewart <F<d-lewart@uiuc.edu>>
1941 Original authors Raphael Manfredi <F<Raphael_Manfredi@pobox.com>> and
1942 Jarkko Hietaniemi <F<jhi@iki.fi>>