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 our($VERSION, @ISA, @EXPORT, %EXPORT_TAGS, $Inf);
16 $Inf = CORE::exp(CORE::exp(30)); # We do want an arithmetic overflow.
17 $! = $e; # Clear ERANGE.
18 undef $Inf unless $Inf =~ /^inf(?:inity)?$/i; # Inf INF inf Infinity
19 $Inf = "Inf" if !defined $Inf || !($Inf > 0); # Desperation.
34 csc cosec sec cot cotan
36 acsc acosec asec acot acotan
38 csch cosech sech coth cotanh
40 acsch acosech asech acoth acotanh
79 my %DISPLAY_FORMAT = ('style' => 'cartesian',
80 'polar_pretty_print' => 1);
81 my $eps = 1e-14; # Epsilon
84 # Object attributes (internal):
85 # cartesian [real, imaginary] -- cartesian form
86 # polar [rho, theta] -- polar form
87 # c_dirty cartesian form not up-to-date
88 # p_dirty polar form not up-to-date
89 # display display format (package's global when not set)
92 # Die on bad *make() arguments.
95 die "@{[(caller(1))[3]]}: Cannot take $_[0] of $_[1].\n";
101 # Create a new complex number (cartesian form)
104 my $self = bless {}, shift;
108 if ( $rre eq ref $self ) {
111 _cannot_make("real part", $rre);
116 if ( $rim eq ref $self ) {
119 _cannot_make("imaginary part", $rim);
122 $self->{'cartesian'} = [ $re, $im ];
123 $self->{c_dirty} = 0;
124 $self->{p_dirty} = 1;
125 $self->display_format('cartesian');
132 # Create a new complex number (exponential form)
135 my $self = bless {}, shift;
136 my ($rho, $theta) = @_;
139 if ( $rrh eq ref $self ) {
142 _cannot_make("rho", $rrh);
145 my $rth = ref $theta;
147 if ( $rth eq ref $self ) {
148 $theta = theta($theta);
150 _cannot_make("theta", $rth);
155 $theta = ($theta <= 0) ? $theta + pi() : $theta - pi();
157 $self->{'polar'} = [$rho, $theta];
158 $self->{p_dirty} = 0;
159 $self->{c_dirty} = 1;
160 $self->display_format('polar');
164 sub new { &make } # For backward compatibility only.
169 # Creates a complex number from a (re, im) tuple.
170 # This avoids the burden of writing Math::Complex->make(re, im).
174 return __PACKAGE__->make($re, defined $im ? $im : 0);
180 # Creates a complex number from a (rho, theta) tuple.
181 # This avoids the burden of writing Math::Complex->emake(rho, theta).
184 my ($rho, $theta) = @_;
185 return __PACKAGE__->emake($rho, defined $theta ? $theta : 0);
191 # The number defined as pi = 180 degrees
193 sub pi () { 4 * CORE::atan2(1, 1) }
200 sub pit2 () { 2 * pi }
207 sub pip2 () { pi / 2 }
212 # One degree in radians, used in stringify_polar.
215 sub deg1 () { pi / 180 }
222 sub uplog10 () { 1 / CORE::log(10) }
227 # The number defined as i*i = -1;
232 $i->{'cartesian'} = [0, 1];
233 $i->{'polar'} = [1, pip2];
247 # Attribute access/set routines
250 sub cartesian {$_[0]->{c_dirty} ?
251 $_[0]->update_cartesian : $_[0]->{'cartesian'}}
252 sub polar {$_[0]->{p_dirty} ?
253 $_[0]->update_polar : $_[0]->{'polar'}}
255 sub set_cartesian { $_[0]->{p_dirty}++; $_[0]->{'cartesian'} = $_[1] }
256 sub set_polar { $_[0]->{c_dirty}++; $_[0]->{'polar'} = $_[1] }
261 # Recompute and return the cartesian form, given accurate polar form.
263 sub update_cartesian {
265 my ($r, $t) = @{$self->{'polar'}};
266 $self->{c_dirty} = 0;
267 return $self->{'cartesian'} = [$r * CORE::cos($t), $r * CORE::sin($t)];
274 # Recompute and return the polar form, given accurate cartesian form.
278 my ($x, $y) = @{$self->{'cartesian'}};
279 $self->{p_dirty} = 0;
280 return $self->{'polar'} = [0, 0] if $x == 0 && $y == 0;
281 return $self->{'polar'} = [CORE::sqrt($x*$x + $y*$y),
282 CORE::atan2($y, $x)];
291 my ($z1, $z2, $regular) = @_;
292 my ($re1, $im1) = @{$z1->cartesian};
293 $z2 = cplx($z2) unless ref $z2;
294 my ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
295 unless (defined $regular) {
296 $z1->set_cartesian([$re1 + $re2, $im1 + $im2]);
299 return (ref $z1)->make($re1 + $re2, $im1 + $im2);
308 my ($z1, $z2, $inverted) = @_;
309 my ($re1, $im1) = @{$z1->cartesian};
310 $z2 = cplx($z2) unless ref $z2;
311 my ($re2, $im2) = @{$z2->cartesian};
312 unless (defined $inverted) {
313 $z1->set_cartesian([$re1 - $re2, $im1 - $im2]);
317 (ref $z1)->make($re2 - $re1, $im2 - $im1) :
318 (ref $z1)->make($re1 - $re2, $im1 - $im2);
328 my ($z1, $z2, $regular) = @_;
329 if ($z1->{p_dirty} == 0 and ref $z2 and $z2->{p_dirty} == 0) {
330 # if both polar better use polar to avoid rounding errors
331 my ($r1, $t1) = @{$z1->polar};
332 my ($r2, $t2) = @{$z2->polar};
334 if ($t > pi()) { $t -= pit2 }
335 elsif ($t <= -pi()) { $t += pit2 }
336 unless (defined $regular) {
337 $z1->set_polar([$r1 * $r2, $t]);
340 return (ref $z1)->emake($r1 * $r2, $t);
342 my ($x1, $y1) = @{$z1->cartesian};
344 my ($x2, $y2) = @{$z2->cartesian};
345 return (ref $z1)->make($x1*$x2-$y1*$y2, $x1*$y2+$y1*$x2);
347 return (ref $z1)->make($x1*$z2, $y1*$z2);
355 # Die on division by zero.
358 my $mess = "$_[0]: Division by zero.\n";
361 $mess .= "(Because in the definition of $_[0], the divisor ";
362 $mess .= "$_[1] " unless ("$_[1]" eq '0');
368 $mess .= "Died at $up[1] line $up[2].\n";
379 my ($z1, $z2, $inverted) = @_;
380 if ($z1->{p_dirty} == 0 and ref $z2 and $z2->{p_dirty} == 0) {
381 # if both polar better use polar to avoid rounding errors
382 my ($r1, $t1) = @{$z1->polar};
383 my ($r2, $t2) = @{$z2->polar};
386 _divbyzero "$z2/0" if ($r1 == 0);
388 if ($t > pi()) { $t -= pit2 }
389 elsif ($t <= -pi()) { $t += pit2 }
390 return (ref $z1)->emake($r2 / $r1, $t);
392 _divbyzero "$z1/0" if ($r2 == 0);
394 if ($t > pi()) { $t -= pit2 }
395 elsif ($t <= -pi()) { $t += pit2 }
396 return (ref $z1)->emake($r1 / $r2, $t);
401 ($x2, $y2) = @{$z1->cartesian};
402 $d = $x2*$x2 + $y2*$y2;
403 _divbyzero "$z2/0" if $d == 0;
404 return (ref $z1)->make(($x2*$z2)/$d, -($y2*$z2)/$d);
406 my ($x1, $y1) = @{$z1->cartesian};
408 ($x2, $y2) = @{$z2->cartesian};
409 $d = $x2*$x2 + $y2*$y2;
410 _divbyzero "$z1/0" if $d == 0;
411 my $u = ($x1*$x2 + $y1*$y2)/$d;
412 my $v = ($y1*$x2 - $x1*$y2)/$d;
413 return (ref $z1)->make($u, $v);
415 _divbyzero "$z1/0" if $z2 == 0;
416 return (ref $z1)->make($x1/$z2, $y1/$z2);
425 # Computes z1**z2 = exp(z2 * log z1)).
428 my ($z1, $z2, $inverted) = @_;
430 return 1 if $z1 == 0 || $z2 == 1;
431 return 0 if $z2 == 0 && Re($z1) > 0;
433 return 1 if $z2 == 0 || $z1 == 1;
434 return 0 if $z1 == 0 && Re($z2) > 0;
436 my $w = $inverted ? &exp($z1 * &log($z2))
437 : &exp($z2 * &log($z1));
438 # If both arguments cartesian, return cartesian, else polar.
439 return $z1->{c_dirty} == 0 &&
440 (not ref $z2 or $z2->{c_dirty} == 0) ?
441 cplx(@{$w->cartesian}) : $w;
447 # Computes z1 <=> z2.
448 # Sorts on the real part first, then on the imaginary part. Thus 2-4i < 3+8i.
451 my ($z1, $z2, $inverted) = @_;
452 my ($re1, $im1) = ref $z1 ? @{$z1->cartesian} : ($z1, 0);
453 my ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
454 my $sgn = $inverted ? -1 : 1;
455 return $sgn * ($re1 <=> $re2) if $re1 != $re2;
456 return $sgn * ($im1 <=> $im2);
464 # (Required in addition to spaceship() because of NaNs.)
466 my ($z1, $z2, $inverted) = @_;
467 my ($re1, $im1) = ref $z1 ? @{$z1->cartesian} : ($z1, 0);
468 my ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
469 return $re1 == $re2 && $im1 == $im2 ? 1 : 0;
480 my ($r, $t) = @{$z->polar};
481 $t = ($t <= 0) ? $t + pi : $t - pi;
482 return (ref $z)->emake($r, $t);
484 my ($re, $im) = @{$z->cartesian};
485 return (ref $z)->make(-$re, -$im);
491 # Compute complex's conjugate.
496 my ($r, $t) = @{$z->polar};
497 return (ref $z)->emake($r, -$t);
499 my ($re, $im) = @{$z->cartesian};
500 return (ref $z)->make($re, -$im);
506 # Compute or set complex's norm (rho).
514 return CORE::abs($z);
518 $z->{'polar'} = [ $rho, ${$z->polar}[1] ];
523 return ${$z->polar}[0];
530 if ($$theta > pi()) { $$theta -= pit2 }
531 elsif ($$theta <= -pi()) { $$theta += pit2 }
537 # Compute or set complex's argument (theta).
540 my ($z, $theta) = @_;
541 return $z unless ref $z;
542 if (defined $theta) {
544 $z->{'polar'} = [ ${$z->polar}[0], $theta ];
548 $theta = ${$z->polar}[1];
559 # It is quite tempting to use wantarray here so that in list context
560 # sqrt() would return the two solutions. This, however, would
563 # print "sqrt(z) = ", sqrt($z), "\n";
565 # The two values would be printed side by side without no intervening
566 # whitespace, quite confusing.
567 # Therefore if you want the two solutions use the root().
571 my ($re, $im) = ref $z ? @{$z->cartesian} : ($z, 0);
572 return $re < 0 ? cplx(0, CORE::sqrt(-$re)) : CORE::sqrt($re)
574 my ($r, $t) = @{$z->polar};
575 return (ref $z)->emake(CORE::sqrt($r), $t/2);
581 # Compute cbrt(z) (cubic root).
583 # Why are we not returning three values? The same answer as for sqrt().
588 -CORE::exp(CORE::log(-$z)/3) :
589 ($z > 0 ? CORE::exp(CORE::log($z)/3): 0)
591 my ($r, $t) = @{$z->polar};
593 return (ref $z)->emake(CORE::exp(CORE::log($r)/3), $t/3);
602 my $mess = "Root $_[0] illegal, root rank must be positive integer.\n";
606 $mess .= "Died at $up[1] line $up[2].\n";
614 # Computes all nth root for z, returning an array whose size is n.
615 # `n' must be a positive integer.
617 # The roots are given by (for k = 0..n-1):
619 # z^(1/n) = r^(1/n) (cos ((t+2 k pi)/n) + i sin ((t+2 k pi)/n))
623 _rootbad($n) if ($n < 1 or int($n) != $n);
624 my ($r, $t) = ref $z ?
625 @{$z->polar} : (CORE::abs($z), $z >= 0 ? 0 : pi);
628 my $theta_inc = pit2 / $n;
629 my $rho = $r ** (1/$n);
631 my $cartesian = ref $z && $z->{c_dirty} == 0;
632 for ($k = 0, $theta = $t / $n; $k < $n; $k++, $theta += $theta_inc) {
633 my $w = cplxe($rho, $theta);
634 # Yes, $cartesian is loop invariant.
635 push @root, $cartesian ? cplx(@{$w->cartesian}) : $w;
643 # Return or set Re(z).
647 return $z unless ref $z;
649 $z->{'cartesian'} = [ $Re, ${$z->cartesian}[1] ];
653 return ${$z->cartesian}[0];
660 # Return or set Im(z).
664 return $z unless ref $z;
666 $z->{'cartesian'} = [ ${$z->cartesian}[0], $Im ];
670 return ${$z->cartesian}[1];
677 # Return or set rho(w).
680 Math::Complex::abs(@_);
686 # Return or set theta(w).
689 Math::Complex::arg(@_);
699 my ($x, $y) = @{$z->cartesian};
700 return (ref $z)->emake(CORE::exp($x), $y);
706 # Die on logarithm of zero.
709 my $mess = "$_[0]: Logarithm of zero.\n";
712 $mess .= "(Because in the definition of $_[0], the argument ";
713 $mess .= "$_[1] " unless ($_[1] eq '0');
719 $mess .= "Died at $up[1] line $up[2].\n";
732 _logofzero("log") if $z == 0;
733 return $z > 0 ? CORE::log($z) : cplx(CORE::log(-$z), pi);
735 my ($r, $t) = @{$z->polar};
736 _logofzero("log") if $r == 0;
737 if ($t > pi()) { $t -= pit2 }
738 elsif ($t <= -pi()) { $t += pit2 }
739 return (ref $z)->make(CORE::log($r), $t);
747 sub ln { Math::Complex::log(@_) }
756 return Math::Complex::log($_[0]) * uplog10;
762 # Compute logn(z,n) = log(z) / log(n)
766 $z = cplx($z, 0) unless ref $z;
767 my $logn = $LOGN{$n};
768 $logn = $LOGN{$n} = CORE::log($n) unless defined $logn; # Cache log(n)
769 return &log($z) / $logn;
775 # Compute cos(z) = (exp(iz) + exp(-iz))/2.
779 return CORE::cos($z) unless ref $z;
780 my ($x, $y) = @{$z->cartesian};
781 my $ey = CORE::exp($y);
782 my $sx = CORE::sin($x);
783 my $cx = CORE::cos($x);
784 my $ey_1 = $ey ? 1 / $ey : $Inf;
785 return (ref $z)->make($cx * ($ey + $ey_1)/2,
786 $sx * ($ey_1 - $ey)/2);
792 # Compute sin(z) = (exp(iz) - exp(-iz))/2.
796 return CORE::sin($z) unless ref $z;
797 my ($x, $y) = @{$z->cartesian};
798 my $ey = CORE::exp($y);
799 my $sx = CORE::sin($x);
800 my $cx = CORE::cos($x);
801 my $ey_1 = $ey ? 1 / $ey : $Inf;
802 return (ref $z)->make($sx * ($ey + $ey_1)/2,
803 $cx * ($ey - $ey_1)/2);
809 # Compute tan(z) = sin(z) / cos(z).
814 _divbyzero "tan($z)", "cos($z)" if $cz == 0;
815 return &sin($z) / $cz;
821 # Computes the secant sec(z) = 1 / cos(z).
826 _divbyzero "sec($z)", "cos($z)" if ($cz == 0);
833 # Computes the cosecant csc(z) = 1 / sin(z).
838 _divbyzero "csc($z)", "sin($z)" if ($sz == 0);
847 sub cosec { Math::Complex::csc(@_) }
852 # Computes cot(z) = cos(z) / sin(z).
857 _divbyzero "cot($z)", "sin($z)" if ($sz == 0);
858 return &cos($z) / $sz;
866 sub cotan { Math::Complex::cot(@_) }
871 # Computes the arc cosine acos(z) = -i log(z + sqrt(z*z-1)).
875 return CORE::atan2(CORE::sqrt(1-$z*$z), $z)
876 if (! ref $z) && CORE::abs($z) <= 1;
877 $z = cplx($z, 0) unless ref $z;
878 my ($x, $y) = @{$z->cartesian};
879 return 0 if $x == 1 && $y == 0;
880 my $t1 = CORE::sqrt(($x+1)*($x+1) + $y*$y);
881 my $t2 = CORE::sqrt(($x-1)*($x-1) + $y*$y);
882 my $alpha = ($t1 + $t2)/2;
883 my $beta = ($t1 - $t2)/2;
884 $alpha = 1 if $alpha < 1;
885 if ($beta > 1) { $beta = 1 }
886 elsif ($beta < -1) { $beta = -1 }
887 my $u = CORE::atan2(CORE::sqrt(1-$beta*$beta), $beta);
888 my $v = CORE::log($alpha + CORE::sqrt($alpha*$alpha-1));
889 $v = -$v if $y > 0 || ($y == 0 && $x < -1);
890 return (ref $z)->make($u, $v);
896 # Computes the arc sine asin(z) = -i log(iz + sqrt(1-z*z)).
900 return CORE::atan2($z, CORE::sqrt(1-$z*$z))
901 if (! ref $z) && CORE::abs($z) <= 1;
902 $z = cplx($z, 0) unless ref $z;
903 my ($x, $y) = @{$z->cartesian};
904 return 0 if $x == 0 && $y == 0;
905 my $t1 = CORE::sqrt(($x+1)*($x+1) + $y*$y);
906 my $t2 = CORE::sqrt(($x-1)*($x-1) + $y*$y);
907 my $alpha = ($t1 + $t2)/2;
908 my $beta = ($t1 - $t2)/2;
909 $alpha = 1 if $alpha < 1;
910 if ($beta > 1) { $beta = 1 }
911 elsif ($beta < -1) { $beta = -1 }
912 my $u = CORE::atan2($beta, CORE::sqrt(1-$beta*$beta));
913 my $v = -CORE::log($alpha + CORE::sqrt($alpha*$alpha-1));
914 $v = -$v if $y > 0 || ($y == 0 && $x < -1);
915 return (ref $z)->make($u, $v);
921 # Computes the arc tangent atan(z) = i/2 log((i+z) / (i-z)).
925 return CORE::atan2($z, 1) unless ref $z;
926 my ($x, $y) = ref $z ? @{$z->cartesian} : ($z, 0);
927 return 0 if $x == 0 && $y == 0;
928 _divbyzero "atan(i)" if ( $z == i);
929 _logofzero "atan(-i)" if (-$z == i); # -i is a bad file test...
930 my $log = &log((i + $z) / (i - $z));
937 # Computes the arc secant asec(z) = acos(1 / z).
941 _divbyzero "asec($z)", $z if ($z == 0);
948 # Computes the arc cosecant acsc(z) = asin(1 / z).
952 _divbyzero "acsc($z)", $z if ($z == 0);
961 sub acosec { Math::Complex::acsc(@_) }
966 # Computes the arc cotangent acot(z) = atan(1 / z)
970 _divbyzero "acot(0)" if $z == 0;
971 return ($z >= 0) ? CORE::atan2(1, $z) : CORE::atan2(-1, -$z)
973 _divbyzero "acot(i)" if ($z - i == 0);
974 _logofzero "acot(-i)" if ($z + i == 0);
983 sub acotan { Math::Complex::acot(@_) }
988 # Computes the hyperbolic cosine cosh(z) = (exp(z) + exp(-z))/2.
995 return $ex ? ($ex + 1/$ex)/2 : $Inf;
997 my ($x, $y) = @{$z->cartesian};
999 my $ex_1 = $ex ? 1 / $ex : $Inf;
1000 return (ref $z)->make(CORE::cos($y) * ($ex + $ex_1)/2,
1001 CORE::sin($y) * ($ex - $ex_1)/2);
1007 # Computes the hyperbolic sine sinh(z) = (exp(z) - exp(-z))/2.
1013 return 0 if $z == 0;
1014 $ex = CORE::exp($z);
1015 return $ex ? ($ex - 1/$ex)/2 : "-$Inf";
1017 my ($x, $y) = @{$z->cartesian};
1018 my $cy = CORE::cos($y);
1019 my $sy = CORE::sin($y);
1020 $ex = CORE::exp($x);
1021 my $ex_1 = $ex ? 1 / $ex : $Inf;
1022 return (ref $z)->make(CORE::cos($y) * ($ex - $ex_1)/2,
1023 CORE::sin($y) * ($ex + $ex_1)/2);
1029 # Computes the hyperbolic tangent tanh(z) = sinh(z) / cosh(z).
1034 _divbyzero "tanh($z)", "cosh($z)" if ($cz == 0);
1035 return sinh($z) / $cz;
1041 # Computes the hyperbolic secant sech(z) = 1 / cosh(z).
1046 _divbyzero "sech($z)", "cosh($z)" if ($cz == 0);
1053 # Computes the hyperbolic cosecant csch(z) = 1 / sinh(z).
1058 _divbyzero "csch($z)", "sinh($z)" if ($sz == 0);
1067 sub cosech { Math::Complex::csch(@_) }
1072 # Computes the hyperbolic cotangent coth(z) = cosh(z) / sinh(z).
1077 _divbyzero "coth($z)", "sinh($z)" if $sz == 0;
1078 return cosh($z) / $sz;
1086 sub cotanh { Math::Complex::coth(@_) }
1091 # Computes the arc hyperbolic cosine acosh(z) = log(z + sqrt(z*z-1)).
1098 my ($re, $im) = @{$z->cartesian};
1100 return CORE::log($re + CORE::sqrt($re*$re - 1))
1102 return cplx(0, CORE::atan2(CORE::sqrt(1 - $re*$re), $re))
1103 if CORE::abs($re) < 1;
1105 my $t = &sqrt($z * $z - 1) + $z;
1106 # Try Taylor if looking bad (this usually means that
1107 # $z was large negative, therefore the sqrt is really
1108 # close to abs(z), summing that with z...)
1109 $t = 1/(2 * $z) - 1/(8 * $z**3) + 1/(16 * $z**5) - 5/(128 * $z**7)
1112 $u->Im(-$u->Im) if $re < 0 && $im == 0;
1113 return $re < 0 ? -$u : $u;
1119 # Computes the arc hyperbolic sine asinh(z) = log(z + sqrt(z*z+1))
1124 my $t = $z + CORE::sqrt($z*$z + 1);
1125 return CORE::log($t) if $t;
1127 my $t = &sqrt($z * $z + 1) + $z;
1128 # Try Taylor if looking bad (this usually means that
1129 # $z was large negative, therefore the sqrt is really
1130 # close to abs(z), summing that with z...)
1131 $t = 1/(2 * $z) - 1/(8 * $z**3) + 1/(16 * $z**5) - 5/(128 * $z**7)
1139 # Computes the arc hyperbolic tangent atanh(z) = 1/2 log((1+z) / (1-z)).
1144 return CORE::log((1 + $z)/(1 - $z))/2 if CORE::abs($z) < 1;
1147 _divbyzero 'atanh(1)', "1 - $z" if (1 - $z == 0);
1148 _logofzero 'atanh(-1)' if (1 + $z == 0);
1149 return 0.5 * &log((1 + $z) / (1 - $z));
1155 # Computes the hyperbolic arc secant asech(z) = acosh(1 / z).
1159 _divbyzero 'asech(0)', "$z" if ($z == 0);
1160 return acosh(1 / $z);
1166 # Computes the hyperbolic arc cosecant acsch(z) = asinh(1 / z).
1170 _divbyzero 'acsch(0)', $z if ($z == 0);
1171 return asinh(1 / $z);
1177 # Alias for acosh().
1179 sub acosech { Math::Complex::acsch(@_) }
1184 # Computes the arc hyperbolic cotangent acoth(z) = 1/2 log((1+z) / (z-1)).
1188 _divbyzero 'acoth(0)' if ($z == 0);
1190 return CORE::log(($z + 1)/($z - 1))/2 if CORE::abs($z) > 1;
1193 _divbyzero 'acoth(1)', "$z - 1" if ($z - 1 == 0);
1194 _logofzero 'acoth(-1)', "1 + $z" if (1 + $z == 0);
1195 return &log((1 + $z) / ($z - 1)) / 2;
1203 sub acotanh { Math::Complex::acoth(@_) }
1208 # Compute atan(z1/z2).
1211 my ($z1, $z2, $inverted) = @_;
1212 my ($re1, $im1, $re2, $im2);
1214 ($re1, $im1) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
1215 ($re2, $im2) = @{$z1->cartesian};
1217 ($re1, $im1) = @{$z1->cartesian};
1218 ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
1221 return CORE::atan2($re1, $re2) if $im1 == 0;
1222 return ($im1<=>0) * pip2 if $re2 == 0;
1224 my $w = atan($z1/$z2);
1225 my ($u, $v) = ref $w ? @{$w->cartesian} : ($w, 0);
1226 $u += pi if $re2 < 0;
1227 $u -= pit2 if $u > pi;
1228 return cplx($u, $v);
1235 # Set (get if no argument) the display format for all complex numbers that
1236 # don't happen to have overridden it via ->display_format
1238 # When called as an object method, this actually sets the display format for
1239 # the current object.
1241 # Valid object formats are 'c' and 'p' for cartesian and polar. The first
1242 # letter is used actually, so the type can be fully spelled out for clarity.
1244 sub display_format {
1246 my %display_format = %DISPLAY_FORMAT;
1248 if (ref $self) { # Called as an object method
1249 if (exists $self->{display_format}) {
1250 my %obj = %{$self->{display_format}};
1251 @display_format{keys %obj} = values %obj;
1254 $display_format{style} = shift;
1257 @display_format{keys %new} = values %new;
1259 } else { # Called as a class method
1261 $display_format{style} = $self;
1264 @display_format{keys %new} = values %new;
1269 if (defined $self) {
1270 $self->{display_format} = { %display_format };
1273 %{$self->{display_format}} :
1274 $self->{display_format}->{style};
1277 %DISPLAY_FORMAT = %display_format;
1281 $DISPLAY_FORMAT{style};
1287 # Show nicely formatted complex number under its cartesian or polar form,
1288 # depending on the current display format:
1290 # . If a specific display format has been recorded for this object, use it.
1291 # . Otherwise, use the generic current default for all complex numbers,
1292 # which is a package global variable.
1297 my $style = $z->display_format;
1299 $style = $DISPLAY_FORMAT{style} unless defined $style;
1301 return $z->stringify_polar if $style =~ /^p/i;
1302 return $z->stringify_cartesian;
1306 # ->stringify_cartesian
1308 # Stringify as a cartesian representation 'a+bi'.
1310 sub stringify_cartesian {
1312 my ($x, $y) = @{$z->cartesian};
1315 my %format = $z->display_format;
1316 my $format = $format{format};
1319 if ($x =~ /^NaN[QS]?$/i) {
1322 if ($x =~ /^-?$Inf$/oi) {
1325 $re = defined $format ? sprintf($format, $x) : $x;
1333 if ($y =~ /^(NaN[QS]?)$/i) {
1336 if ($y =~ /^-?$Inf$/oi) {
1341 sprintf($format, $y) :
1342 ($y == 1 ? "" : ($y == -1 ? "-" : $y));
1355 } elsif ($y > 0 || $im =~ /^NaN[QS]?i$/i) {
1356 $str .= "+" if defined $re;
1359 } elsif (!defined $re) {
1370 # Stringify as a polar representation '[r,t]'.
1372 sub stringify_polar {
1374 my ($r, $t) = @{$z->polar};
1377 my %format = $z->display_format;
1378 my $format = $format{format};
1380 if ($t =~ /^NaN[QS]?$/i || $t =~ /^-?$Inf$/oi) {
1382 } elsif ($t == pi) {
1384 } elsif ($r == 0 || $t == 0) {
1385 $theta = defined $format ? sprintf($format, $t) : $t;
1388 return "[$r,$theta]" if defined $theta;
1391 # Try to identify pi/n and friends.
1394 $t -= int(CORE::abs($t) / pit2) * pit2;
1396 if ($format{polar_pretty_print} && $t) {
1400 if ($b =~ /^-?\d+$/) {
1401 $b = $b < 0 ? "-" : "" if CORE::abs($b) == 1;
1402 $theta = "${b}pi/$a";
1408 if (defined $format) {
1409 $r = sprintf($format, $r);
1410 $theta = sprintf($format, $theta) unless defined $theta;
1412 $theta = $t unless defined $theta;
1415 return "[$r,$theta]";
1423 Math::Complex - complex numbers and associated mathematical functions
1429 $z = Math::Complex->make(5, 6);
1431 $j = cplxe(1, 2*pi/3);
1435 This package lets you create and manipulate complex numbers. By default,
1436 I<Perl> limits itself to real numbers, but an extra C<use> statement brings
1437 full complex support, along with a full set of mathematical functions
1438 typically associated with and/or extended to complex numbers.
1440 If you wonder what complex numbers are, they were invented to be able to solve
1441 the following equation:
1445 and by definition, the solution is noted I<i> (engineers use I<j> instead since
1446 I<i> usually denotes an intensity, but the name does not matter). The number
1447 I<i> is a pure I<imaginary> number.
1449 The arithmetics with pure imaginary numbers works just like you would expect
1450 it with real numbers... you just have to remember that
1456 5i + 7i = i * (5 + 7) = 12i
1457 4i - 3i = i * (4 - 3) = i
1462 Complex numbers are numbers that have both a real part and an imaginary
1463 part, and are usually noted:
1467 where C<a> is the I<real> part and C<b> is the I<imaginary> part. The
1468 arithmetic with complex numbers is straightforward. You have to
1469 keep track of the real and the imaginary parts, but otherwise the
1470 rules used for real numbers just apply:
1472 (4 + 3i) + (5 - 2i) = (4 + 5) + i(3 - 2) = 9 + i
1473 (2 + i) * (4 - i) = 2*4 + 4i -2i -i*i = 8 + 2i + 1 = 9 + 2i
1475 A graphical representation of complex numbers is possible in a plane
1476 (also called the I<complex plane>, but it's really a 2D plane).
1481 is the point whose coordinates are (a, b). Actually, it would
1482 be the vector originating from (0, 0) to (a, b). It follows that the addition
1483 of two complex numbers is a vectorial addition.
1485 Since there is a bijection between a point in the 2D plane and a complex
1486 number (i.e. the mapping is unique and reciprocal), a complex number
1487 can also be uniquely identified with polar coordinates:
1491 where C<rho> is the distance to the origin, and C<theta> the angle between
1492 the vector and the I<x> axis. There is a notation for this using the
1493 exponential form, which is:
1495 rho * exp(i * theta)
1497 where I<i> is the famous imaginary number introduced above. Conversion
1498 between this form and the cartesian form C<a + bi> is immediate:
1500 a = rho * cos(theta)
1501 b = rho * sin(theta)
1503 which is also expressed by this formula:
1505 z = rho * exp(i * theta) = rho * (cos theta + i * sin theta)
1507 In other words, it's the projection of the vector onto the I<x> and I<y>
1508 axes. Mathematicians call I<rho> the I<norm> or I<modulus> and I<theta>
1509 the I<argument> of the complex number. The I<norm> of C<z> will be
1512 The polar notation (also known as the trigonometric
1513 representation) is much more handy for performing multiplications and
1514 divisions of complex numbers, whilst the cartesian notation is better
1515 suited for additions and subtractions. Real numbers are on the I<x>
1516 axis, and therefore I<theta> is zero or I<pi>.
1518 All the common operations that can be performed on a real number have
1519 been defined to work on complex numbers as well, and are merely
1520 I<extensions> of the operations defined on real numbers. This means
1521 they keep their natural meaning when there is no imaginary part, provided
1522 the number is within their definition set.
1524 For instance, the C<sqrt> routine which computes the square root of
1525 its argument is only defined for non-negative real numbers and yields a
1526 non-negative real number (it is an application from B<R+> to B<R+>).
1527 If we allow it to return a complex number, then it can be extended to
1528 negative real numbers to become an application from B<R> to B<C> (the
1529 set of complex numbers):
1531 sqrt(x) = x >= 0 ? sqrt(x) : sqrt(-x)*i
1533 It can also be extended to be an application from B<C> to B<C>,
1534 whilst its restriction to B<R> behaves as defined above by using
1535 the following definition:
1537 sqrt(z = [r,t]) = sqrt(r) * exp(i * t/2)
1539 Indeed, a negative real number can be noted C<[x,pi]> (the modulus
1540 I<x> is always non-negative, so C<[x,pi]> is really C<-x>, a negative
1541 number) and the above definition states that
1543 sqrt([x,pi]) = sqrt(x) * exp(i*pi/2) = [sqrt(x),pi/2] = sqrt(x)*i
1545 which is exactly what we had defined for negative real numbers above.
1546 The C<sqrt> returns only one of the solutions: if you want the both,
1547 use the C<root> function.
1549 All the common mathematical functions defined on real numbers that
1550 are extended to complex numbers share that same property of working
1551 I<as usual> when the imaginary part is zero (otherwise, it would not
1552 be called an extension, would it?).
1554 A I<new> operation possible on a complex number that is
1555 the identity for real numbers is called the I<conjugate>, and is noted
1556 with an horizontal bar above the number, or C<~z> here.
1563 z * ~z = (a + bi) * (a - bi) = a*a + b*b
1565 We saw that the norm of C<z> was noted C<abs(z)> and was defined as the
1566 distance to the origin, also known as:
1568 rho = abs(z) = sqrt(a*a + b*b)
1572 z * ~z = abs(z) ** 2
1574 If z is a pure real number (i.e. C<b == 0>), then the above yields:
1578 which is true (C<abs> has the regular meaning for real number, i.e. stands
1579 for the absolute value). This example explains why the norm of C<z> is
1580 noted C<abs(z)>: it extends the C<abs> function to complex numbers, yet
1581 is the regular C<abs> we know when the complex number actually has no
1582 imaginary part... This justifies I<a posteriori> our use of the C<abs>
1583 notation for the norm.
1587 Given the following notations:
1589 z1 = a + bi = r1 * exp(i * t1)
1590 z2 = c + di = r2 * exp(i * t2)
1591 z = <any complex or real number>
1593 the following (overloaded) operations are supported on complex numbers:
1595 z1 + z2 = (a + c) + i(b + d)
1596 z1 - z2 = (a - c) + i(b - d)
1597 z1 * z2 = (r1 * r2) * exp(i * (t1 + t2))
1598 z1 / z2 = (r1 / r2) * exp(i * (t1 - t2))
1599 z1 ** z2 = exp(z2 * log z1)
1601 abs(z) = r1 = sqrt(a*a + b*b)
1602 sqrt(z) = sqrt(r1) * exp(i * t/2)
1603 exp(z) = exp(a) * exp(i * b)
1604 log(z) = log(r1) + i*t
1605 sin(z) = 1/2i (exp(i * z1) - exp(-i * z))
1606 cos(z) = 1/2 (exp(i * z1) + exp(-i * z))
1607 atan2(z1, z2) = atan(z1/z2)
1609 The following extra operations are supported on both real and complex
1617 cbrt(z) = z ** (1/3)
1618 log10(z) = log(z) / log(10)
1619 logn(z, n) = log(z) / log(n)
1621 tan(z) = sin(z) / cos(z)
1627 asin(z) = -i * log(i*z + sqrt(1-z*z))
1628 acos(z) = -i * log(z + i*sqrt(1-z*z))
1629 atan(z) = i/2 * log((i+z) / (i-z))
1631 acsc(z) = asin(1 / z)
1632 asec(z) = acos(1 / z)
1633 acot(z) = atan(1 / z) = -i/2 * log((i+z) / (z-i))
1635 sinh(z) = 1/2 (exp(z) - exp(-z))
1636 cosh(z) = 1/2 (exp(z) + exp(-z))
1637 tanh(z) = sinh(z) / cosh(z) = (exp(z) - exp(-z)) / (exp(z) + exp(-z))
1639 csch(z) = 1 / sinh(z)
1640 sech(z) = 1 / cosh(z)
1641 coth(z) = 1 / tanh(z)
1643 asinh(z) = log(z + sqrt(z*z+1))
1644 acosh(z) = log(z + sqrt(z*z-1))
1645 atanh(z) = 1/2 * log((1+z) / (1-z))
1647 acsch(z) = asinh(1 / z)
1648 asech(z) = acosh(1 / z)
1649 acoth(z) = atanh(1 / z) = 1/2 * log((1+z) / (z-1))
1651 I<arg>, I<abs>, I<log>, I<csc>, I<cot>, I<acsc>, I<acot>, I<csch>,
1652 I<coth>, I<acosech>, I<acotanh>, have aliases I<rho>, I<theta>, I<ln>,
1653 I<cosec>, I<cotan>, I<acosec>, I<acotan>, I<cosech>, I<cotanh>,
1654 I<acosech>, I<acotanh>, respectively. C<Re>, C<Im>, C<arg>, C<abs>,
1655 C<rho>, and C<theta> can be used also also mutators. The C<cbrt>
1656 returns only one of the solutions: if you want all three, use the
1659 The I<root> function is available to compute all the I<n>
1660 roots of some complex, where I<n> is a strictly positive integer.
1661 There are exactly I<n> such roots, returned as a list. Getting the
1662 number mathematicians call C<j> such that:
1666 is a simple matter of writing:
1668 $j = ((root(1, 3))[1];
1670 The I<k>th root for C<z = [r,t]> is given by:
1672 (root(z, n))[k] = r**(1/n) * exp(i * (t + 2*k*pi)/n)
1674 The I<spaceship> comparison operator, E<lt>=E<gt>, is also defined. In
1675 order to ensure its restriction to real numbers is conform to what you
1676 would expect, the comparison is run on the real part of the complex
1677 number first, and imaginary parts are compared only when the real
1682 To create a complex number, use either:
1684 $z = Math::Complex->make(3, 4);
1687 if you know the cartesian form of the number, or
1691 if you like. To create a number using the polar form, use either:
1693 $z = Math::Complex->emake(5, pi/3);
1694 $x = cplxe(5, pi/3);
1696 instead. The first argument is the modulus, the second is the angle
1697 (in radians, the full circle is 2*pi). (Mnemonic: C<e> is used as a
1698 notation for complex numbers in the polar form).
1700 It is possible to write:
1702 $x = cplxe(-3, pi/4);
1704 but that will be silently converted into C<[3,-3pi/4]>, since the
1705 modulus must be non-negative (it represents the distance to the origin
1706 in the complex plane).
1708 It is also possible to have a complex number as either argument of
1709 either the C<make> or C<emake>: the appropriate component of
1710 the argument will be used.
1715 =head1 STRINGIFICATION
1717 When printed, a complex number is usually shown under its cartesian
1718 style I<a+bi>, but there are legitimate cases where the polar style
1719 I<[r,t]> is more appropriate.
1721 By calling the class method C<Math::Complex::display_format> and
1722 supplying either C<"polar"> or C<"cartesian"> as an argument, you
1723 override the default display style, which is C<"cartesian">. Not
1724 supplying any argument returns the current settings.
1726 This default can be overridden on a per-number basis by calling the
1727 C<display_format> method instead. As before, not supplying any argument
1728 returns the current display style for this number. Otherwise whatever you
1729 specify will be the new display style for I<this> particular number.
1735 Math::Complex::display_format('polar');
1736 $j = (root(1, 3))[1];
1737 print "j = $j\n"; # Prints "j = [1,2pi/3]"
1738 $j->display_format('cartesian');
1739 print "j = $j\n"; # Prints "j = -0.5+0.866025403784439i"
1741 The polar style attempts to emphasize arguments like I<k*pi/n>
1742 (where I<n> is a positive integer and I<k> an integer within [-9, +9]),
1743 this is called I<polar pretty-printing>.
1745 =head2 CHANGED IN PERL 5.6
1747 The C<display_format> class method and the corresponding
1748 C<display_format> object method can now be called using
1749 a parameter hash instead of just a one parameter.
1751 The old display format style, which can have values C<"cartesian"> or
1752 C<"polar">, can be changed using the C<"style"> parameter.
1754 $j->display_format(style => "polar");
1756 The one parameter calling convention also still works.
1758 $j->display_format("polar");
1760 There are two new display parameters.
1762 The first one is C<"format">, which is a sprintf()-style format string
1763 to be used for both numeric parts of the complex number(s). The is
1764 somewhat system-dependent but most often it corresponds to C<"%.15g">.
1765 You can revert to the default by setting the C<format> to C<undef>.
1767 # the $j from the above example
1769 $j->display_format('format' => '%.5f');
1770 print "j = $j\n"; # Prints "j = -0.50000+0.86603i"
1771 $j->display_format('format' => undef);
1772 print "j = $j\n"; # Prints "j = -0.5+0.86603i"
1774 Notice that this affects also the return values of the
1775 C<display_format> methods: in list context the whole parameter hash
1776 will be returned, as opposed to only the style parameter value.
1777 This is a potential incompatibility with earlier versions if you
1778 have been calling the C<display_format> method in list context.
1780 The second new display parameter is C<"polar_pretty_print">, which can
1781 be set to true or false, the default being true. See the previous
1782 section for what this means.
1786 Thanks to overloading, the handling of arithmetics with complex numbers
1787 is simple and almost transparent.
1789 Here are some examples:
1793 $j = cplxe(1, 2*pi/3); # $j ** 3 == 1
1794 print "j = $j, j**3 = ", $j ** 3, "\n";
1795 print "1 + j + j**2 = ", 1 + $j + $j**2, "\n";
1797 $z = -16 + 0*i; # Force it to be a complex
1798 print "sqrt($z) = ", sqrt($z), "\n";
1800 $k = exp(i * 2*pi/3);
1801 print "$j - $k = ", $j - $k, "\n";
1803 $z->Re(3); # Re, Im, arg, abs,
1804 $j->arg(2); # (the last two aka rho, theta)
1805 # can be used also as mutators.
1807 =head1 ERRORS DUE TO DIVISION BY ZERO OR LOGARITHM OF ZERO
1809 The division (/) and the following functions
1815 atanh asech acsch acoth
1817 cannot be computed for all arguments because that would mean dividing
1818 by zero or taking logarithm of zero. These situations cause fatal
1819 runtime errors looking like this
1821 cot(0): Division by zero.
1822 (Because in the definition of cot(0), the divisor sin(0) is 0)
1827 atanh(-1): Logarithm of zero.
1830 For the C<csc>, C<cot>, C<asec>, C<acsc>, C<acot>, C<csch>, C<coth>,
1831 C<asech>, C<acsch>, the argument cannot be C<0> (zero). For the the
1832 logarithmic functions and the C<atanh>, C<acoth>, the argument cannot
1833 be C<1> (one). For the C<atanh>, C<acoth>, the argument cannot be
1834 C<-1> (minus one). For the C<atan>, C<acot>, the argument cannot be
1835 C<i> (the imaginary unit). For the C<atan>, C<acoth>, the argument
1836 cannot be C<-i> (the negative imaginary unit). For the C<tan>,
1837 C<sec>, C<tanh>, the argument cannot be I<pi/2 + k * pi>, where I<k>
1840 Note that because we are operating on approximations of real numbers,
1841 these errors can happen when merely `too close' to the singularities
1844 =head1 ERRORS DUE TO INDIGESTIBLE ARGUMENTS
1846 The C<make> and C<emake> accept both real and complex arguments.
1847 When they cannot recognize the arguments they will die with error
1848 messages like the following
1850 Math::Complex::make: Cannot take real part of ...
1851 Math::Complex::make: Cannot take real part of ...
1852 Math::Complex::emake: Cannot take rho of ...
1853 Math::Complex::emake: Cannot take theta of ...
1857 Saying C<use Math::Complex;> exports many mathematical routines in the
1858 caller environment and even overrides some (C<sqrt>, C<log>).
1859 This is construed as a feature by the Authors, actually... ;-)
1861 All routines expect to be given real or complex numbers. Don't attempt to
1862 use BigFloat, since Perl has currently no rule to disambiguate a '+'
1863 operation (for instance) between two overloaded entities.
1865 In Cray UNICOS there is some strange numerical instability that results
1866 in root(), cos(), sin(), cosh(), sinh(), losing accuracy fast. Beware.
1867 The bug may be in UNICOS math libs, in UNICOS C compiler, in Math::Complex.
1868 Whatever it is, it does not manifest itself anywhere else where Perl runs.
1872 Raphael Manfredi <F<Raphael_Manfredi@pobox.com>> and
1873 Jarkko Hietaniemi <F<jhi@iki.fi>>.
1875 Extensive patches by Daniel S. Lewart <F<d-lewart@uiuc.edu>>.