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);
15 eval { require POSIX; import POSIX 'HUGE_VAL' };
16 if (exists &HUGE_VAL) {
17 $Inf = sprintf "%g", &HUGE_VAL;
20 $Inf = CORE::exp(CORE::exp(30));
21 $! = $e; # Clear ERANGE.
22 undef $Inf unless $Inf =~ /^inf$/; # Inf INF inf
24 $Inf = "Inf" if !defined $Inf || !($Inf > 0);
39 csc cosec sec cot cotan
41 acsc acosec asec acot acotan
43 csch cosech sech coth cotanh
45 acsch acosech asech acoth acotanh
84 my %DISPLAY_FORMAT = ('style' => 'cartesian',
85 'polar_pretty_print' => 1);
86 my $eps = 1e-14; # Epsilon
89 # Object attributes (internal):
90 # cartesian [real, imaginary] -- cartesian form
91 # polar [rho, theta] -- polar form
92 # c_dirty cartesian form not up-to-date
93 # p_dirty polar form not up-to-date
94 # display display format (package's global when not set)
97 # Die on bad *make() arguments.
100 die "@{[(caller(1))[3]]}: Cannot take $_[0] of $_[1].\n";
106 # Create a new complex number (cartesian form)
109 my $self = bless {}, shift;
113 if ( $rre eq ref $self ) {
116 _cannot_make("real part", $rre);
121 if ( $rim eq ref $self ) {
124 _cannot_make("imaginary part", $rim);
127 $self->{'cartesian'} = [ $re, $im ];
128 $self->{c_dirty} = 0;
129 $self->{p_dirty} = 1;
130 $self->display_format('cartesian');
137 # Create a new complex number (exponential form)
140 my $self = bless {}, shift;
141 my ($rho, $theta) = @_;
144 if ( $rrh eq ref $self ) {
147 _cannot_make("rho", $rrh);
150 my $rth = ref $theta;
152 if ( $rth eq ref $self ) {
153 $theta = theta($theta);
155 _cannot_make("theta", $rth);
160 $theta = ($theta <= 0) ? $theta + pi() : $theta - pi();
162 $self->{'polar'} = [$rho, $theta];
163 $self->{p_dirty} = 0;
164 $self->{c_dirty} = 1;
165 $self->display_format('polar');
169 sub new { &make } # For backward compatibility only.
174 # Creates a complex number from a (re, im) tuple.
175 # This avoids the burden of writing Math::Complex->make(re, im).
179 return __PACKAGE__->make($re, defined $im ? $im : 0);
185 # Creates a complex number from a (rho, theta) tuple.
186 # This avoids the burden of writing Math::Complex->emake(rho, theta).
189 my ($rho, $theta) = @_;
190 return __PACKAGE__->emake($rho, defined $theta ? $theta : 0);
196 # The number defined as pi = 180 degrees
198 sub pi () { 4 * CORE::atan2(1, 1) }
205 sub pit2 () { 2 * pi }
212 sub pip2 () { pi / 2 }
217 # One degree in radians, used in stringify_polar.
220 sub deg1 () { pi / 180 }
227 sub uplog10 () { 1 / CORE::log(10) }
232 # The number defined as i*i = -1;
237 $i->{'cartesian'} = [0, 1];
238 $i->{'polar'} = [1, pip2];
252 # Attribute access/set routines
255 sub cartesian {$_[0]->{c_dirty} ?
256 $_[0]->update_cartesian : $_[0]->{'cartesian'}}
257 sub polar {$_[0]->{p_dirty} ?
258 $_[0]->update_polar : $_[0]->{'polar'}}
260 sub set_cartesian { $_[0]->{p_dirty}++; $_[0]->{'cartesian'} = $_[1] }
261 sub set_polar { $_[0]->{c_dirty}++; $_[0]->{'polar'} = $_[1] }
266 # Recompute and return the cartesian form, given accurate polar form.
268 sub update_cartesian {
270 my ($r, $t) = @{$self->{'polar'}};
271 $self->{c_dirty} = 0;
272 return $self->{'cartesian'} = [$r * CORE::cos($t), $r * CORE::sin($t)];
279 # Recompute and return the polar form, given accurate cartesian form.
283 my ($x, $y) = @{$self->{'cartesian'}};
284 $self->{p_dirty} = 0;
285 return $self->{'polar'} = [0, 0] if $x == 0 && $y == 0;
286 return $self->{'polar'} = [CORE::sqrt($x*$x + $y*$y),
287 CORE::atan2($y, $x)];
296 my ($z1, $z2, $regular) = @_;
297 my ($re1, $im1) = @{$z1->cartesian};
298 $z2 = cplx($z2) unless ref $z2;
299 my ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
300 unless (defined $regular) {
301 $z1->set_cartesian([$re1 + $re2, $im1 + $im2]);
304 return (ref $z1)->make($re1 + $re2, $im1 + $im2);
313 my ($z1, $z2, $inverted) = @_;
314 my ($re1, $im1) = @{$z1->cartesian};
315 $z2 = cplx($z2) unless ref $z2;
316 my ($re2, $im2) = @{$z2->cartesian};
317 unless (defined $inverted) {
318 $z1->set_cartesian([$re1 - $re2, $im1 - $im2]);
322 (ref $z1)->make($re2 - $re1, $im2 - $im1) :
323 (ref $z1)->make($re1 - $re2, $im1 - $im2);
333 my ($z1, $z2, $regular) = @_;
334 if ($z1->{p_dirty} == 0 and ref $z2 and $z2->{p_dirty} == 0) {
335 # if both polar better use polar to avoid rounding errors
336 my ($r1, $t1) = @{$z1->polar};
337 my ($r2, $t2) = @{$z2->polar};
339 if ($t > pi()) { $t -= pit2 }
340 elsif ($t <= -pi()) { $t += pit2 }
341 unless (defined $regular) {
342 $z1->set_polar([$r1 * $r2, $t]);
345 return (ref $z1)->emake($r1 * $r2, $t);
347 my ($x1, $y1) = @{$z1->cartesian};
349 my ($x2, $y2) = @{$z2->cartesian};
350 return (ref $z1)->make($x1*$x2-$y1*$y2, $x1*$y2+$y1*$x2);
352 return (ref $z1)->make($x1*$z2, $y1*$z2);
360 # Die on division by zero.
363 my $mess = "$_[0]: Division by zero.\n";
366 $mess .= "(Because in the definition of $_[0], the divisor ";
367 $mess .= "$_[1] " unless ("$_[1]" eq '0');
373 $mess .= "Died at $up[1] line $up[2].\n";
384 my ($z1, $z2, $inverted) = @_;
385 if ($z1->{p_dirty} == 0 and ref $z2 and $z2->{p_dirty} == 0) {
386 # if both polar better use polar to avoid rounding errors
387 my ($r1, $t1) = @{$z1->polar};
388 my ($r2, $t2) = @{$z2->polar};
391 _divbyzero "$z2/0" if ($r1 == 0);
393 if ($t > pi()) { $t -= pit2 }
394 elsif ($t <= -pi()) { $t += pit2 }
395 return (ref $z1)->emake($r2 / $r1, $t);
397 _divbyzero "$z1/0" if ($r2 == 0);
399 if ($t > pi()) { $t -= pit2 }
400 elsif ($t <= -pi()) { $t += pit2 }
401 return (ref $z1)->emake($r1 / $r2, $t);
406 ($x2, $y2) = @{$z1->cartesian};
407 $d = $x2*$x2 + $y2*$y2;
408 _divbyzero "$z2/0" if $d == 0;
409 return (ref $z1)->make(($x2*$z2)/$d, -($y2*$z2)/$d);
411 my ($x1, $y1) = @{$z1->cartesian};
413 ($x2, $y2) = @{$z2->cartesian};
414 $d = $x2*$x2 + $y2*$y2;
415 _divbyzero "$z1/0" if $d == 0;
416 my $u = ($x1*$x2 + $y1*$y2)/$d;
417 my $v = ($y1*$x2 - $x1*$y2)/$d;
418 return (ref $z1)->make($u, $v);
420 _divbyzero "$z1/0" if $z2 == 0;
421 return (ref $z1)->make($x1/$z2, $y1/$z2);
430 # Computes z1**z2 = exp(z2 * log z1)).
433 my ($z1, $z2, $inverted) = @_;
435 return 1 if $z1 == 0 || $z2 == 1;
436 return 0 if $z2 == 0 && Re($z1) > 0;
438 return 1 if $z2 == 0 || $z1 == 1;
439 return 0 if $z1 == 0 && Re($z2) > 0;
441 my $w = $inverted ? &exp($z1 * &log($z2))
442 : &exp($z2 * &log($z1));
443 # If both arguments cartesian, return cartesian, else polar.
444 return $z1->{c_dirty} == 0 &&
445 (not ref $z2 or $z2->{c_dirty} == 0) ?
446 cplx(@{$w->cartesian}) : $w;
452 # Computes z1 <=> z2.
453 # Sorts on the real part first, then on the imaginary part. Thus 2-4i < 3+8i.
456 my ($z1, $z2, $inverted) = @_;
457 my ($re1, $im1) = ref $z1 ? @{$z1->cartesian} : ($z1, 0);
458 my ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
459 my $sgn = $inverted ? -1 : 1;
460 return $sgn * ($re1 <=> $re2) if $re1 != $re2;
461 return $sgn * ($im1 <=> $im2);
469 # (Required in addition to spaceship() because of NaNs.)
471 my ($z1, $z2, $inverted) = @_;
472 my ($re1, $im1) = ref $z1 ? @{$z1->cartesian} : ($z1, 0);
473 my ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
474 return $re1 == $re2 && $im1 == $im2 ? 1 : 0;
485 my ($r, $t) = @{$z->polar};
486 $t = ($t <= 0) ? $t + pi : $t - pi;
487 return (ref $z)->emake($r, $t);
489 my ($re, $im) = @{$z->cartesian};
490 return (ref $z)->make(-$re, -$im);
496 # Compute complex's conjugate.
501 my ($r, $t) = @{$z->polar};
502 return (ref $z)->emake($r, -$t);
504 my ($re, $im) = @{$z->cartesian};
505 return (ref $z)->make($re, -$im);
511 # Compute or set complex's norm (rho).
519 return CORE::abs($z);
523 $z->{'polar'} = [ $rho, ${$z->polar}[1] ];
528 return ${$z->polar}[0];
535 if ($$theta > pi()) { $$theta -= pit2 }
536 elsif ($$theta <= -pi()) { $$theta += pit2 }
542 # Compute or set complex's argument (theta).
545 my ($z, $theta) = @_;
546 return $z unless ref $z;
547 if (defined $theta) {
549 $z->{'polar'} = [ ${$z->polar}[0], $theta ];
553 $theta = ${$z->polar}[1];
564 # It is quite tempting to use wantarray here so that in list context
565 # sqrt() would return the two solutions. This, however, would
568 # print "sqrt(z) = ", sqrt($z), "\n";
570 # The two values would be printed side by side without no intervening
571 # whitespace, quite confusing.
572 # Therefore if you want the two solutions use the root().
576 my ($re, $im) = ref $z ? @{$z->cartesian} : ($z, 0);
577 return $re < 0 ? cplx(0, CORE::sqrt(-$re)) : CORE::sqrt($re)
579 my ($r, $t) = @{$z->polar};
580 return (ref $z)->emake(CORE::sqrt($r), $t/2);
586 # Compute cbrt(z) (cubic root).
588 # Why are we not returning three values? The same answer as for sqrt().
593 -CORE::exp(CORE::log(-$z)/3) :
594 ($z > 0 ? CORE::exp(CORE::log($z)/3): 0)
596 my ($r, $t) = @{$z->polar};
598 return (ref $z)->emake(CORE::exp(CORE::log($r)/3), $t/3);
607 my $mess = "Root $_[0] illegal, root rank must be positive integer.\n";
611 $mess .= "Died at $up[1] line $up[2].\n";
619 # Computes all nth root for z, returning an array whose size is n.
620 # `n' must be a positive integer.
622 # The roots are given by (for k = 0..n-1):
624 # z^(1/n) = r^(1/n) (cos ((t+2 k pi)/n) + i sin ((t+2 k pi)/n))
628 _rootbad($n) if ($n < 1 or int($n) != $n);
629 my ($r, $t) = ref $z ?
630 @{$z->polar} : (CORE::abs($z), $z >= 0 ? 0 : pi);
633 my $theta_inc = pit2 / $n;
634 my $rho = $r ** (1/$n);
636 my $cartesian = ref $z && $z->{c_dirty} == 0;
637 for ($k = 0, $theta = $t / $n; $k < $n; $k++, $theta += $theta_inc) {
638 my $w = cplxe($rho, $theta);
639 # Yes, $cartesian is loop invariant.
640 push @root, $cartesian ? cplx(@{$w->cartesian}) : $w;
648 # Return or set Re(z).
652 return $z unless ref $z;
654 $z->{'cartesian'} = [ $Re, ${$z->cartesian}[1] ];
658 return ${$z->cartesian}[0];
665 # Return or set Im(z).
669 return $z unless ref $z;
671 $z->{'cartesian'} = [ ${$z->cartesian}[0], $Im ];
675 return ${$z->cartesian}[1];
682 # Return or set rho(w).
685 Math::Complex::abs(@_);
691 # Return or set theta(w).
694 Math::Complex::arg(@_);
704 my ($x, $y) = @{$z->cartesian};
705 return (ref $z)->emake(CORE::exp($x), $y);
711 # Die on logarithm of zero.
714 my $mess = "$_[0]: Logarithm of zero.\n";
717 $mess .= "(Because in the definition of $_[0], the argument ";
718 $mess .= "$_[1] " unless ($_[1] eq '0');
724 $mess .= "Died at $up[1] line $up[2].\n";
737 _logofzero("log") if $z == 0;
738 return $z > 0 ? CORE::log($z) : cplx(CORE::log(-$z), pi);
740 my ($r, $t) = @{$z->polar};
741 _logofzero("log") if $r == 0;
742 if ($t > pi()) { $t -= pit2 }
743 elsif ($t <= -pi()) { $t += pit2 }
744 return (ref $z)->make(CORE::log($r), $t);
752 sub ln { Math::Complex::log(@_) }
761 return Math::Complex::log($_[0]) * uplog10;
767 # Compute logn(z,n) = log(z) / log(n)
771 $z = cplx($z, 0) unless ref $z;
772 my $logn = $LOGN{$n};
773 $logn = $LOGN{$n} = CORE::log($n) unless defined $logn; # Cache log(n)
774 return &log($z) / $logn;
780 # Compute cos(z) = (exp(iz) + exp(-iz))/2.
784 return CORE::cos($z) unless ref $z;
785 my ($x, $y) = @{$z->cartesian};
786 my $ey = CORE::exp($y);
787 my $sx = CORE::sin($x);
788 my $cx = CORE::cos($x);
789 my $ey_1 = $ey ? 1 / $ey : $Inf;
790 return (ref $z)->make($cx * ($ey + $ey_1)/2,
791 $sx * ($ey_1 - $ey)/2);
797 # Compute sin(z) = (exp(iz) - exp(-iz))/2.
801 return CORE::sin($z) unless ref $z;
802 my ($x, $y) = @{$z->cartesian};
803 my $ey = CORE::exp($y);
804 my $sx = CORE::sin($x);
805 my $cx = CORE::cos($x);
806 my $ey_1 = $ey ? 1 / $ey : $Inf;
807 return (ref $z)->make($sx * ($ey + $ey_1)/2,
808 $cx * ($ey - $ey_1)/2);
814 # Compute tan(z) = sin(z) / cos(z).
819 _divbyzero "tan($z)", "cos($z)" if $cz == 0;
820 return &sin($z) / $cz;
826 # Computes the secant sec(z) = 1 / cos(z).
831 _divbyzero "sec($z)", "cos($z)" if ($cz == 0);
838 # Computes the cosecant csc(z) = 1 / sin(z).
843 _divbyzero "csc($z)", "sin($z)" if ($sz == 0);
852 sub cosec { Math::Complex::csc(@_) }
857 # Computes cot(z) = cos(z) / sin(z).
862 _divbyzero "cot($z)", "sin($z)" if ($sz == 0);
863 return &cos($z) / $sz;
871 sub cotan { Math::Complex::cot(@_) }
876 # Computes the arc cosine acos(z) = -i log(z + sqrt(z*z-1)).
880 return CORE::atan2(CORE::sqrt(1-$z*$z), $z)
881 if (! ref $z) && CORE::abs($z) <= 1;
882 my ($x, $y) = ref $z ? @{$z->cartesian} : ($z, 0);
883 return 0 if $x == 1 && $y == 0;
884 my $t1 = CORE::sqrt(($x+1)*($x+1) + $y*$y);
885 my $t2 = CORE::sqrt(($x-1)*($x-1) + $y*$y);
886 my $alpha = ($t1 + $t2)/2;
887 my $beta = ($t1 - $t2)/2;
888 $alpha = 1 if $alpha < 1;
889 if ($beta > 1) { $beta = 1 }
890 elsif ($beta < -1) { $beta = -1 }
891 my $u = CORE::atan2(CORE::sqrt(1-$beta*$beta), $beta);
892 my $v = CORE::log($alpha + CORE::sqrt($alpha*$alpha-1));
893 $v = -$v if $y > 0 || ($y == 0 && $x < -1);
894 return __PACKAGE__->make($u, $v);
900 # Computes the arc sine asin(z) = -i log(iz + sqrt(1-z*z)).
904 return CORE::atan2($z, CORE::sqrt(1-$z*$z))
905 if (! ref $z) && CORE::abs($z) <= 1;
906 my ($x, $y) = ref $z ? @{$z->cartesian} : ($z, 0);
907 return 0 if $x == 0 && $y == 0;
908 my $t1 = CORE::sqrt(($x+1)*($x+1) + $y*$y);
909 my $t2 = CORE::sqrt(($x-1)*($x-1) + $y*$y);
910 my $alpha = ($t1 + $t2)/2;
911 my $beta = ($t1 - $t2)/2;
912 $alpha = 1 if $alpha < 1;
913 if ($beta > 1) { $beta = 1 }
914 elsif ($beta < -1) { $beta = -1 }
915 my $u = CORE::atan2($beta, CORE::sqrt(1-$beta*$beta));
916 my $v = -CORE::log($alpha + CORE::sqrt($alpha*$alpha-1));
917 $v = -$v if $y > 0 || ($y == 0 && $x < -1);
918 return __PACKAGE__->make($u, $v);
924 # Computes the arc tangent atan(z) = i/2 log((i+z) / (i-z)).
928 return CORE::atan2($z, 1) unless ref $z;
929 my ($x, $y) = ref $z ? @{$z->cartesian} : ($z, 0);
930 return 0 if $x == 0 && $y == 0;
931 _divbyzero "atan(i)" if ( $z == i);
932 _logofzero "atan(-i)" if (-$z == i); # -i is a bad file test...
933 my $log = &log((i + $z) / (i - $z));
940 # Computes the arc secant asec(z) = acos(1 / z).
944 _divbyzero "asec($z)", $z if ($z == 0);
951 # Computes the arc cosecant acsc(z) = asin(1 / z).
955 _divbyzero "acsc($z)", $z if ($z == 0);
964 sub acosec { Math::Complex::acsc(@_) }
969 # Computes the arc cotangent acot(z) = atan(1 / z)
973 _divbyzero "acot(0)" if $z == 0;
974 return ($z >= 0) ? CORE::atan2(1, $z) : CORE::atan2(-1, -$z)
976 _divbyzero "acot(i)" if ($z - i == 0);
977 _logofzero "acot(-i)" if ($z + i == 0);
986 sub acotan { Math::Complex::acot(@_) }
991 # Computes the hyperbolic cosine cosh(z) = (exp(z) + exp(-z))/2.
998 return $ex ? ($ex + 1/$ex)/2 : $Inf;
1000 my ($x, $y) = @{$z->cartesian};
1001 my $cy = CORE::cos($y);
1002 my $sy = CORE::cos($y);
1003 $ex = CORE::exp($x);
1004 my $ex_1 = $ex ? 1 / $ex : $Inf;
1005 return (ref $z)->make(CORE::cos($y) * ($ex + $ex_1)/2,
1006 CORE::sin($y) * ($ex - $ex_1)/2);
1012 # Computes the hyperbolic sine sinh(z) = (exp(z) - exp(-z))/2.
1018 return 0 if $z == 0;
1019 $ex = CORE::exp($z);
1020 return $ex ? ($ex - 1/$ex)/2 : "-$Inf";
1022 my ($x, $y) = @{$z->cartesian};
1023 my $cy = CORE::cos($y);
1024 my $sy = CORE::sin($y);
1025 $ex = CORE::exp($x);
1026 my $ex_1 = $ex ? 1 / $ex : $Inf;
1027 return (ref $z)->make($cy * ($ex - $ex_1)/2,
1028 $sy * ($ex + $ex_1)/2);
1034 # Computes the hyperbolic tangent tanh(z) = sinh(z) / cosh(z).
1039 _divbyzero "tanh($z)", "cosh($z)" if ($cz == 0);
1040 return sinh($z) / $cz;
1046 # Computes the hyperbolic secant sech(z) = 1 / cosh(z).
1051 _divbyzero "sech($z)", "cosh($z)" if ($cz == 0);
1058 # Computes the hyperbolic cosecant csch(z) = 1 / sinh(z).
1063 _divbyzero "csch($z)", "sinh($z)" if ($sz == 0);
1072 sub cosech { Math::Complex::csch(@_) }
1077 # Computes the hyperbolic cotangent coth(z) = cosh(z) / sinh(z).
1082 _divbyzero "coth($z)", "sinh($z)" if $sz == 0;
1083 return cosh($z) / $sz;
1091 sub cotanh { Math::Complex::coth(@_) }
1096 # Computes the arc hyperbolic cosine acosh(z) = log(z + sqrt(z*z-1)).
1103 my ($re, $im) = @{$z->cartesian};
1105 return CORE::log($re + CORE::sqrt($re*$re - 1))
1107 return cplx(0, CORE::atan2(CORE::sqrt(1 - $re*$re), $re))
1108 if CORE::abs($re) < 1;
1110 my $s = &sqrt($z*$z - 1);
1112 $t = 1/(2*$s) if $t == 0 || $t && &abs(cosh(&log($t)) - $z) > $eps;
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 $s = &sqrt($z*$z + 1);
1129 # Try Taylor series if looking bad.
1130 $t = 1/(2*$s) if $t == 0 || $t && &abs(sinh(&log($t)) - $z) > $eps;
1137 # Computes the arc hyperbolic tangent atanh(z) = 1/2 log((1+z) / (1-z)).
1142 return CORE::log((1 + $z)/(1 - $z))/2 if CORE::abs($z) < 1;
1145 _divbyzero 'atanh(1)', "1 - $z" if (1 - $z == 0);
1146 _logofzero 'atanh(-1)' if (1 + $z == 0);
1147 return 0.5 * &log((1 + $z) / (1 - $z));
1153 # Computes the hyperbolic arc secant asech(z) = acosh(1 / z).
1157 _divbyzero 'asech(0)', "$z" if ($z == 0);
1158 return acosh(1 / $z);
1164 # Computes the hyperbolic arc cosecant acsch(z) = asinh(1 / z).
1168 _divbyzero 'acsch(0)', $z if ($z == 0);
1169 return asinh(1 / $z);
1175 # Alias for acosh().
1177 sub acosech { Math::Complex::acsch(@_) }
1182 # Computes the arc hyperbolic cotangent acoth(z) = 1/2 log((1+z) / (z-1)).
1186 _divbyzero 'acoth(0)' if ($z == 0);
1188 return CORE::log(($z + 1)/($z - 1))/2 if CORE::abs($z) > 1;
1191 _divbyzero 'acoth(1)', "$z - 1" if ($z - 1 == 0);
1192 _logofzero 'acoth(-1)', "1 + $z" if (1 + $z == 0);
1193 return &log((1 + $z) / ($z - 1)) / 2;
1201 sub acotanh { Math::Complex::acoth(@_) }
1206 # Compute atan(z1/z2).
1209 my ($z1, $z2, $inverted) = @_;
1210 my ($re1, $im1, $re2, $im2);
1212 ($re1, $im1) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
1213 ($re2, $im2) = @{$z1->cartesian};
1215 ($re1, $im1) = @{$z1->cartesian};
1216 ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
1219 return CORE::atan2($re1, $re2) if $im1 == 0;
1220 return ($im1<=>0) * pip2 if $re2 == 0;
1222 my $w = atan($z1/$z2);
1223 my ($u, $v) = ref $w ? @{$w->cartesian} : ($w, 0);
1224 $u += pi if $re2 < 0;
1225 $u -= pit2 if $u > pi;
1226 return cplx($u, $v);
1233 # Set (get if no argument) the display format for all complex numbers that
1234 # don't happen to have overridden it via ->display_format
1236 # When called as an object method, this actually sets the display format for
1237 # the current object.
1239 # Valid object formats are 'c' and 'p' for cartesian and polar. The first
1240 # letter is used actually, so the type can be fully spelled out for clarity.
1242 sub display_format {
1244 my %display_format = %DISPLAY_FORMAT;
1246 if (ref $self) { # Called as an object method
1247 if (exists $self->{display_format}) {
1248 my %obj = %{$self->{display_format}};
1249 @display_format{keys %obj} = values %obj;
1252 $display_format{style} = shift;
1255 @display_format{keys %new} = values %new;
1257 } else { # Called as a class method
1259 $display_format{style} = $self;
1262 @display_format{keys %new} = values %new;
1267 if (defined $self) {
1268 $self->{display_format} = { %display_format };
1271 %{$self->{display_format}} :
1272 $self->{display_format}->{style};
1275 %DISPLAY_FORMAT = %display_format;
1279 $DISPLAY_FORMAT{style};
1285 # Show nicely formatted complex number under its cartesian or polar form,
1286 # depending on the current display format:
1288 # . If a specific display format has been recorded for this object, use it.
1289 # . Otherwise, use the generic current default for all complex numbers,
1290 # which is a package global variable.
1295 my $style = $z->display_format;
1297 $style = $DISPLAY_FORMAT{style} unless defined $style;
1299 return $z->stringify_polar if $style =~ /^p/i;
1300 return $z->stringify_cartesian;
1304 # ->stringify_cartesian
1306 # Stringify as a cartesian representation 'a+bi'.
1308 sub stringify_cartesian {
1310 my ($x, $y) = @{$z->cartesian};
1313 my %format = $z->display_format;
1314 my $format = $format{format};
1317 if ($x =~ /^NaN[QS]?$/i) {
1320 if ($x =~ /^-?$Inf$/oi) {
1323 $re = defined $format ? sprintf($format, $x) : $x;
1331 if ($y == 1) { $im = "" }
1332 elsif ($y == -1) { $im = "-" }
1333 elsif ($y =~ /^(NaN[QS]?)$/i) {
1336 if ($y =~ /^-?$Inf$/oi) {
1339 $im = defined $format ? sprintf($format, $y) : $y;
1352 } elsif ($y > 0 || $im =~ /^NaN[QS]?i$/i) {
1353 $str .= "+" if defined $re;
1356 } elsif (!defined $re) {
1367 # Stringify as a polar representation '[r,t]'.
1369 sub stringify_polar {
1371 my ($r, $t) = @{$z->polar};
1374 my %format = $z->display_format;
1375 my $format = $format{format};
1377 if ($t =~ /^NaN[QS]?$/i || $t =~ /^-?$Inf$/oi) {
1379 } elsif ($t == pi) {
1381 } elsif ($r == 0 || $t == 0) {
1382 $theta = defined $format ? sprintf($format, $t) : $t;
1385 return "[$r,$theta]" if defined $theta;
1388 # Try to identify pi/n and friends.
1391 $t -= int(CORE::abs($t) / pit2) * pit2;
1393 if ($format{polar_pretty_print}) {
1395 for $a (2, 3, 4, 6, 8, 12, 16, 24, 30, 32, 36, 48, 60, 64, 72) {
1397 if (int($b) == $b) {
1398 $b = $b < 0 ? "-" : "" if CORE::abs($b) == 1;
1399 $theta = "${b}pi/$a";
1405 if (defined $format) {
1406 $r = sprintf($format, $r);
1407 $theta = sprintf($format, $theta) unless defined $theta;
1409 $theta = $t unless defined $theta;
1412 return "[$r,$theta]";
1420 Math::Complex - complex numbers and associated mathematical functions
1426 $z = Math::Complex->make(5, 6);
1428 $j = cplxe(1, 2*pi/3);
1432 This package lets you create and manipulate complex numbers. By default,
1433 I<Perl> limits itself to real numbers, but an extra C<use> statement brings
1434 full complex support, along with a full set of mathematical functions
1435 typically associated with and/or extended to complex numbers.
1437 If you wonder what complex numbers are, they were invented to be able to solve
1438 the following equation:
1442 and by definition, the solution is noted I<i> (engineers use I<j> instead since
1443 I<i> usually denotes an intensity, but the name does not matter). The number
1444 I<i> is a pure I<imaginary> number.
1446 The arithmetics with pure imaginary numbers works just like you would expect
1447 it with real numbers... you just have to remember that
1453 5i + 7i = i * (5 + 7) = 12i
1454 4i - 3i = i * (4 - 3) = i
1459 Complex numbers are numbers that have both a real part and an imaginary
1460 part, and are usually noted:
1464 where C<a> is the I<real> part and C<b> is the I<imaginary> part. The
1465 arithmetic with complex numbers is straightforward. You have to
1466 keep track of the real and the imaginary parts, but otherwise the
1467 rules used for real numbers just apply:
1469 (4 + 3i) + (5 - 2i) = (4 + 5) + i(3 - 2) = 9 + i
1470 (2 + i) * (4 - i) = 2*4 + 4i -2i -i*i = 8 + 2i + 1 = 9 + 2i
1472 A graphical representation of complex numbers is possible in a plane
1473 (also called the I<complex plane>, but it's really a 2D plane).
1478 is the point whose coordinates are (a, b). Actually, it would
1479 be the vector originating from (0, 0) to (a, b). It follows that the addition
1480 of two complex numbers is a vectorial addition.
1482 Since there is a bijection between a point in the 2D plane and a complex
1483 number (i.e. the mapping is unique and reciprocal), a complex number
1484 can also be uniquely identified with polar coordinates:
1488 where C<rho> is the distance to the origin, and C<theta> the angle between
1489 the vector and the I<x> axis. There is a notation for this using the
1490 exponential form, which is:
1492 rho * exp(i * theta)
1494 where I<i> is the famous imaginary number introduced above. Conversion
1495 between this form and the cartesian form C<a + bi> is immediate:
1497 a = rho * cos(theta)
1498 b = rho * sin(theta)
1500 which is also expressed by this formula:
1502 z = rho * exp(i * theta) = rho * (cos theta + i * sin theta)
1504 In other words, it's the projection of the vector onto the I<x> and I<y>
1505 axes. Mathematicians call I<rho> the I<norm> or I<modulus> and I<theta>
1506 the I<argument> of the complex number. The I<norm> of C<z> will be
1509 The polar notation (also known as the trigonometric
1510 representation) is much more handy for performing multiplications and
1511 divisions of complex numbers, whilst the cartesian notation is better
1512 suited for additions and subtractions. Real numbers are on the I<x>
1513 axis, and therefore I<theta> is zero or I<pi>.
1515 All the common operations that can be performed on a real number have
1516 been defined to work on complex numbers as well, and are merely
1517 I<extensions> of the operations defined on real numbers. This means
1518 they keep their natural meaning when there is no imaginary part, provided
1519 the number is within their definition set.
1521 For instance, the C<sqrt> routine which computes the square root of
1522 its argument is only defined for non-negative real numbers and yields a
1523 non-negative real number (it is an application from B<R+> to B<R+>).
1524 If we allow it to return a complex number, then it can be extended to
1525 negative real numbers to become an application from B<R> to B<C> (the
1526 set of complex numbers):
1528 sqrt(x) = x >= 0 ? sqrt(x) : sqrt(-x)*i
1530 It can also be extended to be an application from B<C> to B<C>,
1531 whilst its restriction to B<R> behaves as defined above by using
1532 the following definition:
1534 sqrt(z = [r,t]) = sqrt(r) * exp(i * t/2)
1536 Indeed, a negative real number can be noted C<[x,pi]> (the modulus
1537 I<x> is always non-negative, so C<[x,pi]> is really C<-x>, a negative
1538 number) and the above definition states that
1540 sqrt([x,pi]) = sqrt(x) * exp(i*pi/2) = [sqrt(x),pi/2] = sqrt(x)*i
1542 which is exactly what we had defined for negative real numbers above.
1543 The C<sqrt> returns only one of the solutions: if you want the both,
1544 use the C<root> function.
1546 All the common mathematical functions defined on real numbers that
1547 are extended to complex numbers share that same property of working
1548 I<as usual> when the imaginary part is zero (otherwise, it would not
1549 be called an extension, would it?).
1551 A I<new> operation possible on a complex number that is
1552 the identity for real numbers is called the I<conjugate>, and is noted
1553 with an horizontal bar above the number, or C<~z> here.
1560 z * ~z = (a + bi) * (a - bi) = a*a + b*b
1562 We saw that the norm of C<z> was noted C<abs(z)> and was defined as the
1563 distance to the origin, also known as:
1565 rho = abs(z) = sqrt(a*a + b*b)
1569 z * ~z = abs(z) ** 2
1571 If z is a pure real number (i.e. C<b == 0>), then the above yields:
1575 which is true (C<abs> has the regular meaning for real number, i.e. stands
1576 for the absolute value). This example explains why the norm of C<z> is
1577 noted C<abs(z)>: it extends the C<abs> function to complex numbers, yet
1578 is the regular C<abs> we know when the complex number actually has no
1579 imaginary part... This justifies I<a posteriori> our use of the C<abs>
1580 notation for the norm.
1584 Given the following notations:
1586 z1 = a + bi = r1 * exp(i * t1)
1587 z2 = c + di = r2 * exp(i * t2)
1588 z = <any complex or real number>
1590 the following (overloaded) operations are supported on complex numbers:
1592 z1 + z2 = (a + c) + i(b + d)
1593 z1 - z2 = (a - c) + i(b - d)
1594 z1 * z2 = (r1 * r2) * exp(i * (t1 + t2))
1595 z1 / z2 = (r1 / r2) * exp(i * (t1 - t2))
1596 z1 ** z2 = exp(z2 * log z1)
1598 abs(z) = r1 = sqrt(a*a + b*b)
1599 sqrt(z) = sqrt(r1) * exp(i * t/2)
1600 exp(z) = exp(a) * exp(i * b)
1601 log(z) = log(r1) + i*t
1602 sin(z) = 1/2i (exp(i * z1) - exp(-i * z))
1603 cos(z) = 1/2 (exp(i * z1) + exp(-i * z))
1604 atan2(z1, z2) = atan(z1/z2)
1606 The following extra operations are supported on both real and complex
1614 cbrt(z) = z ** (1/3)
1615 log10(z) = log(z) / log(10)
1616 logn(z, n) = log(z) / log(n)
1618 tan(z) = sin(z) / cos(z)
1624 asin(z) = -i * log(i*z + sqrt(1-z*z))
1625 acos(z) = -i * log(z + i*sqrt(1-z*z))
1626 atan(z) = i/2 * log((i+z) / (i-z))
1628 acsc(z) = asin(1 / z)
1629 asec(z) = acos(1 / z)
1630 acot(z) = atan(1 / z) = -i/2 * log((i+z) / (z-i))
1632 sinh(z) = 1/2 (exp(z) - exp(-z))
1633 cosh(z) = 1/2 (exp(z) + exp(-z))
1634 tanh(z) = sinh(z) / cosh(z) = (exp(z) - exp(-z)) / (exp(z) + exp(-z))
1636 csch(z) = 1 / sinh(z)
1637 sech(z) = 1 / cosh(z)
1638 coth(z) = 1 / tanh(z)
1640 asinh(z) = log(z + sqrt(z*z+1))
1641 acosh(z) = log(z + sqrt(z*z-1))
1642 atanh(z) = 1/2 * log((1+z) / (1-z))
1644 acsch(z) = asinh(1 / z)
1645 asech(z) = acosh(1 / z)
1646 acoth(z) = atanh(1 / z) = 1/2 * log((1+z) / (z-1))
1648 I<arg>, I<abs>, I<log>, I<csc>, I<cot>, I<acsc>, I<acot>, I<csch>,
1649 I<coth>, I<acosech>, I<acotanh>, have aliases I<rho>, I<theta>, I<ln>,
1650 I<cosec>, I<cotan>, I<acosec>, I<acotan>, I<cosech>, I<cotanh>,
1651 I<acosech>, I<acotanh>, respectively. C<Re>, C<Im>, C<arg>, C<abs>,
1652 C<rho>, and C<theta> can be used also also mutators. The C<cbrt>
1653 returns only one of the solutions: if you want all three, use the
1656 The I<root> function is available to compute all the I<n>
1657 roots of some complex, where I<n> is a strictly positive integer.
1658 There are exactly I<n> such roots, returned as a list. Getting the
1659 number mathematicians call C<j> such that:
1663 is a simple matter of writing:
1665 $j = ((root(1, 3))[1];
1667 The I<k>th root for C<z = [r,t]> is given by:
1669 (root(z, n))[k] = r**(1/n) * exp(i * (t + 2*k*pi)/n)
1671 The I<spaceship> comparison operator, E<lt>=E<gt>, is also defined. In
1672 order to ensure its restriction to real numbers is conform to what you
1673 would expect, the comparison is run on the real part of the complex
1674 number first, and imaginary parts are compared only when the real
1679 To create a complex number, use either:
1681 $z = Math::Complex->make(3, 4);
1684 if you know the cartesian form of the number, or
1688 if you like. To create a number using the polar form, use either:
1690 $z = Math::Complex->emake(5, pi/3);
1691 $x = cplxe(5, pi/3);
1693 instead. The first argument is the modulus, the second is the angle
1694 (in radians, the full circle is 2*pi). (Mnemonic: C<e> is used as a
1695 notation for complex numbers in the polar form).
1697 It is possible to write:
1699 $x = cplxe(-3, pi/4);
1701 but that will be silently converted into C<[3,-3pi/4]>, since the
1702 modulus must be non-negative (it represents the distance to the origin
1703 in the complex plane).
1705 It is also possible to have a complex number as either argument of
1706 either the C<make> or C<emake>: the appropriate component of
1707 the argument will be used.
1712 =head1 STRINGIFICATION
1714 When printed, a complex number is usually shown under its cartesian
1715 style I<a+bi>, but there are legitimate cases where the polar style
1716 I<[r,t]> is more appropriate.
1718 By calling the class method C<Math::Complex::display_format> and
1719 supplying either C<"polar"> or C<"cartesian"> as an argument, you
1720 override the default display style, which is C<"cartesian">. Not
1721 supplying any argument returns the current settings.
1723 This default can be overridden on a per-number basis by calling the
1724 C<display_format> method instead. As before, not supplying any argument
1725 returns the current display style for this number. Otherwise whatever you
1726 specify will be the new display style for I<this> particular number.
1732 Math::Complex::display_format('polar');
1733 $j = (root(1, 3))[1];
1734 print "j = $j\n"; # Prints "j = [1,2pi/3]"
1735 $j->display_format('cartesian');
1736 print "j = $j\n"; # Prints "j = -0.5+0.866025403784439i"
1738 The polar style attempts to emphasize arguments like I<k*pi/n>
1739 (where I<n> is a positive integer and I<k> an integer within [-9,+9]),
1740 this is called I<polar pretty-printing>.
1742 =head2 CHANGED IN PERL 5.6
1744 The C<display_format> class method and the corresponding
1745 C<display_format> object method can now be called using
1746 a parameter hash instead of just a one parameter.
1748 The old display format style, which can have values C<"cartesian"> or
1749 C<"polar">, can be changed using the C<"style"> parameter. (The one
1750 parameter calling convention also still works.)
1752 There are two new display parameters.
1754 The first one is C<"format">, which is a sprintf()-style format
1755 string to be used for both parts of the complex number(s). The
1756 default is C<undef>, which corresponds usually (this is somewhat
1757 system-dependent) to C<"%.15g">. You can revert to the default by
1758 setting the format string to C<undef>.
1760 # the $j from the above example
1762 $j->display_format('format' => '%.5f');
1763 print "j = $j\n"; # Prints "j = -0.50000+0.86603i"
1764 $j->display_format('format' => '%.6f');
1765 print "j = $j\n"; # Prints "j = -0.5+0.86603i"
1767 Notice that this affects also the return values of the
1768 C<display_format> methods: in list context the whole parameter hash
1769 will be returned, as opposed to only the style parameter value. If
1770 you want to know the whole truth for a complex number, you must call
1771 both the class method and the object method:
1773 The second new display parameter is C<"polar_pretty_print">, which can
1774 be set to true or false, the default being true. See the previous
1775 section for what this means.
1779 Thanks to overloading, the handling of arithmetics with complex numbers
1780 is simple and almost transparent.
1782 Here are some examples:
1786 $j = cplxe(1, 2*pi/3); # $j ** 3 == 1
1787 print "j = $j, j**3 = ", $j ** 3, "\n";
1788 print "1 + j + j**2 = ", 1 + $j + $j**2, "\n";
1790 $z = -16 + 0*i; # Force it to be a complex
1791 print "sqrt($z) = ", sqrt($z), "\n";
1793 $k = exp(i * 2*pi/3);
1794 print "$j - $k = ", $j - $k, "\n";
1796 $z->Re(3); # Re, Im, arg, abs,
1797 $j->arg(2); # (the last two aka rho, theta)
1798 # can be used also as mutators.
1800 =head1 ERRORS DUE TO DIVISION BY ZERO OR LOGARITHM OF ZERO
1802 The division (/) and the following functions
1808 atanh asech acsch acoth
1810 cannot be computed for all arguments because that would mean dividing
1811 by zero or taking logarithm of zero. These situations cause fatal
1812 runtime errors looking like this
1814 cot(0): Division by zero.
1815 (Because in the definition of cot(0), the divisor sin(0) is 0)
1820 atanh(-1): Logarithm of zero.
1823 For the C<csc>, C<cot>, C<asec>, C<acsc>, C<acot>, C<csch>, C<coth>,
1824 C<asech>, C<acsch>, the argument cannot be C<0> (zero). For the the
1825 logarithmic functions and the C<atanh>, C<acoth>, the argument cannot
1826 be C<1> (one). For the C<atanh>, C<acoth>, the argument cannot be
1827 C<-1> (minus one). For the C<atan>, C<acot>, the argument cannot be
1828 C<i> (the imaginary unit). For the C<atan>, C<acoth>, the argument
1829 cannot be C<-i> (the negative imaginary unit). For the C<tan>,
1830 C<sec>, C<tanh>, the argument cannot be I<pi/2 + k * pi>, where I<k>
1833 Note that because we are operating on approximations of real numbers,
1834 these errors can happen when merely `too close' to the singularities
1835 listed above. For example C<tan(2*atan2(1,1)+1e-15)> will die of
1838 =head1 ERRORS DUE TO INDIGESTIBLE ARGUMENTS
1840 The C<make> and C<emake> accept both real and complex arguments.
1841 When they cannot recognize the arguments they will die with error
1842 messages like the following
1844 Math::Complex::make: Cannot take real part of ...
1845 Math::Complex::make: Cannot take real part of ...
1846 Math::Complex::emake: Cannot take rho of ...
1847 Math::Complex::emake: Cannot take theta of ...
1851 Saying C<use Math::Complex;> exports many mathematical routines in the
1852 caller environment and even overrides some (C<sqrt>, C<log>).
1853 This is construed as a feature by the Authors, actually... ;-)
1855 All routines expect to be given real or complex numbers. Don't attempt to
1856 use BigFloat, since Perl has currently no rule to disambiguate a '+'
1857 operation (for instance) between two overloaded entities.
1859 In Cray UNICOS there is some strange numerical instability that results
1860 in root(), cos(), sin(), cosh(), sinh(), losing accuracy fast. Beware.
1861 The bug may be in UNICOS math libs, in UNICOS C compiler, in Math::Complex.
1862 Whatever it is, it does not manifest itself anywhere else where Perl runs.
1866 Raphael Manfredi <F<Raphael_Manfredi@pobox.com>> and
1867 Jarkko Hietaniemi <F<jhi@iki.fi>>.
1869 Extensive patches by Daniel S. Lewart <F<d-lewart@uiuc.edu>>.