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
14 our($VERSION, @ISA, @EXPORT, %EXPORT_TAGS);
18 $VERSION = sprintf("%s", q$Id: Complex.pm,v 1.26 1998/11/01 00:00:00 dsl Exp $ =~ /(\d+\.\d+)/);
25 csc cosec sec cot cotan
27 acsc acosec asec acot acotan
29 csch cosech sech coth cotanh
31 acsch acosech asech acoth acotanh
70 my $package = 'Math::Complex'; # Package name
71 my %DISPLAY_FORMAT = ('style' => 'cartesian',
72 'polar_pretty_print' => 1);
73 my $eps = 1e-14; # Epsilon
75 my $Inf = CORE::exp(CORE::exp(30));
76 $Inf = "Inf" if !defined $Inf || !$Inf > 0;
79 # Object attributes (internal):
80 # cartesian [real, imaginary] -- cartesian form
81 # polar [rho, theta] -- polar form
82 # c_dirty cartesian form not up-to-date
83 # p_dirty polar form not up-to-date
84 # display display format (package's global when not set)
87 # Die on bad *make() arguments.
90 die "@{[(caller(1))[3]]}: Cannot take $_[0] of $_[1].\n";
96 # Create a new complex number (cartesian form)
99 my $self = bless {}, shift;
103 if ( $rre eq ref $self ) {
106 _cannot_make("real part", $rre);
111 if ( $rim eq ref $self ) {
114 _cannot_make("imaginary part", $rim);
117 $self->{'cartesian'} = [ $re, $im ];
118 $self->{c_dirty} = 0;
119 $self->{p_dirty} = 1;
120 $self->display_format('cartesian');
127 # Create a new complex number (exponential form)
130 my $self = bless {}, shift;
131 my ($rho, $theta) = @_;
134 if ( $rrh eq ref $self ) {
137 _cannot_make("rho", $rrh);
140 my $rth = ref $theta;
142 if ( $rth eq ref $self ) {
143 $theta = theta($theta);
145 _cannot_make("theta", $rth);
150 $theta = ($theta <= 0) ? $theta + pi() : $theta - pi();
152 $self->{'polar'} = [$rho, $theta];
153 $self->{p_dirty} = 0;
154 $self->{c_dirty} = 1;
155 $self->display_format('polar');
159 sub new { &make } # For backward compatibility only.
164 # Creates a complex number from a (re, im) tuple.
165 # This avoids the burden of writing Math::Complex->make(re, im).
169 return __PACKAGE__->make($re, defined $im ? $im : 0);
175 # Creates a complex number from a (rho, theta) tuple.
176 # This avoids the burden of writing Math::Complex->emake(rho, theta).
179 my ($rho, $theta) = @_;
180 return __PACKAGE__->emake($rho, defined $theta ? $theta : 0);
186 # The number defined as pi = 180 degrees
188 sub pi () { 4 * CORE::atan2(1, 1) }
195 sub pit2 () { 2 * pi }
202 sub pip2 () { pi / 2 }
207 # One degree in radians, used in stringify_polar.
210 sub deg1 () { pi / 180 }
217 sub uplog10 () { 1 / CORE::log(10) }
222 # The number defined as i*i = -1;
227 $i->{'cartesian'} = [0, 1];
228 $i->{'polar'} = [1, pip2];
242 # Attribute access/set routines
245 sub cartesian {$_[0]->{c_dirty} ?
246 $_[0]->update_cartesian : $_[0]->{'cartesian'}}
247 sub polar {$_[0]->{p_dirty} ?
248 $_[0]->update_polar : $_[0]->{'polar'}}
250 sub set_cartesian { $_[0]->{p_dirty}++; $_[0]->{'cartesian'} = $_[1] }
251 sub set_polar { $_[0]->{c_dirty}++; $_[0]->{'polar'} = $_[1] }
256 # Recompute and return the cartesian form, given accurate polar form.
258 sub update_cartesian {
260 my ($r, $t) = @{$self->{'polar'}};
261 $self->{c_dirty} = 0;
262 return $self->{'cartesian'} = [$r * CORE::cos($t), $r * CORE::sin($t)];
269 # Recompute and return the polar form, given accurate cartesian form.
273 my ($x, $y) = @{$self->{'cartesian'}};
274 $self->{p_dirty} = 0;
275 return $self->{'polar'} = [0, 0] if $x == 0 && $y == 0;
276 return $self->{'polar'} = [CORE::sqrt($x*$x + $y*$y),
277 CORE::atan2($y, $x)];
286 my ($z1, $z2, $regular) = @_;
287 my ($re1, $im1) = @{$z1->cartesian};
288 $z2 = cplx($z2) unless ref $z2;
289 my ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
290 unless (defined $regular) {
291 $z1->set_cartesian([$re1 + $re2, $im1 + $im2]);
294 return (ref $z1)->make($re1 + $re2, $im1 + $im2);
303 my ($z1, $z2, $inverted) = @_;
304 my ($re1, $im1) = @{$z1->cartesian};
305 $z2 = cplx($z2) unless ref $z2;
306 my ($re2, $im2) = @{$z2->cartesian};
307 unless (defined $inverted) {
308 $z1->set_cartesian([$re1 - $re2, $im1 - $im2]);
312 (ref $z1)->make($re2 - $re1, $im2 - $im1) :
313 (ref $z1)->make($re1 - $re2, $im1 - $im2);
323 my ($z1, $z2, $regular) = @_;
324 if ($z1->{p_dirty} == 0 and ref $z2 and $z2->{p_dirty} == 0) {
325 # if both polar better use polar to avoid rounding errors
326 my ($r1, $t1) = @{$z1->polar};
327 my ($r2, $t2) = @{$z2->polar};
329 if ($t > pi()) { $t -= pit2 }
330 elsif ($t <= -pi()) { $t += pit2 }
331 unless (defined $regular) {
332 $z1->set_polar([$r1 * $r2, $t]);
335 return (ref $z1)->emake($r1 * $r2, $t);
337 my ($x1, $y1) = @{$z1->cartesian};
339 my ($x2, $y2) = @{$z2->cartesian};
340 return (ref $z1)->make($x1*$x2-$y1*$y2, $x1*$y2+$y1*$x2);
342 return (ref $z1)->make($x1*$z2, $y1*$z2);
350 # Die on division by zero.
353 my $mess = "$_[0]: Division by zero.\n";
356 $mess .= "(Because in the definition of $_[0], the divisor ";
357 $mess .= "$_[1] " unless ("$_[1]" eq '0');
363 $mess .= "Died at $up[1] line $up[2].\n";
374 my ($z1, $z2, $inverted) = @_;
375 if ($z1->{p_dirty} == 0 and ref $z2 and $z2->{p_dirty} == 0) {
376 # if both polar better use polar to avoid rounding errors
377 my ($r1, $t1) = @{$z1->polar};
378 my ($r2, $t2) = @{$z2->polar};
381 _divbyzero "$z2/0" if ($r1 == 0);
383 if ($t > pi()) { $t -= pit2 }
384 elsif ($t <= -pi()) { $t += pit2 }
385 return (ref $z1)->emake($r2 / $r1, $t);
387 _divbyzero "$z1/0" if ($r2 == 0);
389 if ($t > pi()) { $t -= pit2 }
390 elsif ($t <= -pi()) { $t += pit2 }
391 return (ref $z1)->emake($r1 / $r2, $t);
396 ($x2, $y2) = @{$z1->cartesian};
397 $d = $x2*$x2 + $y2*$y2;
398 _divbyzero "$z2/0" if $d == 0;
399 return (ref $z1)->make(($x2*$z2)/$d, -($y2*$z2)/$d);
401 my ($x1, $y1) = @{$z1->cartesian};
403 ($x2, $y2) = @{$z2->cartesian};
404 $d = $x2*$x2 + $y2*$y2;
405 _divbyzero "$z1/0" if $d == 0;
406 my $u = ($x1*$x2 + $y1*$y2)/$d;
407 my $v = ($y1*$x2 - $x1*$y2)/$d;
408 return (ref $z1)->make($u, $v);
410 _divbyzero "$z1/0" if $z2 == 0;
411 return (ref $z1)->make($x1/$z2, $y1/$z2);
420 # Computes z1**z2 = exp(z2 * log z1)).
423 my ($z1, $z2, $inverted) = @_;
425 return 1 if $z1 == 0 || $z2 == 1;
426 return 0 if $z2 == 0 && Re($z1) > 0;
428 return 1 if $z2 == 0 || $z1 == 1;
429 return 0 if $z1 == 0 && Re($z2) > 0;
431 my $w = $inverted ? &exp($z1 * &log($z2))
432 : &exp($z2 * &log($z1));
433 # If both arguments cartesian, return cartesian, else polar.
434 return $z1->{c_dirty} == 0 &&
435 (not ref $z2 or $z2->{c_dirty} == 0) ?
436 cplx(@{$w->cartesian}) : $w;
442 # Computes z1 <=> z2.
443 # Sorts on the real part first, then on the imaginary part. Thus 2-4i < 3+8i.
446 my ($z1, $z2, $inverted) = @_;
447 my ($re1, $im1) = ref $z1 ? @{$z1->cartesian} : ($z1, 0);
448 my ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
449 my $sgn = $inverted ? -1 : 1;
450 return $sgn * ($re1 <=> $re2) if $re1 != $re2;
451 return $sgn * ($im1 <=> $im2);
459 # (Required in addition to spaceship() because of NaNs.)
461 my ($z1, $z2, $inverted) = @_;
462 my ($re1, $im1) = ref $z1 ? @{$z1->cartesian} : ($z1, 0);
463 my ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
464 return $re1 == $re2 && $im1 == $im2 ? 1 : 0;
475 my ($r, $t) = @{$z->polar};
476 $t = ($t <= 0) ? $t + pi : $t - pi;
477 return (ref $z)->emake($r, $t);
479 my ($re, $im) = @{$z->cartesian};
480 return (ref $z)->make(-$re, -$im);
486 # Compute complex's conjugate.
491 my ($r, $t) = @{$z->polar};
492 return (ref $z)->emake($r, -$t);
494 my ($re, $im) = @{$z->cartesian};
495 return (ref $z)->make($re, -$im);
501 # Compute or set complex's norm (rho).
509 return CORE::abs($z);
513 $z->{'polar'} = [ $rho, ${$z->polar}[1] ];
518 return ${$z->polar}[0];
525 if ($$theta > pi()) { $$theta -= pit2 }
526 elsif ($$theta <= -pi()) { $$theta += pit2 }
532 # Compute or set complex's argument (theta).
535 my ($z, $theta) = @_;
536 return $z unless ref $z;
537 if (defined $theta) {
539 $z->{'polar'} = [ ${$z->polar}[0], $theta ];
543 $theta = ${$z->polar}[1];
554 # It is quite tempting to use wantarray here so that in list context
555 # sqrt() would return the two solutions. This, however, would
558 # print "sqrt(z) = ", sqrt($z), "\n";
560 # The two values would be printed side by side without no intervening
561 # whitespace, quite confusing.
562 # Therefore if you want the two solutions use the root().
566 my ($re, $im) = ref $z ? @{$z->cartesian} : ($z, 0);
567 return $re < 0 ? cplx(0, CORE::sqrt(-$re)) : CORE::sqrt($re)
569 my ($r, $t) = @{$z->polar};
570 return (ref $z)->emake(CORE::sqrt($r), $t/2);
576 # Compute cbrt(z) (cubic root).
578 # Why are we not returning three values? The same answer as for sqrt().
583 -CORE::exp(CORE::log(-$z)/3) :
584 ($z > 0 ? CORE::exp(CORE::log($z)/3): 0)
586 my ($r, $t) = @{$z->polar};
588 return (ref $z)->emake(CORE::exp(CORE::log($r)/3), $t/3);
597 my $mess = "Root $_[0] illegal, root rank must be positive integer.\n";
601 $mess .= "Died at $up[1] line $up[2].\n";
609 # Computes all nth root for z, returning an array whose size is n.
610 # `n' must be a positive integer.
612 # The roots are given by (for k = 0..n-1):
614 # z^(1/n) = r^(1/n) (cos ((t+2 k pi)/n) + i sin ((t+2 k pi)/n))
618 _rootbad($n) if ($n < 1 or int($n) != $n);
619 my ($r, $t) = ref $z ?
620 @{$z->polar} : (CORE::abs($z), $z >= 0 ? 0 : pi);
623 my $theta_inc = pit2 / $n;
624 my $rho = $r ** (1/$n);
626 my $cartesian = ref $z && $z->{c_dirty} == 0;
627 for ($k = 0, $theta = $t / $n; $k < $n; $k++, $theta += $theta_inc) {
628 my $w = cplxe($rho, $theta);
629 # Yes, $cartesian is loop invariant.
630 push @root, $cartesian ? cplx(@{$w->cartesian}) : $w;
638 # Return or set Re(z).
642 return $z unless ref $z;
644 $z->{'cartesian'} = [ $Re, ${$z->cartesian}[1] ];
648 return ${$z->cartesian}[0];
655 # Return or set Im(z).
659 return $z unless ref $z;
661 $z->{'cartesian'} = [ ${$z->cartesian}[0], $Im ];
665 return ${$z->cartesian}[1];
672 # Return or set rho(w).
675 Math::Complex::abs(@_);
681 # Return or set theta(w).
684 Math::Complex::arg(@_);
694 my ($x, $y) = @{$z->cartesian};
695 return (ref $z)->emake(CORE::exp($x), $y);
701 # Die on logarithm of zero.
704 my $mess = "$_[0]: Logarithm of zero.\n";
707 $mess .= "(Because in the definition of $_[0], the argument ";
708 $mess .= "$_[1] " unless ($_[1] eq '0');
714 $mess .= "Died at $up[1] line $up[2].\n";
727 _logofzero("log") if $z == 0;
728 return $z > 0 ? CORE::log($z) : cplx(CORE::log(-$z), pi);
730 my ($r, $t) = @{$z->polar};
731 _logofzero("log") if $r == 0;
732 if ($t > pi()) { $t -= pit2 }
733 elsif ($t <= -pi()) { $t += pit2 }
734 return (ref $z)->make(CORE::log($r), $t);
742 sub ln { Math::Complex::log(@_) }
751 return Math::Complex::log($_[0]) * uplog10;
757 # Compute logn(z,n) = log(z) / log(n)
761 $z = cplx($z, 0) unless ref $z;
762 my $logn = $logn{$n};
763 $logn = $logn{$n} = CORE::log($n) unless defined $logn; # Cache log(n)
764 return &log($z) / $logn;
770 # Compute cos(z) = (exp(iz) + exp(-iz))/2.
774 return CORE::cos($z) unless ref $z;
775 my ($x, $y) = @{$z->cartesian};
776 my $ey = CORE::exp($y);
777 my $sx = CORE::sin($x);
778 my $cx = CORE::cos($x);
779 my $ey_1 = $ey ? 1 / $ey : $Inf;
780 return (ref $z)->make($cx * ($ey + $ey_1)/2,
781 $sx * ($ey_1 - $ey)/2);
787 # Compute sin(z) = (exp(iz) - exp(-iz))/2.
791 return CORE::sin($z) unless ref $z;
792 my ($x, $y) = @{$z->cartesian};
793 my $ey = CORE::exp($y);
794 my $sx = CORE::sin($x);
795 my $cx = CORE::cos($x);
796 my $ey_1 = $ey ? 1 / $ey : $Inf;
797 return (ref $z)->make($sx * ($ey + $ey_1)/2,
798 $cx * ($ey - $ey_1)/2);
804 # Compute tan(z) = sin(z) / cos(z).
809 _divbyzero "tan($z)", "cos($z)" if $cz == 0;
810 return &sin($z) / $cz;
816 # Computes the secant sec(z) = 1 / cos(z).
821 _divbyzero "sec($z)", "cos($z)" if ($cz == 0);
828 # Computes the cosecant csc(z) = 1 / sin(z).
833 _divbyzero "csc($z)", "sin($z)" if ($sz == 0);
842 sub cosec { Math::Complex::csc(@_) }
847 # Computes cot(z) = cos(z) / sin(z).
852 _divbyzero "cot($z)", "sin($z)" if ($sz == 0);
853 return &cos($z) / $sz;
861 sub cotan { Math::Complex::cot(@_) }
866 # Computes the arc cosine acos(z) = -i log(z + sqrt(z*z-1)).
870 return CORE::atan2(CORE::sqrt(1-$z*$z), $z)
871 if (! ref $z) && CORE::abs($z) <= 1;
872 my ($x, $y) = ref $z ? @{$z->cartesian} : ($z, 0);
873 return 0 if $x == 1 && $y == 0;
874 my $t1 = CORE::sqrt(($x+1)*($x+1) + $y*$y);
875 my $t2 = CORE::sqrt(($x-1)*($x-1) + $y*$y);
876 my $alpha = ($t1 + $t2)/2;
877 my $beta = ($t1 - $t2)/2;
878 $alpha = 1 if $alpha < 1;
879 if ($beta > 1) { $beta = 1 }
880 elsif ($beta < -1) { $beta = -1 }
881 my $u = CORE::atan2(CORE::sqrt(1-$beta*$beta), $beta);
882 my $v = CORE::log($alpha + CORE::sqrt($alpha*$alpha-1));
883 $v = -$v if $y > 0 || ($y == 0 && $x < -1);
884 return __PACKAGE__->make($u, $v);
890 # Computes the arc sine asin(z) = -i log(iz + sqrt(1-z*z)).
894 return CORE::atan2($z, CORE::sqrt(1-$z*$z))
895 if (! ref $z) && CORE::abs($z) <= 1;
896 my ($x, $y) = ref $z ? @{$z->cartesian} : ($z, 0);
897 return 0 if $x == 0 && $y == 0;
898 my $t1 = CORE::sqrt(($x+1)*($x+1) + $y*$y);
899 my $t2 = CORE::sqrt(($x-1)*($x-1) + $y*$y);
900 my $alpha = ($t1 + $t2)/2;
901 my $beta = ($t1 - $t2)/2;
902 $alpha = 1 if $alpha < 1;
903 if ($beta > 1) { $beta = 1 }
904 elsif ($beta < -1) { $beta = -1 }
905 my $u = CORE::atan2($beta, CORE::sqrt(1-$beta*$beta));
906 my $v = -CORE::log($alpha + CORE::sqrt($alpha*$alpha-1));
907 $v = -$v if $y > 0 || ($y == 0 && $x < -1);
908 return __PACKAGE__->make($u, $v);
914 # Computes the arc tangent atan(z) = i/2 log((i+z) / (i-z)).
918 return CORE::atan2($z, 1) unless ref $z;
919 my ($x, $y) = ref $z ? @{$z->cartesian} : ($z, 0);
920 return 0 if $x == 0 && $y == 0;
921 _divbyzero "atan(i)" if ( $z == i);
922 _logofzero "atan(-i)" if (-$z == i); # -i is a bad file test...
923 my $log = &log((i + $z) / (i - $z));
930 # Computes the arc secant asec(z) = acos(1 / z).
934 _divbyzero "asec($z)", $z if ($z == 0);
941 # Computes the arc cosecant acsc(z) = asin(1 / z).
945 _divbyzero "acsc($z)", $z if ($z == 0);
954 sub acosec { Math::Complex::acsc(@_) }
959 # Computes the arc cotangent acot(z) = atan(1 / z)
963 _divbyzero "acot(0)" if $z == 0;
964 return ($z >= 0) ? CORE::atan2(1, $z) : CORE::atan2(-1, -$z)
966 _divbyzero "acot(i)" if ($z - i == 0);
967 _logofzero "acot(-i)" if ($z + i == 0);
976 sub acotan { Math::Complex::acot(@_) }
981 # Computes the hyperbolic cosine cosh(z) = (exp(z) + exp(-z))/2.
988 return $ex ? ($ex + 1/$ex)/2 : $Inf;
990 my ($x, $y) = @{$z->cartesian};
991 my $cy = CORE::cos($y);
992 my $sy = CORE::cos($y);
994 my $ex_1 = $ex ? 1 / $ex : $Inf;
995 return (ref $z)->make(CORE::cos($y) * ($ex + $ex_1)/2,
996 CORE::sin($y) * ($ex - $ex_1)/2);
1002 # Computes the hyperbolic sine sinh(z) = (exp(z) - exp(-z))/2.
1008 return 0 if $z == 0;
1009 $ex = CORE::exp($z);
1010 return $ex ? ($ex - 1/$ex)/2 : "-$Inf";
1012 my ($x, $y) = @{$z->cartesian};
1013 my $cy = CORE::cos($y);
1014 my $sy = CORE::sin($y);
1015 $ex = CORE::exp($x);
1016 my $ex_1 = $ex ? 1 / $ex : $Inf;
1017 return (ref $z)->make($cy * ($ex - $ex_1)/2,
1018 $sy * ($ex + $ex_1)/2);
1024 # Computes the hyperbolic tangent tanh(z) = sinh(z) / cosh(z).
1029 _divbyzero "tanh($z)", "cosh($z)" if ($cz == 0);
1030 return sinh($z) / $cz;
1036 # Computes the hyperbolic secant sech(z) = 1 / cosh(z).
1041 _divbyzero "sech($z)", "cosh($z)" if ($cz == 0);
1048 # Computes the hyperbolic cosecant csch(z) = 1 / sinh(z).
1053 _divbyzero "csch($z)", "sinh($z)" if ($sz == 0);
1062 sub cosech { Math::Complex::csch(@_) }
1067 # Computes the hyperbolic cotangent coth(z) = cosh(z) / sinh(z).
1072 _divbyzero "coth($z)", "sinh($z)" if $sz == 0;
1073 return cosh($z) / $sz;
1081 sub cotanh { Math::Complex::coth(@_) }
1086 # Computes the arc hyperbolic cosine acosh(z) = log(z + sqrt(z*z-1)).
1093 my ($re, $im) = @{$z->cartesian};
1095 return CORE::log($re + CORE::sqrt($re*$re - 1))
1097 return cplx(0, CORE::atan2(CORE::sqrt(1 - $re*$re), $re))
1098 if CORE::abs($re) < 1;
1100 my $s = &sqrt($z*$z - 1);
1102 $t = 1/(2*$s) if $t == 0 || $t && &abs(cosh(&log($t)) - $z) > $eps;
1109 # Computes the arc hyperbolic sine asinh(z) = log(z + sqrt(z*z+1))
1114 my $t = $z + CORE::sqrt($z*$z + 1);
1115 return CORE::log($t) if $t;
1117 my $s = &sqrt($z*$z + 1);
1119 # Try Taylor series if looking bad.
1120 $t = 1/(2*$s) if $t == 0 || $t && &abs(sinh(&log($t)) - $z) > $eps;
1127 # Computes the arc hyperbolic tangent atanh(z) = 1/2 log((1+z) / (1-z)).
1132 return CORE::log((1 + $z)/(1 - $z))/2 if CORE::abs($z) < 1;
1135 _divbyzero 'atanh(1)', "1 - $z" if (1 - $z == 0);
1136 _logofzero 'atanh(-1)' if (1 + $z == 0);
1137 return 0.5 * &log((1 + $z) / (1 - $z));
1143 # Computes the hyperbolic arc secant asech(z) = acosh(1 / z).
1147 _divbyzero 'asech(0)', "$z" if ($z == 0);
1148 return acosh(1 / $z);
1154 # Computes the hyperbolic arc cosecant acsch(z) = asinh(1 / z).
1158 _divbyzero 'acsch(0)', $z if ($z == 0);
1159 return asinh(1 / $z);
1165 # Alias for acosh().
1167 sub acosech { Math::Complex::acsch(@_) }
1172 # Computes the arc hyperbolic cotangent acoth(z) = 1/2 log((1+z) / (z-1)).
1176 _divbyzero 'acoth(0)' if ($z == 0);
1178 return CORE::log(($z + 1)/($z - 1))/2 if CORE::abs($z) > 1;
1181 _divbyzero 'acoth(1)', "$z - 1" if ($z - 1 == 0);
1182 _logofzero 'acoth(-1)', "1 + $z" if (1 + $z == 0);
1183 return &log((1 + $z) / ($z - 1)) / 2;
1191 sub acotanh { Math::Complex::acoth(@_) }
1196 # Compute atan(z1/z2).
1199 my ($z1, $z2, $inverted) = @_;
1200 my ($re1, $im1, $re2, $im2);
1202 ($re1, $im1) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
1203 ($re2, $im2) = @{$z1->cartesian};
1205 ($re1, $im1) = @{$z1->cartesian};
1206 ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
1209 return CORE::atan2($re1, $re2) if $im1 == 0;
1210 return ($im1<=>0) * pip2 if $re2 == 0;
1212 my $w = atan($z1/$z2);
1213 my ($u, $v) = ref $w ? @{$w->cartesian} : ($w, 0);
1214 $u += pi if $re2 < 0;
1215 $u -= pit2 if $u > pi;
1216 return cplx($u, $v);
1223 # Set (get if no argument) the display format for all complex numbers that
1224 # don't happen to have overridden it via ->display_format
1226 # When called as an object method, this actually sets the display format for
1227 # the current object.
1229 # Valid object formats are 'c' and 'p' for cartesian and polar. The first
1230 # letter is used actually, so the type can be fully spelled out for clarity.
1232 sub display_format {
1234 my %display_format = %DISPLAY_FORMAT;
1236 if (ref $self) { # Called as an object method
1237 if (exists $self->{display_format}) {
1238 my %obj = %{$self->{display_format}};
1239 @display_format{keys %obj} = values %obj;
1242 $display_format{style} = shift;
1245 @display_format{keys %new} = values %new;
1247 } else { # Called as a class method
1249 $display_format{style} = $self;
1252 @display_format{keys %new} = values %new;
1257 if (defined $self) {
1258 $self->{display_format} = { %display_format };
1261 %{$self->{display_format}} :
1262 $self->{display_format}->{style};
1265 %DISPLAY_FORMAT = %display_format;
1269 $DISPLAY_FORMAT{style};
1275 # Show nicely formatted complex number under its cartesian or polar form,
1276 # depending on the current display format:
1278 # . If a specific display format has been recorded for this object, use it.
1279 # . Otherwise, use the generic current default for all complex numbers,
1280 # which is a package global variable.
1285 my $style = $z->display_format;
1287 $style = $DISPLAY_FORMAT{style} unless defined $style;
1289 return $z->stringify_polar if $style =~ /^p/i;
1290 return $z->stringify_cartesian;
1294 # ->stringify_cartesian
1296 # Stringify as a cartesian representation 'a+bi'.
1298 sub stringify_cartesian {
1300 my ($x, $y) = @{$z->cartesian};
1303 my %format = $z->display_format;
1304 my $format = $format{format};
1307 if ($x =~ /^NaN[QS]?$/i) {
1310 if ($x =~ /^-?$Inf$/oi) {
1313 $re = defined $format ? sprintf($format, $x) : $x;
1321 if ($y == 1) { $im = "" }
1322 elsif ($y == -1) { $im = "-" }
1323 elsif ($y =~ /^(NaN[QS]?)$/i) {
1326 if ($y =~ /^-?$Inf$/oi) {
1329 $im = defined $format ? sprintf($format, $y) : $y;
1342 } elsif ($y > 0 || $im =~ /^NaN[QS]?i$/i) {
1343 $str .= "+" if defined $re;
1346 } elsif (!defined $re) {
1357 # Stringify as a polar representation '[r,t]'.
1359 sub stringify_polar {
1361 my ($r, $t) = @{$z->polar};
1364 my %format = $z->display_format;
1365 my $format = $format{format};
1367 if ($t =~ /^NaN[QS]?$/i || $t =~ /^-?$Inf$/oi) {
1369 } elsif ($t == pi) {
1371 } elsif ($r == 0 || $t == 0) {
1372 $theta = defined $format ? sprintf($format, $t) : $t;
1375 return "[$r,$theta]" if defined $theta;
1378 # Try to identify pi/n and friends.
1381 $t -= int(CORE::abs($t) / pit2) * pit2;
1383 if ($format{polar_pretty_print}) {
1385 for $a (2, 3, 4, 6, 8, 12, 16, 24, 30, 32, 36, 48, 60, 64, 72) {
1387 if (int($b) == $b) {
1388 $b = $b < 0 ? "-" : "" if CORE::abs($b) == 1;
1389 $theta = "${b}pi/$a";
1395 if (defined $format) {
1396 $r = sprintf($format, $r);
1397 $theta = sprintf($format, $theta) unless defined $theta;
1399 $theta = $t unless defined $theta;
1402 return "[$r,$theta]";
1410 Math::Complex - complex numbers and associated mathematical functions
1416 $z = Math::Complex->make(5, 6);
1418 $j = cplxe(1, 2*pi/3);
1422 This package lets you create and manipulate complex numbers. By default,
1423 I<Perl> limits itself to real numbers, but an extra C<use> statement brings
1424 full complex support, along with a full set of mathematical functions
1425 typically associated with and/or extended to complex numbers.
1427 If you wonder what complex numbers are, they were invented to be able to solve
1428 the following equation:
1432 and by definition, the solution is noted I<i> (engineers use I<j> instead since
1433 I<i> usually denotes an intensity, but the name does not matter). The number
1434 I<i> is a pure I<imaginary> number.
1436 The arithmetics with pure imaginary numbers works just like you would expect
1437 it with real numbers... you just have to remember that
1443 5i + 7i = i * (5 + 7) = 12i
1444 4i - 3i = i * (4 - 3) = i
1449 Complex numbers are numbers that have both a real part and an imaginary
1450 part, and are usually noted:
1454 where C<a> is the I<real> part and C<b> is the I<imaginary> part. The
1455 arithmetic with complex numbers is straightforward. You have to
1456 keep track of the real and the imaginary parts, but otherwise the
1457 rules used for real numbers just apply:
1459 (4 + 3i) + (5 - 2i) = (4 + 5) + i(3 - 2) = 9 + i
1460 (2 + i) * (4 - i) = 2*4 + 4i -2i -i*i = 8 + 2i + 1 = 9 + 2i
1462 A graphical representation of complex numbers is possible in a plane
1463 (also called the I<complex plane>, but it's really a 2D plane).
1468 is the point whose coordinates are (a, b). Actually, it would
1469 be the vector originating from (0, 0) to (a, b). It follows that the addition
1470 of two complex numbers is a vectorial addition.
1472 Since there is a bijection between a point in the 2D plane and a complex
1473 number (i.e. the mapping is unique and reciprocal), a complex number
1474 can also be uniquely identified with polar coordinates:
1478 where C<rho> is the distance to the origin, and C<theta> the angle between
1479 the vector and the I<x> axis. There is a notation for this using the
1480 exponential form, which is:
1482 rho * exp(i * theta)
1484 where I<i> is the famous imaginary number introduced above. Conversion
1485 between this form and the cartesian form C<a + bi> is immediate:
1487 a = rho * cos(theta)
1488 b = rho * sin(theta)
1490 which is also expressed by this formula:
1492 z = rho * exp(i * theta) = rho * (cos theta + i * sin theta)
1494 In other words, it's the projection of the vector onto the I<x> and I<y>
1495 axes. Mathematicians call I<rho> the I<norm> or I<modulus> and I<theta>
1496 the I<argument> of the complex number. The I<norm> of C<z> will be
1499 The polar notation (also known as the trigonometric
1500 representation) is much more handy for performing multiplications and
1501 divisions of complex numbers, whilst the cartesian notation is better
1502 suited for additions and subtractions. Real numbers are on the I<x>
1503 axis, and therefore I<theta> is zero or I<pi>.
1505 All the common operations that can be performed on a real number have
1506 been defined to work on complex numbers as well, and are merely
1507 I<extensions> of the operations defined on real numbers. This means
1508 they keep their natural meaning when there is no imaginary part, provided
1509 the number is within their definition set.
1511 For instance, the C<sqrt> routine which computes the square root of
1512 its argument is only defined for non-negative real numbers and yields a
1513 non-negative real number (it is an application from B<R+> to B<R+>).
1514 If we allow it to return a complex number, then it can be extended to
1515 negative real numbers to become an application from B<R> to B<C> (the
1516 set of complex numbers):
1518 sqrt(x) = x >= 0 ? sqrt(x) : sqrt(-x)*i
1520 It can also be extended to be an application from B<C> to B<C>,
1521 whilst its restriction to B<R> behaves as defined above by using
1522 the following definition:
1524 sqrt(z = [r,t]) = sqrt(r) * exp(i * t/2)
1526 Indeed, a negative real number can be noted C<[x,pi]> (the modulus
1527 I<x> is always non-negative, so C<[x,pi]> is really C<-x>, a negative
1528 number) and the above definition states that
1530 sqrt([x,pi]) = sqrt(x) * exp(i*pi/2) = [sqrt(x),pi/2] = sqrt(x)*i
1532 which is exactly what we had defined for negative real numbers above.
1533 The C<sqrt> returns only one of the solutions: if you want the both,
1534 use the C<root> function.
1536 All the common mathematical functions defined on real numbers that
1537 are extended to complex numbers share that same property of working
1538 I<as usual> when the imaginary part is zero (otherwise, it would not
1539 be called an extension, would it?).
1541 A I<new> operation possible on a complex number that is
1542 the identity for real numbers is called the I<conjugate>, and is noted
1543 with an horizontal bar above the number, or C<~z> here.
1550 z * ~z = (a + bi) * (a - bi) = a*a + b*b
1552 We saw that the norm of C<z> was noted C<abs(z)> and was defined as the
1553 distance to the origin, also known as:
1555 rho = abs(z) = sqrt(a*a + b*b)
1559 z * ~z = abs(z) ** 2
1561 If z is a pure real number (i.e. C<b == 0>), then the above yields:
1565 which is true (C<abs> has the regular meaning for real number, i.e. stands
1566 for the absolute value). This example explains why the norm of C<z> is
1567 noted C<abs(z)>: it extends the C<abs> function to complex numbers, yet
1568 is the regular C<abs> we know when the complex number actually has no
1569 imaginary part... This justifies I<a posteriori> our use of the C<abs>
1570 notation for the norm.
1574 Given the following notations:
1576 z1 = a + bi = r1 * exp(i * t1)
1577 z2 = c + di = r2 * exp(i * t2)
1578 z = <any complex or real number>
1580 the following (overloaded) operations are supported on complex numbers:
1582 z1 + z2 = (a + c) + i(b + d)
1583 z1 - z2 = (a - c) + i(b - d)
1584 z1 * z2 = (r1 * r2) * exp(i * (t1 + t2))
1585 z1 / z2 = (r1 / r2) * exp(i * (t1 - t2))
1586 z1 ** z2 = exp(z2 * log z1)
1588 abs(z) = r1 = sqrt(a*a + b*b)
1589 sqrt(z) = sqrt(r1) * exp(i * t/2)
1590 exp(z) = exp(a) * exp(i * b)
1591 log(z) = log(r1) + i*t
1592 sin(z) = 1/2i (exp(i * z1) - exp(-i * z))
1593 cos(z) = 1/2 (exp(i * z1) + exp(-i * z))
1594 atan2(z1, z2) = atan(z1/z2)
1596 The following extra operations are supported on both real and complex
1604 cbrt(z) = z ** (1/3)
1605 log10(z) = log(z) / log(10)
1606 logn(z, n) = log(z) / log(n)
1608 tan(z) = sin(z) / cos(z)
1614 asin(z) = -i * log(i*z + sqrt(1-z*z))
1615 acos(z) = -i * log(z + i*sqrt(1-z*z))
1616 atan(z) = i/2 * log((i+z) / (i-z))
1618 acsc(z) = asin(1 / z)
1619 asec(z) = acos(1 / z)
1620 acot(z) = atan(1 / z) = -i/2 * log((i+z) / (z-i))
1622 sinh(z) = 1/2 (exp(z) - exp(-z))
1623 cosh(z) = 1/2 (exp(z) + exp(-z))
1624 tanh(z) = sinh(z) / cosh(z) = (exp(z) - exp(-z)) / (exp(z) + exp(-z))
1626 csch(z) = 1 / sinh(z)
1627 sech(z) = 1 / cosh(z)
1628 coth(z) = 1 / tanh(z)
1630 asinh(z) = log(z + sqrt(z*z+1))
1631 acosh(z) = log(z + sqrt(z*z-1))
1632 atanh(z) = 1/2 * log((1+z) / (1-z))
1634 acsch(z) = asinh(1 / z)
1635 asech(z) = acosh(1 / z)
1636 acoth(z) = atanh(1 / z) = 1/2 * log((1+z) / (z-1))
1638 I<arg>, I<abs>, I<log>, I<csc>, I<cot>, I<acsc>, I<acot>, I<csch>,
1639 I<coth>, I<acosech>, I<acotanh>, have aliases I<rho>, I<theta>, I<ln>,
1640 I<cosec>, I<cotan>, I<acosec>, I<acotan>, I<cosech>, I<cotanh>,
1641 I<acosech>, I<acotanh>, respectively. C<Re>, C<Im>, C<arg>, C<abs>,
1642 C<rho>, and C<theta> can be used also also mutators. The C<cbrt>
1643 returns only one of the solutions: if you want all three, use the
1646 The I<root> function is available to compute all the I<n>
1647 roots of some complex, where I<n> is a strictly positive integer.
1648 There are exactly I<n> such roots, returned as a list. Getting the
1649 number mathematicians call C<j> such that:
1653 is a simple matter of writing:
1655 $j = ((root(1, 3))[1];
1657 The I<k>th root for C<z = [r,t]> is given by:
1659 (root(z, n))[k] = r**(1/n) * exp(i * (t + 2*k*pi)/n)
1661 The I<spaceship> comparison operator, E<lt>=E<gt>, is also defined. In
1662 order to ensure its restriction to real numbers is conform to what you
1663 would expect, the comparison is run on the real part of the complex
1664 number first, and imaginary parts are compared only when the real
1669 To create a complex number, use either:
1671 $z = Math::Complex->make(3, 4);
1674 if you know the cartesian form of the number, or
1678 if you like. To create a number using the polar form, use either:
1680 $z = Math::Complex->emake(5, pi/3);
1681 $x = cplxe(5, pi/3);
1683 instead. The first argument is the modulus, the second is the angle
1684 (in radians, the full circle is 2*pi). (Mnemonic: C<e> is used as a
1685 notation for complex numbers in the polar form).
1687 It is possible to write:
1689 $x = cplxe(-3, pi/4);
1691 but that will be silently converted into C<[3,-3pi/4]>, since the
1692 modulus must be non-negative (it represents the distance to the origin
1693 in the complex plane).
1695 It is also possible to have a complex number as either argument of
1696 either the C<make> or C<emake>: the appropriate component of
1697 the argument will be used.
1702 =head1 STRINGIFICATION
1704 When printed, a complex number is usually shown under its cartesian
1705 style I<a+bi>, but there are legitimate cases where the polar style
1706 I<[r,t]> is more appropriate.
1708 By calling the class method C<Math::Complex::display_format> and
1709 supplying either C<"polar"> or C<"cartesian"> as an argument, you
1710 override the default display style, which is C<"cartesian">. Not
1711 supplying any argument returns the current settings.
1713 This default can be overridden on a per-number basis by calling the
1714 C<display_format> method instead. As before, not supplying any argument
1715 returns the current display style for this number. Otherwise whatever you
1716 specify will be the new display style for I<this> particular number.
1722 Math::Complex::display_format('polar');
1723 $j = (root(1, 3))[1];
1724 print "j = $j\n"; # Prints "j = [1,2pi/3]"
1725 $j->display_format('cartesian');
1726 print "j = $j\n"; # Prints "j = -0.5+0.866025403784439i"
1728 The polar style attempts to emphasize arguments like I<k*pi/n>
1729 (where I<n> is a positive integer and I<k> an integer within [-9,+9]),
1730 this is called I<polar pretty-printing>.
1732 =head2 CHANGED IN PERL 5.6
1734 The C<display_format> class method and the corresponding
1735 C<display_format> object method can now be called using
1736 a parameter hash instead of just a one parameter.
1738 The old display format style, which can have values C<"cartesian"> or
1739 C<"polar">, can be changed using the C<"style"> parameter. (The one
1740 parameter calling convention also still works.)
1742 There are two new display parameters.
1744 The first one is C<"format">, which is a sprintf()-style format
1745 string to be used for both parts of the complex number(s). The
1746 default is C<undef>, which corresponds usually (this is somewhat
1747 system-dependent) to C<"%.15g">. You can revert to the default by
1748 setting the format string to C<undef>.
1750 # the $j from the above example
1752 $j->display_format('format' => '%.5f');
1753 print "j = $j\n"; # Prints "j = -0.50000+0.86603i"
1754 $j->display_format('format' => '%.6f');
1755 print "j = $j\n"; # Prints "j = -0.5+0.86603i"
1757 Notice that this affects also the return values of the
1758 C<display_format> methods: in list context the whole parameter hash
1759 will be returned, as opposed to only the style parameter value. If
1760 you want to know the whole truth for a complex number, you must call
1761 both the class method and the object method:
1763 The second new display parameter is C<"polar_pretty_print">, which can
1764 be set to true or false, the default being true. See the previous
1765 section for what this means.
1769 Thanks to overloading, the handling of arithmetics with complex numbers
1770 is simple and almost transparent.
1772 Here are some examples:
1776 $j = cplxe(1, 2*pi/3); # $j ** 3 == 1
1777 print "j = $j, j**3 = ", $j ** 3, "\n";
1778 print "1 + j + j**2 = ", 1 + $j + $j**2, "\n";
1780 $z = -16 + 0*i; # Force it to be a complex
1781 print "sqrt($z) = ", sqrt($z), "\n";
1783 $k = exp(i * 2*pi/3);
1784 print "$j - $k = ", $j - $k, "\n";
1786 $z->Re(3); # Re, Im, arg, abs,
1787 $j->arg(2); # (the last two aka rho, theta)
1788 # can be used also as mutators.
1790 =head1 ERRORS DUE TO DIVISION BY ZERO OR LOGARITHM OF ZERO
1792 The division (/) and the following functions
1798 atanh asech acsch acoth
1800 cannot be computed for all arguments because that would mean dividing
1801 by zero or taking logarithm of zero. These situations cause fatal
1802 runtime errors looking like this
1804 cot(0): Division by zero.
1805 (Because in the definition of cot(0), the divisor sin(0) is 0)
1810 atanh(-1): Logarithm of zero.
1813 For the C<csc>, C<cot>, C<asec>, C<acsc>, C<acot>, C<csch>, C<coth>,
1814 C<asech>, C<acsch>, the argument cannot be C<0> (zero). For the the
1815 logarithmic functions and the C<atanh>, C<acoth>, the argument cannot
1816 be C<1> (one). For the C<atanh>, C<acoth>, the argument cannot be
1817 C<-1> (minus one). For the C<atan>, C<acot>, the argument cannot be
1818 C<i> (the imaginary unit). For the C<atan>, C<acoth>, the argument
1819 cannot be C<-i> (the negative imaginary unit). For the C<tan>,
1820 C<sec>, C<tanh>, the argument cannot be I<pi/2 + k * pi>, where I<k>
1823 Note that because we are operating on approximations of real numbers,
1824 these errors can happen when merely `too close' to the singularities
1825 listed above. For example C<tan(2*atan2(1,1)+1e-15)> will die of
1828 =head1 ERRORS DUE TO INDIGESTIBLE ARGUMENTS
1830 The C<make> and C<emake> accept both real and complex arguments.
1831 When they cannot recognize the arguments they will die with error
1832 messages like the following
1834 Math::Complex::make: Cannot take real part of ...
1835 Math::Complex::make: Cannot take real part of ...
1836 Math::Complex::emake: Cannot take rho of ...
1837 Math::Complex::emake: Cannot take theta of ...
1841 Saying C<use Math::Complex;> exports many mathematical routines in the
1842 caller environment and even overrides some (C<sqrt>, C<log>).
1843 This is construed as a feature by the Authors, actually... ;-)
1845 All routines expect to be given real or complex numbers. Don't attempt to
1846 use BigFloat, since Perl has currently no rule to disambiguate a '+'
1847 operation (for instance) between two overloaded entities.
1849 In Cray UNICOS there is some strange numerical instability that results
1850 in root(), cos(), sin(), cosh(), sinh(), losing accuracy fast. Beware.
1851 The bug may be in UNICOS math libs, in UNICOS C compiler, in Math::Complex.
1852 Whatever it is, it does not manifest itself anywhere else where Perl runs.
1856 Raphael Manfredi <F<Raphael_Manfredi@pobox.com>> and
1857 Jarkko Hietaniemi <F<jhi@iki.fi>>.
1859 Extensive patches by Daniel S. Lewart <F<d-lewart@uiuc.edu>>.