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 my ($x, $y) = ref $z ? @{$z->cartesian} : ($z, 0);
878 return 0 if $x == 1 && $y == 0;
879 my $t1 = CORE::sqrt(($x+1)*($x+1) + $y*$y);
880 my $t2 = CORE::sqrt(($x-1)*($x-1) + $y*$y);
881 my $alpha = ($t1 + $t2)/2;
882 my $beta = ($t1 - $t2)/2;
883 $alpha = 1 if $alpha < 1;
884 if ($beta > 1) { $beta = 1 }
885 elsif ($beta < -1) { $beta = -1 }
886 my $u = CORE::atan2(CORE::sqrt(1-$beta*$beta), $beta);
887 my $v = CORE::log($alpha + CORE::sqrt($alpha*$alpha-1));
888 $v = -$v if $y > 0 || ($y == 0 && $x < -1);
889 return __PACKAGE__->make($u, $v);
895 # Computes the arc sine asin(z) = -i log(iz + sqrt(1-z*z)).
899 return CORE::atan2($z, CORE::sqrt(1-$z*$z))
900 if (! ref $z) && CORE::abs($z) <= 1;
901 my ($x, $y) = ref $z ? @{$z->cartesian} : ($z, 0);
902 return 0 if $x == 0 && $y == 0;
903 my $t1 = CORE::sqrt(($x+1)*($x+1) + $y*$y);
904 my $t2 = CORE::sqrt(($x-1)*($x-1) + $y*$y);
905 my $alpha = ($t1 + $t2)/2;
906 my $beta = ($t1 - $t2)/2;
907 $alpha = 1 if $alpha < 1;
908 if ($beta > 1) { $beta = 1 }
909 elsif ($beta < -1) { $beta = -1 }
910 my $u = CORE::atan2($beta, CORE::sqrt(1-$beta*$beta));
911 my $v = -CORE::log($alpha + CORE::sqrt($alpha*$alpha-1));
912 $v = -$v if $y > 0 || ($y == 0 && $x < -1);
913 return __PACKAGE__->make($u, $v);
919 # Computes the arc tangent atan(z) = i/2 log((i+z) / (i-z)).
923 return CORE::atan2($z, 1) unless ref $z;
924 my ($x, $y) = ref $z ? @{$z->cartesian} : ($z, 0);
925 return 0 if $x == 0 && $y == 0;
926 _divbyzero "atan(i)" if ( $z == i);
927 _logofzero "atan(-i)" if (-$z == i); # -i is a bad file test...
928 my $log = &log((i + $z) / (i - $z));
935 # Computes the arc secant asec(z) = acos(1 / z).
939 _divbyzero "asec($z)", $z if ($z == 0);
946 # Computes the arc cosecant acsc(z) = asin(1 / z).
950 _divbyzero "acsc($z)", $z if ($z == 0);
959 sub acosec { Math::Complex::acsc(@_) }
964 # Computes the arc cotangent acot(z) = atan(1 / z)
968 _divbyzero "acot(0)" if $z == 0;
969 return ($z >= 0) ? CORE::atan2(1, $z) : CORE::atan2(-1, -$z)
971 _divbyzero "acot(i)" if ($z - i == 0);
972 _logofzero "acot(-i)" if ($z + i == 0);
981 sub acotan { Math::Complex::acot(@_) }
986 # Computes the hyperbolic cosine cosh(z) = (exp(z) + exp(-z))/2.
993 return $ex ? ($ex + 1/$ex)/2 : $Inf;
995 my ($x, $y) = @{$z->cartesian};
997 my $ex_1 = $ex ? 1 / $ex : $Inf;
998 return (ref $z)->make(CORE::cos($y) * ($ex + $ex_1)/2,
999 CORE::sin($y) * ($ex - $ex_1)/2);
1005 # Computes the hyperbolic sine sinh(z) = (exp(z) - exp(-z))/2.
1011 return 0 if $z == 0;
1012 $ex = CORE::exp($z);
1013 return $ex ? ($ex - 1/$ex)/2 : "-$Inf";
1015 my ($x, $y) = @{$z->cartesian};
1016 my $cy = CORE::cos($y);
1017 my $sy = CORE::sin($y);
1018 $ex = CORE::exp($x);
1019 my $ex_1 = $ex ? 1 / $ex : $Inf;
1020 return (ref $z)->make(CORE::cos($y) * ($ex - $ex_1)/2,
1021 CORE::sin($y) * ($ex + $ex_1)/2);
1027 # Computes the hyperbolic tangent tanh(z) = sinh(z) / cosh(z).
1032 _divbyzero "tanh($z)", "cosh($z)" if ($cz == 0);
1033 return sinh($z) / $cz;
1039 # Computes the hyperbolic secant sech(z) = 1 / cosh(z).
1044 _divbyzero "sech($z)", "cosh($z)" if ($cz == 0);
1051 # Computes the hyperbolic cosecant csch(z) = 1 / sinh(z).
1056 _divbyzero "csch($z)", "sinh($z)" if ($sz == 0);
1065 sub cosech { Math::Complex::csch(@_) }
1070 # Computes the hyperbolic cotangent coth(z) = cosh(z) / sinh(z).
1075 _divbyzero "coth($z)", "sinh($z)" if $sz == 0;
1076 return cosh($z) / $sz;
1084 sub cotanh { Math::Complex::coth(@_) }
1089 # Computes the arc hyperbolic cosine acosh(z) = log(z + sqrt(z*z-1)).
1096 my ($re, $im) = @{$z->cartesian};
1098 return CORE::log($re + CORE::sqrt($re*$re - 1))
1100 return cplx(0, CORE::atan2(CORE::sqrt(1 - $re*$re), $re))
1101 if CORE::abs($re) < 1;
1103 my $t = &sqrt($z * $z - 1) + $z;
1104 # Try MacLaurin if looking bad.
1105 $t = 1/(2 * $z) - 1/(8 * $z**3) + 1/(16 * $z**5) - 5/(128 * $z**7)
1108 $u->Im(-$u->Im) if $im == 0;
1109 return $re < 0 ? -$u : $u;
1115 # Computes the arc hyperbolic sine asinh(z) = log(z + sqrt(z*z+1))
1120 my $t = $z + CORE::sqrt($z*$z + 1);
1121 return CORE::log($t) if $t;
1123 my $t = &sqrt($z * $z + 1) + $z;
1124 # Try MacLaurin if looking bad.
1125 $t = 1/(2 * $z) - 1/(8 * $z**3) + 1/(16 * $z**5) - 5/(128 * $z**7)
1133 # Computes the arc hyperbolic tangent atanh(z) = 1/2 log((1+z) / (1-z)).
1138 return CORE::log((1 + $z)/(1 - $z))/2 if CORE::abs($z) < 1;
1141 _divbyzero 'atanh(1)', "1 - $z" if (1 - $z == 0);
1142 _logofzero 'atanh(-1)' if (1 + $z == 0);
1143 return 0.5 * &log((1 + $z) / (1 - $z));
1149 # Computes the hyperbolic arc secant asech(z) = acosh(1 / z).
1153 _divbyzero 'asech(0)', "$z" if ($z == 0);
1154 return acosh(1 / $z);
1160 # Computes the hyperbolic arc cosecant acsch(z) = asinh(1 / z).
1164 _divbyzero 'acsch(0)', $z if ($z == 0);
1165 return asinh(1 / $z);
1171 # Alias for acosh().
1173 sub acosech { Math::Complex::acsch(@_) }
1178 # Computes the arc hyperbolic cotangent acoth(z) = 1/2 log((1+z) / (z-1)).
1182 _divbyzero 'acoth(0)' if ($z == 0);
1184 return CORE::log(($z + 1)/($z - 1))/2 if CORE::abs($z) > 1;
1187 _divbyzero 'acoth(1)', "$z - 1" if ($z - 1 == 0);
1188 _logofzero 'acoth(-1)', "1 + $z" if (1 + $z == 0);
1189 return &log((1 + $z) / ($z - 1)) / 2;
1197 sub acotanh { Math::Complex::acoth(@_) }
1202 # Compute atan(z1/z2).
1205 my ($z1, $z2, $inverted) = @_;
1206 my ($re1, $im1, $re2, $im2);
1208 ($re1, $im1) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
1209 ($re2, $im2) = @{$z1->cartesian};
1211 ($re1, $im1) = @{$z1->cartesian};
1212 ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
1215 return CORE::atan2($re1, $re2) if $im1 == 0;
1216 return ($im1<=>0) * pip2 if $re2 == 0;
1218 my $w = atan($z1/$z2);
1219 my ($u, $v) = ref $w ? @{$w->cartesian} : ($w, 0);
1220 $u += pi if $re2 < 0;
1221 $u -= pit2 if $u > pi;
1222 return cplx($u, $v);
1229 # Set (get if no argument) the display format for all complex numbers that
1230 # don't happen to have overridden it via ->display_format
1232 # When called as an object method, this actually sets the display format for
1233 # the current object.
1235 # Valid object formats are 'c' and 'p' for cartesian and polar. The first
1236 # letter is used actually, so the type can be fully spelled out for clarity.
1238 sub display_format {
1240 my %display_format = %DISPLAY_FORMAT;
1242 if (ref $self) { # Called as an object method
1243 if (exists $self->{display_format}) {
1244 my %obj = %{$self->{display_format}};
1245 @display_format{keys %obj} = values %obj;
1248 $display_format{style} = shift;
1251 @display_format{keys %new} = values %new;
1253 } else { # Called as a class method
1255 $display_format{style} = $self;
1258 @display_format{keys %new} = values %new;
1263 if (defined $self) {
1264 $self->{display_format} = { %display_format };
1267 %{$self->{display_format}} :
1268 $self->{display_format}->{style};
1271 %DISPLAY_FORMAT = %display_format;
1275 $DISPLAY_FORMAT{style};
1281 # Show nicely formatted complex number under its cartesian or polar form,
1282 # depending on the current display format:
1284 # . If a specific display format has been recorded for this object, use it.
1285 # . Otherwise, use the generic current default for all complex numbers,
1286 # which is a package global variable.
1291 my $style = $z->display_format;
1293 $style = $DISPLAY_FORMAT{style} unless defined $style;
1295 return $z->stringify_polar if $style =~ /^p/i;
1296 return $z->stringify_cartesian;
1300 # ->stringify_cartesian
1302 # Stringify as a cartesian representation 'a+bi'.
1304 sub stringify_cartesian {
1306 my ($x, $y) = @{$z->cartesian};
1309 my %format = $z->display_format;
1310 my $format = $format{format};
1313 if ($x =~ /^NaN[QS]?$/i) {
1316 if ($x =~ /^-?$Inf$/oi) {
1319 $re = defined $format ? sprintf($format, $x) : $x;
1327 if ($y == 1) { $im = "" }
1328 elsif ($y == -1) { $im = "-" }
1329 elsif ($y =~ /^(NaN[QS]?)$/i) {
1332 if ($y =~ /^-?$Inf$/oi) {
1335 $im = defined $format ? sprintf($format, $y) : $y;
1348 } elsif ($y > 0 || $im =~ /^NaN[QS]?i$/i) {
1349 $str .= "+" if defined $re;
1352 } elsif (!defined $re) {
1363 # Stringify as a polar representation '[r,t]'.
1365 sub stringify_polar {
1367 my ($r, $t) = @{$z->polar};
1370 my %format = $z->display_format;
1371 my $format = $format{format};
1373 if ($t =~ /^NaN[QS]?$/i || $t =~ /^-?$Inf$/oi) {
1375 } elsif ($t == pi) {
1377 } elsif ($r == 0 || $t == 0) {
1378 $theta = defined $format ? sprintf($format, $t) : $t;
1381 return "[$r,$theta]" if defined $theta;
1384 # Try to identify pi/n and friends.
1387 $t -= int(CORE::abs($t) / pit2) * pit2;
1389 if ($format{polar_pretty_print}) {
1393 if (int($b) == $b) {
1394 $b = $b < 0 ? "-" : "" if CORE::abs($b) == 1;
1395 $theta = "${b}pi/$a";
1401 if (defined $format) {
1402 $r = sprintf($format, $r);
1403 $theta = sprintf($format, $theta) unless defined $theta;
1405 $theta = $t unless defined $theta;
1408 return "[$r,$theta]";
1416 Math::Complex - complex numbers and associated mathematical functions
1422 $z = Math::Complex->make(5, 6);
1424 $j = cplxe(1, 2*pi/3);
1428 This package lets you create and manipulate complex numbers. By default,
1429 I<Perl> limits itself to real numbers, but an extra C<use> statement brings
1430 full complex support, along with a full set of mathematical functions
1431 typically associated with and/or extended to complex numbers.
1433 If you wonder what complex numbers are, they were invented to be able to solve
1434 the following equation:
1438 and by definition, the solution is noted I<i> (engineers use I<j> instead since
1439 I<i> usually denotes an intensity, but the name does not matter). The number
1440 I<i> is a pure I<imaginary> number.
1442 The arithmetics with pure imaginary numbers works just like you would expect
1443 it with real numbers... you just have to remember that
1449 5i + 7i = i * (5 + 7) = 12i
1450 4i - 3i = i * (4 - 3) = i
1455 Complex numbers are numbers that have both a real part and an imaginary
1456 part, and are usually noted:
1460 where C<a> is the I<real> part and C<b> is the I<imaginary> part. The
1461 arithmetic with complex numbers is straightforward. You have to
1462 keep track of the real and the imaginary parts, but otherwise the
1463 rules used for real numbers just apply:
1465 (4 + 3i) + (5 - 2i) = (4 + 5) + i(3 - 2) = 9 + i
1466 (2 + i) * (4 - i) = 2*4 + 4i -2i -i*i = 8 + 2i + 1 = 9 + 2i
1468 A graphical representation of complex numbers is possible in a plane
1469 (also called the I<complex plane>, but it's really a 2D plane).
1474 is the point whose coordinates are (a, b). Actually, it would
1475 be the vector originating from (0, 0) to (a, b). It follows that the addition
1476 of two complex numbers is a vectorial addition.
1478 Since there is a bijection between a point in the 2D plane and a complex
1479 number (i.e. the mapping is unique and reciprocal), a complex number
1480 can also be uniquely identified with polar coordinates:
1484 where C<rho> is the distance to the origin, and C<theta> the angle between
1485 the vector and the I<x> axis. There is a notation for this using the
1486 exponential form, which is:
1488 rho * exp(i * theta)
1490 where I<i> is the famous imaginary number introduced above. Conversion
1491 between this form and the cartesian form C<a + bi> is immediate:
1493 a = rho * cos(theta)
1494 b = rho * sin(theta)
1496 which is also expressed by this formula:
1498 z = rho * exp(i * theta) = rho * (cos theta + i * sin theta)
1500 In other words, it's the projection of the vector onto the I<x> and I<y>
1501 axes. Mathematicians call I<rho> the I<norm> or I<modulus> and I<theta>
1502 the I<argument> of the complex number. The I<norm> of C<z> will be
1505 The polar notation (also known as the trigonometric
1506 representation) is much more handy for performing multiplications and
1507 divisions of complex numbers, whilst the cartesian notation is better
1508 suited for additions and subtractions. Real numbers are on the I<x>
1509 axis, and therefore I<theta> is zero or I<pi>.
1511 All the common operations that can be performed on a real number have
1512 been defined to work on complex numbers as well, and are merely
1513 I<extensions> of the operations defined on real numbers. This means
1514 they keep their natural meaning when there is no imaginary part, provided
1515 the number is within their definition set.
1517 For instance, the C<sqrt> routine which computes the square root of
1518 its argument is only defined for non-negative real numbers and yields a
1519 non-negative real number (it is an application from B<R+> to B<R+>).
1520 If we allow it to return a complex number, then it can be extended to
1521 negative real numbers to become an application from B<R> to B<C> (the
1522 set of complex numbers):
1524 sqrt(x) = x >= 0 ? sqrt(x) : sqrt(-x)*i
1526 It can also be extended to be an application from B<C> to B<C>,
1527 whilst its restriction to B<R> behaves as defined above by using
1528 the following definition:
1530 sqrt(z = [r,t]) = sqrt(r) * exp(i * t/2)
1532 Indeed, a negative real number can be noted C<[x,pi]> (the modulus
1533 I<x> is always non-negative, so C<[x,pi]> is really C<-x>, a negative
1534 number) and the above definition states that
1536 sqrt([x,pi]) = sqrt(x) * exp(i*pi/2) = [sqrt(x),pi/2] = sqrt(x)*i
1538 which is exactly what we had defined for negative real numbers above.
1539 The C<sqrt> returns only one of the solutions: if you want the both,
1540 use the C<root> function.
1542 All the common mathematical functions defined on real numbers that
1543 are extended to complex numbers share that same property of working
1544 I<as usual> when the imaginary part is zero (otherwise, it would not
1545 be called an extension, would it?).
1547 A I<new> operation possible on a complex number that is
1548 the identity for real numbers is called the I<conjugate>, and is noted
1549 with an horizontal bar above the number, or C<~z> here.
1556 z * ~z = (a + bi) * (a - bi) = a*a + b*b
1558 We saw that the norm of C<z> was noted C<abs(z)> and was defined as the
1559 distance to the origin, also known as:
1561 rho = abs(z) = sqrt(a*a + b*b)
1565 z * ~z = abs(z) ** 2
1567 If z is a pure real number (i.e. C<b == 0>), then the above yields:
1571 which is true (C<abs> has the regular meaning for real number, i.e. stands
1572 for the absolute value). This example explains why the norm of C<z> is
1573 noted C<abs(z)>: it extends the C<abs> function to complex numbers, yet
1574 is the regular C<abs> we know when the complex number actually has no
1575 imaginary part... This justifies I<a posteriori> our use of the C<abs>
1576 notation for the norm.
1580 Given the following notations:
1582 z1 = a + bi = r1 * exp(i * t1)
1583 z2 = c + di = r2 * exp(i * t2)
1584 z = <any complex or real number>
1586 the following (overloaded) operations are supported on complex numbers:
1588 z1 + z2 = (a + c) + i(b + d)
1589 z1 - z2 = (a - c) + i(b - d)
1590 z1 * z2 = (r1 * r2) * exp(i * (t1 + t2))
1591 z1 / z2 = (r1 / r2) * exp(i * (t1 - t2))
1592 z1 ** z2 = exp(z2 * log z1)
1594 abs(z) = r1 = sqrt(a*a + b*b)
1595 sqrt(z) = sqrt(r1) * exp(i * t/2)
1596 exp(z) = exp(a) * exp(i * b)
1597 log(z) = log(r1) + i*t
1598 sin(z) = 1/2i (exp(i * z1) - exp(-i * z))
1599 cos(z) = 1/2 (exp(i * z1) + exp(-i * z))
1600 atan2(z1, z2) = atan(z1/z2)
1602 The following extra operations are supported on both real and complex
1610 cbrt(z) = z ** (1/3)
1611 log10(z) = log(z) / log(10)
1612 logn(z, n) = log(z) / log(n)
1614 tan(z) = sin(z) / cos(z)
1620 asin(z) = -i * log(i*z + sqrt(1-z*z))
1621 acos(z) = -i * log(z + i*sqrt(1-z*z))
1622 atan(z) = i/2 * log((i+z) / (i-z))
1624 acsc(z) = asin(1 / z)
1625 asec(z) = acos(1 / z)
1626 acot(z) = atan(1 / z) = -i/2 * log((i+z) / (z-i))
1628 sinh(z) = 1/2 (exp(z) - exp(-z))
1629 cosh(z) = 1/2 (exp(z) + exp(-z))
1630 tanh(z) = sinh(z) / cosh(z) = (exp(z) - exp(-z)) / (exp(z) + exp(-z))
1632 csch(z) = 1 / sinh(z)
1633 sech(z) = 1 / cosh(z)
1634 coth(z) = 1 / tanh(z)
1636 asinh(z) = log(z + sqrt(z*z+1))
1637 acosh(z) = log(z + sqrt(z*z-1))
1638 atanh(z) = 1/2 * log((1+z) / (1-z))
1640 acsch(z) = asinh(1 / z)
1641 asech(z) = acosh(1 / z)
1642 acoth(z) = atanh(1 / z) = 1/2 * log((1+z) / (z-1))
1644 I<arg>, I<abs>, I<log>, I<csc>, I<cot>, I<acsc>, I<acot>, I<csch>,
1645 I<coth>, I<acosech>, I<acotanh>, have aliases I<rho>, I<theta>, I<ln>,
1646 I<cosec>, I<cotan>, I<acosec>, I<acotan>, I<cosech>, I<cotanh>,
1647 I<acosech>, I<acotanh>, respectively. C<Re>, C<Im>, C<arg>, C<abs>,
1648 C<rho>, and C<theta> can be used also also mutators. The C<cbrt>
1649 returns only one of the solutions: if you want all three, use the
1652 The I<root> function is available to compute all the I<n>
1653 roots of some complex, where I<n> is a strictly positive integer.
1654 There are exactly I<n> such roots, returned as a list. Getting the
1655 number mathematicians call C<j> such that:
1659 is a simple matter of writing:
1661 $j = ((root(1, 3))[1];
1663 The I<k>th root for C<z = [r,t]> is given by:
1665 (root(z, n))[k] = r**(1/n) * exp(i * (t + 2*k*pi)/n)
1667 The I<spaceship> comparison operator, E<lt>=E<gt>, is also defined. In
1668 order to ensure its restriction to real numbers is conform to what you
1669 would expect, the comparison is run on the real part of the complex
1670 number first, and imaginary parts are compared only when the real
1675 To create a complex number, use either:
1677 $z = Math::Complex->make(3, 4);
1680 if you know the cartesian form of the number, or
1684 if you like. To create a number using the polar form, use either:
1686 $z = Math::Complex->emake(5, pi/3);
1687 $x = cplxe(5, pi/3);
1689 instead. The first argument is the modulus, the second is the angle
1690 (in radians, the full circle is 2*pi). (Mnemonic: C<e> is used as a
1691 notation for complex numbers in the polar form).
1693 It is possible to write:
1695 $x = cplxe(-3, pi/4);
1697 but that will be silently converted into C<[3,-3pi/4]>, since the
1698 modulus must be non-negative (it represents the distance to the origin
1699 in the complex plane).
1701 It is also possible to have a complex number as either argument of
1702 either the C<make> or C<emake>: the appropriate component of
1703 the argument will be used.
1708 =head1 STRINGIFICATION
1710 When printed, a complex number is usually shown under its cartesian
1711 style I<a+bi>, but there are legitimate cases where the polar style
1712 I<[r,t]> is more appropriate.
1714 By calling the class method C<Math::Complex::display_format> and
1715 supplying either C<"polar"> or C<"cartesian"> as an argument, you
1716 override the default display style, which is C<"cartesian">. Not
1717 supplying any argument returns the current settings.
1719 This default can be overridden on a per-number basis by calling the
1720 C<display_format> method instead. As before, not supplying any argument
1721 returns the current display style for this number. Otherwise whatever you
1722 specify will be the new display style for I<this> particular number.
1728 Math::Complex::display_format('polar');
1729 $j = (root(1, 3))[1];
1730 print "j = $j\n"; # Prints "j = [1,2pi/3]"
1731 $j->display_format('cartesian');
1732 print "j = $j\n"; # Prints "j = -0.5+0.866025403784439i"
1734 The polar style attempts to emphasize arguments like I<k*pi/n>
1735 (where I<n> is a positive integer and I<k> an integer within [-9, +9]),
1736 this is called I<polar pretty-printing>.
1738 =head2 CHANGED IN PERL 5.6
1740 The C<display_format> class method and the corresponding
1741 C<display_format> object method can now be called using
1742 a parameter hash instead of just a one parameter.
1744 The old display format style, which can have values C<"cartesian"> or
1745 C<"polar">, can be changed using the C<"style"> parameter. (The one
1746 parameter calling convention also still works.)
1748 There are two new display parameters.
1750 The first one is C<"format">, which is a sprintf()-style format
1751 string to be used for both parts of the complex number(s). The
1752 default is C<undef>, which corresponds usually (this is somewhat
1753 system-dependent) to C<"%.15g">. You can revert to the default by
1754 setting the format string to C<undef>.
1756 # the $j from the above example
1758 $j->display_format('format' => '%.5f');
1759 print "j = $j\n"; # Prints "j = -0.50000+0.86603i"
1760 $j->display_format('format' => '%.6f');
1761 print "j = $j\n"; # Prints "j = -0.5+0.86603i"
1763 Notice that this affects also the return values of the
1764 C<display_format> methods: in list context the whole parameter hash
1765 will be returned, as opposed to only the style parameter value. If
1766 you want to know the whole truth for a complex number, you must call
1767 both the class method and the object method:
1769 The second new display parameter is C<"polar_pretty_print">, which can
1770 be set to true or false, the default being true. See the previous
1771 section for what this means.
1775 Thanks to overloading, the handling of arithmetics with complex numbers
1776 is simple and almost transparent.
1778 Here are some examples:
1782 $j = cplxe(1, 2*pi/3); # $j ** 3 == 1
1783 print "j = $j, j**3 = ", $j ** 3, "\n";
1784 print "1 + j + j**2 = ", 1 + $j + $j**2, "\n";
1786 $z = -16 + 0*i; # Force it to be a complex
1787 print "sqrt($z) = ", sqrt($z), "\n";
1789 $k = exp(i * 2*pi/3);
1790 print "$j - $k = ", $j - $k, "\n";
1792 $z->Re(3); # Re, Im, arg, abs,
1793 $j->arg(2); # (the last two aka rho, theta)
1794 # can be used also as mutators.
1796 =head1 ERRORS DUE TO DIVISION BY ZERO OR LOGARITHM OF ZERO
1798 The division (/) and the following functions
1804 atanh asech acsch acoth
1806 cannot be computed for all arguments because that would mean dividing
1807 by zero or taking logarithm of zero. These situations cause fatal
1808 runtime errors looking like this
1810 cot(0): Division by zero.
1811 (Because in the definition of cot(0), the divisor sin(0) is 0)
1816 atanh(-1): Logarithm of zero.
1819 For the C<csc>, C<cot>, C<asec>, C<acsc>, C<acot>, C<csch>, C<coth>,
1820 C<asech>, C<acsch>, the argument cannot be C<0> (zero). For the the
1821 logarithmic functions and the C<atanh>, C<acoth>, the argument cannot
1822 be C<1> (one). For the C<atanh>, C<acoth>, the argument cannot be
1823 C<-1> (minus one). For the C<atan>, C<acot>, the argument cannot be
1824 C<i> (the imaginary unit). For the C<atan>, C<acoth>, the argument
1825 cannot be C<-i> (the negative imaginary unit). For the C<tan>,
1826 C<sec>, C<tanh>, the argument cannot be I<pi/2 + k * pi>, where I<k>
1829 Note that because we are operating on approximations of real numbers,
1830 these errors can happen when merely `too close' to the singularities
1831 listed above. For example C<tan(2*atan2(1,1)+1e-15)> will die of
1834 =head1 ERRORS DUE TO INDIGESTIBLE ARGUMENTS
1836 The C<make> and C<emake> accept both real and complex arguments.
1837 When they cannot recognize the arguments they will die with error
1838 messages like the following
1840 Math::Complex::make: Cannot take real part of ...
1841 Math::Complex::make: Cannot take real part of ...
1842 Math::Complex::emake: Cannot take rho of ...
1843 Math::Complex::emake: Cannot take theta of ...
1847 Saying C<use Math::Complex;> exports many mathematical routines in the
1848 caller environment and even overrides some (C<sqrt>, C<log>).
1849 This is construed as a feature by the Authors, actually... ;-)
1851 All routines expect to be given real or complex numbers. Don't attempt to
1852 use BigFloat, since Perl has currently no rule to disambiguate a '+'
1853 operation (for instance) between two overloaded entities.
1855 In Cray UNICOS there is some strange numerical instability that results
1856 in root(), cos(), sin(), cosh(), sinh(), losing accuracy fast. Beware.
1857 The bug may be in UNICOS math libs, in UNICOS C compiler, in Math::Complex.
1858 Whatever it is, it does not manifest itself anywhere else where Perl runs.
1862 Raphael Manfredi <F<Raphael_Manfredi@pobox.com>> and
1863 Jarkko Hietaniemi <F<jhi@iki.fi>>.
1865 Extensive patches by Daniel S. Lewart <F<d-lewart@uiuc.edu>>.