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
13 use vars qw($VERSION @ISA @EXPORT %EXPORT_TAGS);
15 my ( $i, $ip2, %logn );
17 $VERSION = sprintf("%s", q$Id: Complex.pm,v 1.25 1998/02/05 16:07:37 jhi Exp $ =~ /(\d+\.\d+)/);
24 csc cosec sec cot cotan
26 acsc acosec asec acot acotan
28 csch cosech sech coth cotanh
30 acsch acosech asech acoth acotanh
68 my $package = 'Math::Complex'; # Package name
69 my $display = 'cartesian'; # Default display format
70 my $eps = 1e-14; # Epsilon
73 # Object attributes (internal):
74 # cartesian [real, imaginary] -- cartesian form
75 # polar [rho, theta] -- polar form
76 # c_dirty cartesian form not up-to-date
77 # p_dirty polar form not up-to-date
78 # display display format (package's global when not set)
81 # Die on bad *make() arguments.
84 die "@{[(caller(1))[3]]}: Cannot take $_[0] of $_[1].\n";
90 # Create a new complex number (cartesian form)
93 my $self = bless {}, shift;
97 if ( $rre eq ref $self ) {
100 _cannot_make("real part", $rre);
105 if ( $rim eq ref $self ) {
108 _cannot_make("imaginary part", $rim);
111 $self->{'cartesian'} = [ $re, $im ];
112 $self->{c_dirty} = 0;
113 $self->{p_dirty} = 1;
114 $self->display_format('cartesian');
121 # Create a new complex number (exponential form)
124 my $self = bless {}, shift;
125 my ($rho, $theta) = @_;
128 if ( $rrh eq ref $self ) {
131 _cannot_make("rho", $rrh);
134 my $rth = ref $theta;
136 if ( $rth eq ref $self ) {
137 $theta = theta($theta);
139 _cannot_make("theta", $rth);
144 $theta = ($theta <= 0) ? $theta + pi() : $theta - pi();
146 $self->{'polar'} = [$rho, $theta];
147 $self->{p_dirty} = 0;
148 $self->{c_dirty} = 1;
149 $self->display_format('polar');
153 sub new { &make } # For backward compatibility only.
158 # Creates a complex number from a (re, im) tuple.
159 # This avoids the burden of writing Math::Complex->make(re, im).
163 return $package->make($re, defined $im ? $im : 0);
169 # Creates a complex number from a (rho, theta) tuple.
170 # This avoids the burden of writing Math::Complex->emake(rho, theta).
173 my ($rho, $theta) = @_;
174 return $package->emake($rho, defined $theta ? $theta : 0);
180 # The number defined as pi = 180 degrees
182 use constant pi => 4 * atan2(1, 1);
189 use constant pit2 => 2 * pi;
196 use constant pip2 => pi / 2;
201 # One degree in radians, used in stringify_polar.
204 use constant deg1 => pi / 180;
211 use constant uplog10 => 1 / log(10);
216 # The number defined as i*i = -1;
221 $i->{'cartesian'} = [0, 1];
222 $i->{'polar'} = [1, pip2];
229 # Attribute access/set routines
232 sub cartesian {$_[0]->{c_dirty} ?
233 $_[0]->update_cartesian : $_[0]->{'cartesian'}}
234 sub polar {$_[0]->{p_dirty} ?
235 $_[0]->update_polar : $_[0]->{'polar'}}
237 sub set_cartesian { $_[0]->{p_dirty}++; $_[0]->{'cartesian'} = $_[1] }
238 sub set_polar { $_[0]->{c_dirty}++; $_[0]->{'polar'} = $_[1] }
243 # Recompute and return the cartesian form, given accurate polar form.
245 sub update_cartesian {
247 my ($r, $t) = @{$self->{'polar'}};
248 $self->{c_dirty} = 0;
249 return $self->{'cartesian'} = [$r * cos $t, $r * sin $t];
256 # Recompute and return the polar form, given accurate cartesian form.
260 my ($x, $y) = @{$self->{'cartesian'}};
261 $self->{p_dirty} = 0;
262 return $self->{'polar'} = [0, 0] if $x == 0 && $y == 0;
263 return $self->{'polar'} = [sqrt($x*$x + $y*$y), atan2($y, $x)];
272 my ($z1, $z2, $regular) = @_;
273 my ($re1, $im1) = @{$z1->cartesian};
274 $z2 = cplx($z2) unless ref $z2;
275 my ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
276 unless (defined $regular) {
277 $z1->set_cartesian([$re1 + $re2, $im1 + $im2]);
280 return (ref $z1)->make($re1 + $re2, $im1 + $im2);
289 my ($z1, $z2, $inverted) = @_;
290 my ($re1, $im1) = @{$z1->cartesian};
291 $z2 = cplx($z2) unless ref $z2;
292 my ($re2, $im2) = @{$z2->cartesian};
293 unless (defined $inverted) {
294 $z1->set_cartesian([$re1 - $re2, $im1 - $im2]);
298 (ref $z1)->make($re2 - $re1, $im2 - $im1) :
299 (ref $z1)->make($re1 - $re2, $im1 - $im2);
309 my ($z1, $z2, $regular) = @_;
310 if ($z1->{p_dirty} == 0 and ref $z2 and $z2->{p_dirty} == 0) {
311 # if both polar better use polar to avoid rounding errors
312 my ($r1, $t1) = @{$z1->polar};
313 my ($r2, $t2) = @{$z2->polar};
315 if ($t > pi()) { $t -= pit2 }
316 elsif ($t <= -pi()) { $t += pit2 }
317 unless (defined $regular) {
318 $z1->set_polar([$r1 * $r2, $t]);
321 return (ref $z1)->emake($r1 * $r2, $t);
323 my ($x1, $y1) = @{$z1->cartesian};
325 my ($x2, $y2) = @{$z2->cartesian};
326 return (ref $z1)->make($x1*$x2-$y1*$y2, $x1*$y2+$y1*$x2);
328 return (ref $z1)->make($x1*$z2, $y1*$z2);
336 # Die on division by zero.
339 my $mess = "$_[0]: Division by zero.\n";
342 $mess .= "(Because in the definition of $_[0], the divisor ";
343 $mess .= "$_[1] " unless ($_[1] eq '0');
349 $mess .= "Died at $up[1] line $up[2].\n";
360 my ($z1, $z2, $inverted) = @_;
361 if ($z1->{p_dirty} == 0 and ref $z2 and $z2->{p_dirty} == 0) {
362 # if both polar better use polar to avoid rounding errors
363 my ($r1, $t1) = @{$z1->polar};
364 my ($r2, $t2) = @{$z2->polar};
367 _divbyzero "$z2/0" if ($r1 == 0);
369 if ($t > pi()) { $t -= pit2 }
370 elsif ($t <= -pi()) { $t += pit2 }
371 return (ref $z1)->emake($r2 / $r1, $t);
373 _divbyzero "$z1/0" if ($r2 == 0);
375 if ($t > pi()) { $t -= pit2 }
376 elsif ($t <= -pi()) { $t += pit2 }
377 return (ref $z1)->emake($r1 / $r2, $t);
382 ($x2, $y2) = @{$z1->cartesian};
383 $d = $x2*$x2 + $y2*$y2;
384 _divbyzero "$z2/0" if $d == 0;
385 return (ref $z1)->make(($x2*$z2)/$d, -($y2*$z2)/$d);
387 my ($x1, $y1) = @{$z1->cartesian};
389 ($x2, $y2) = @{$z2->cartesian};
390 $d = $x2*$x2 + $y2*$y2;
391 _divbyzero "$z1/0" if $d == 0;
392 my $u = ($x1*$x2 + $y1*$y2)/$d;
393 my $v = ($y1*$x2 - $x1*$y2)/$d;
394 return (ref $z1)->make($u, $v);
396 _divbyzero "$z1/0" if $z2 == 0;
397 return (ref $z1)->make($x1/$z2, $y1/$z2);
406 # Die on zero raised to the zeroth.
409 my $mess = "The zero raised to the zeroth power is not defined.\n";
413 $mess .= "Died at $up[1] line $up[2].\n";
421 # Computes z1**z2 = exp(z2 * log z1)).
424 my ($z1, $z2, $inverted) = @_;
427 _zerotozero if ($z1z and $z2z);
430 return 1 if ($z1z or $z2 == 1);
433 return 1 if ($z2z or $z1 == 1);
435 my $w = $inverted ? exp($z1 * log $z2) : exp($z2 * log $z1);
436 # If both arguments cartesian, return cartesian, else polar.
437 return $z1->{c_dirty} == 0 &&
438 (not ref $z2 or $z2->{c_dirty} == 0) ?
439 cplx(@{$w->cartesian}) : $w;
445 # Computes z1 <=> z2.
446 # Sorts on the real part first, then on the imaginary part. Thus 2-4i > 3+8i.
449 my ($z1, $z2, $inverted) = @_;
450 my ($re1, $im1) = ref $z1 ? @{$z1->cartesian} : ($z1, 0);
451 my ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
452 my $sgn = $inverted ? -1 : 1;
453 return $sgn * ($re1 <=> $re2) if $re1 != $re2;
454 return $sgn * ($im1 <=> $im2);
465 my ($r, $t) = @{$z->polar};
466 $t = ($t <= 0) ? $t + pi : $t - pi;
467 return (ref $z)->emake($r, $t);
469 my ($re, $im) = @{$z->cartesian};
470 return (ref $z)->make(-$re, -$im);
476 # Compute complex's conjugate.
481 my ($r, $t) = @{$z->polar};
482 return (ref $z)->emake($r, -$t);
484 my ($re, $im) = @{$z->cartesian};
485 return (ref $z)->make($re, -$im);
491 # Compute or set complex's norm (rho).
495 return $z unless ref $z;
497 $z->{'polar'} = [ $rho, ${$z->polar}[1] ];
502 return ${$z->polar}[0];
509 if ($$theta > pi()) { $$theta -= pit2 }
510 elsif ($$theta <= -pi()) { $$theta += pit2 }
516 # Compute or set complex's argument (theta).
519 my ($z, $theta) = @_;
520 return $z unless ref $z;
521 if (defined $theta) {
523 $z->{'polar'} = [ ${$z->polar}[0], $theta ];
527 $theta = ${$z->polar}[1];
538 # It is quite tempting to use wantarray here so that in list context
539 # sqrt() would return the two solutions. This, however, would
542 # print "sqrt(z) = ", sqrt($z), "\n";
544 # The two values would be printed side by side without no intervening
545 # whitespace, quite confusing.
546 # Therefore if you want the two solutions use the root().
550 my ($re, $im) = ref $z ? @{$z->cartesian} : ($z, 0);
551 return $re < 0 ? cplx(0, sqrt(-$re)) : sqrt($re) if $im == 0;
552 my ($r, $t) = @{$z->polar};
553 return (ref $z)->emake(sqrt($r), $t/2);
559 # Compute cbrt(z) (cubic root).
561 # Why are we not returning three values? The same answer as for sqrt().
565 return $z < 0 ? -exp(log(-$z)/3) : ($z > 0 ? exp(log($z)/3): 0)
567 my ($r, $t) = @{$z->polar};
568 return (ref $z)->emake(exp(log($r)/3), $t/3);
577 my $mess = "Root $_[0] not defined, root must be positive integer.\n";
581 $mess .= "Died at $up[1] line $up[2].\n";
589 # Computes all nth root for z, returning an array whose size is n.
590 # `n' must be a positive integer.
592 # The roots are given by (for k = 0..n-1):
594 # z^(1/n) = r^(1/n) (cos ((t+2 k pi)/n) + i sin ((t+2 k pi)/n))
598 _rootbad($n) if ($n < 1 or int($n) != $n);
599 my ($r, $t) = ref $z ? @{$z->polar} : (abs($z), $z >= 0 ? 0 : pi);
602 my $theta_inc = pit2 / $n;
603 my $rho = $r ** (1/$n);
605 my $cartesian = ref $z && $z->{c_dirty} == 0;
606 for ($k = 0, $theta = $t / $n; $k < $n; $k++, $theta += $theta_inc) {
607 my $w = cplxe($rho, $theta);
608 # Yes, $cartesian is loop invariant.
609 push @root, $cartesian ? cplx(@{$w->cartesian}) : $w;
617 # Return or set Re(z).
621 return $z unless ref $z;
623 $z->{'cartesian'} = [ $Re, ${$z->cartesian}[1] ];
627 return ${$z->cartesian}[0];
634 # Return or set Im(z).
638 return $z unless ref $z;
640 $z->{'cartesian'} = [ ${$z->cartesian}[0], $Im ];
644 return ${$z->cartesian}[1];
651 # Return or set rho(w).
654 Math::Complex::abs(@_);
660 # Return or set theta(w).
663 Math::Complex::arg(@_);
673 my ($x, $y) = @{$z->cartesian};
674 return (ref $z)->emake(exp($x), $y);
680 # Die on logarithm of zero.
683 my $mess = "$_[0]: Logarithm of zero.\n";
686 $mess .= "(Because in the definition of $_[0], the argument ";
687 $mess .= "$_[1] " unless ($_[1] eq '0');
693 $mess .= "Died at $up[1] line $up[2].\n";
706 _logofzero("log") if $z == 0;
707 return $z > 0 ? log($z) : cplx(log(-$z), pi);
709 my ($r, $t) = @{$z->polar};
710 _logofzero("log") if $r == 0;
711 if ($t > pi()) { $t -= pit2 }
712 elsif ($t <= -pi()) { $t += pit2 }
713 return (ref $z)->make(log($r), $t);
721 sub ln { Math::Complex::log(@_) }
730 return Math::Complex::log($_[0]) * uplog10;
736 # Compute logn(z,n) = log(z) / log(n)
740 $z = cplx($z, 0) unless ref $z;
741 my $logn = $logn{$n};
742 $logn = $logn{$n} = log($n) unless defined $logn; # Cache log(n)
743 return log($z) / $logn;
749 # Compute cos(z) = (exp(iz) + exp(-iz))/2.
753 my ($x, $y) = @{$z->cartesian};
756 return (ref $z)->make(cos($x) * ($ey + $ey_1)/2,
757 sin($x) * ($ey_1 - $ey)/2);
763 # Compute sin(z) = (exp(iz) - exp(-iz))/2.
767 my ($x, $y) = @{$z->cartesian};
770 return (ref $z)->make(sin($x) * ($ey + $ey_1)/2,
771 cos($x) * ($ey - $ey_1)/2);
777 # Compute tan(z) = sin(z) / cos(z).
782 _divbyzero "tan($z)", "cos($z)" if (abs($cz) < $eps);
783 return sin($z) / $cz;
789 # Computes the secant sec(z) = 1 / cos(z).
794 _divbyzero "sec($z)", "cos($z)" if ($cz == 0);
801 # Computes the cosecant csc(z) = 1 / sin(z).
806 _divbyzero "csc($z)", "sin($z)" if ($sz == 0);
815 sub cosec { Math::Complex::csc(@_) }
820 # Computes cot(z) = cos(z) / sin(z).
825 _divbyzero "cot($z)", "sin($z)" if ($sz == 0);
826 return cos($z) / $sz;
834 sub cotan { Math::Complex::cot(@_) }
839 # Computes the arc cosine acos(z) = -i log(z + sqrt(z*z-1)).
843 return atan2(sqrt(1-$z*$z), $z) if (! ref $z) && abs($z) <= 1;
844 my ($x, $y) = ref $z ? @{$z->cartesian} : ($z, 0);
845 my $t1 = sqrt(($x+1)*($x+1) + $y*$y);
846 my $t2 = sqrt(($x-1)*($x-1) + $y*$y);
847 my $alpha = ($t1 + $t2)/2;
848 my $beta = ($t1 - $t2)/2;
849 $alpha = 1 if $alpha < 1;
850 if ($beta > 1) { $beta = 1 }
851 elsif ($beta < -1) { $beta = -1 }
852 my $u = atan2(sqrt(1-$beta*$beta), $beta);
853 my $v = log($alpha + sqrt($alpha*$alpha-1));
854 $v = -$v if $y > 0 || ($y == 0 && $x < -1);
855 return $package->make($u, $v);
861 # Computes the arc sine asin(z) = -i log(iz + sqrt(1-z*z)).
865 return atan2($z, sqrt(1-$z*$z)) if (! ref $z) && abs($z) <= 1;
866 my ($x, $y) = ref $z ? @{$z->cartesian} : ($z, 0);
867 my $t1 = sqrt(($x+1)*($x+1) + $y*$y);
868 my $t2 = sqrt(($x-1)*($x-1) + $y*$y);
869 my $alpha = ($t1 + $t2)/2;
870 my $beta = ($t1 - $t2)/2;
871 $alpha = 1 if $alpha < 1;
872 if ($beta > 1) { $beta = 1 }
873 elsif ($beta < -1) { $beta = -1 }
874 my $u = atan2($beta, sqrt(1-$beta*$beta));
875 my $v = -log($alpha + sqrt($alpha*$alpha-1));
876 $v = -$v if $y > 0 || ($y == 0 && $x < -1);
877 return $package->make($u, $v);
883 # Computes the arc tangent atan(z) = i/2 log((i+z) / (i-z)).
887 return atan2($z, 1) unless ref $z;
888 _divbyzero "atan(i)" if ( $z == i);
889 _divbyzero "atan(-i)" if (-$z == i);
890 my $log = log((i + $z) / (i - $z));
891 $ip2 = 0.5 * i unless defined $ip2;
898 # Computes the arc secant asec(z) = acos(1 / z).
902 _divbyzero "asec($z)", $z if ($z == 0);
909 # Computes the arc cosecant acsc(z) = asin(1 / z).
913 _divbyzero "acsc($z)", $z if ($z == 0);
922 sub acosec { Math::Complex::acsc(@_) }
927 # Computes the arc cotangent acot(z) = atan(1 / z)
931 _divbyzero "acot(0)" if (abs($z) < $eps);
932 return ($z >= 0) ? atan2(1, $z) : atan2(-1, -$z) unless ref $z;
933 _divbyzero "acot(i)" if (abs($z - i) < $eps);
934 _logofzero "acot(-i)" if (abs($z + i) < $eps);
943 sub acotan { Math::Complex::acot(@_) }
948 # Computes the hyperbolic cosine cosh(z) = (exp(z) + exp(-z))/2.
955 return ($ex + 1/$ex)/2;
957 my ($x, $y) = @{$z->cartesian};
960 return (ref $z)->make(cos($y) * ($ex + $ex_1)/2,
961 sin($y) * ($ex - $ex_1)/2);
967 # Computes the hyperbolic sine sinh(z) = (exp(z) - exp(-z))/2.
974 return ($ex - 1/$ex)/2;
976 my ($x, $y) = @{$z->cartesian};
979 return (ref $z)->make(cos($y) * ($ex - $ex_1)/2,
980 sin($y) * ($ex + $ex_1)/2);
986 # Computes the hyperbolic tangent tanh(z) = sinh(z) / cosh(z).
991 _divbyzero "tanh($z)", "cosh($z)" if ($cz == 0);
992 return sinh($z) / $cz;
998 # Computes the hyperbolic secant sech(z) = 1 / cosh(z).
1003 _divbyzero "sech($z)", "cosh($z)" if ($cz == 0);
1010 # Computes the hyperbolic cosecant csch(z) = 1 / sinh(z).
1015 _divbyzero "csch($z)", "sinh($z)" if ($sz == 0);
1024 sub cosech { Math::Complex::csch(@_) }
1029 # Computes the hyperbolic cotangent coth(z) = cosh(z) / sinh(z).
1034 _divbyzero "coth($z)", "sinh($z)" if ($sz == 0);
1035 return cosh($z) / $sz;
1043 sub cotanh { Math::Complex::coth(@_) }
1048 # Computes the arc hyperbolic cosine acosh(z) = log(z + sqrt(z*z-1)).
1053 return log($z + sqrt($z*$z-1)) if $z >= 1;
1056 my ($re, $im) = @{$z->cartesian};
1058 return cplx(log($re + sqrt($re*$re - 1)), 0) if $re >= 1;
1059 return cplx(0, atan2(sqrt(1-$re*$re), $re)) if abs($re) <= 1;
1061 return log($z + sqrt($z*$z - 1));
1067 # Computes the arc hyperbolic sine asinh(z) = log(z + sqrt(z*z-1))
1071 return log($z + sqrt($z*$z + 1));
1077 # Computes the arc hyperbolic tangent atanh(z) = 1/2 log((1+z) / (1-z)).
1082 return log((1 + $z)/(1 - $z))/2 if abs($z) < 1;
1085 _divbyzero 'atanh(1)', "1 - $z" if ($z == 1);
1086 _logofzero 'atanh(-1)' if ($z == -1);
1087 return 0.5 * log((1 + $z) / (1 - $z));
1093 # Computes the hyperbolic arc secant asech(z) = acosh(1 / z).
1097 _divbyzero 'asech(0)', $z if ($z == 0);
1098 return acosh(1 / $z);
1104 # Computes the hyperbolic arc cosecant acsch(z) = asinh(1 / z).
1108 _divbyzero 'acsch(0)', $z if ($z == 0);
1109 return asinh(1 / $z);
1115 # Alias for acosh().
1117 sub acosech { Math::Complex::acsch(@_) }
1122 # Computes the arc hyperbolic cotangent acoth(z) = 1/2 log((1+z) / (z-1)).
1126 _divbyzero 'acoth(0)' if (abs($z) < $eps);
1128 return log(($z + 1)/($z - 1))/2 if abs($z) > 1;
1131 _divbyzero 'acoth(1)', "$z - 1" if (abs($z - 1) < $eps);
1132 _logofzero 'acoth(-1)', "1 / $z" if (abs($z + 1) < $eps);
1133 return log((1 + $z) / ($z - 1)) / 2;
1141 sub acotanh { Math::Complex::acoth(@_) }
1146 # Compute atan(z1/z2).
1149 my ($z1, $z2, $inverted) = @_;
1150 my ($re1, $im1, $re2, $im2);
1152 ($re1, $im1) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
1153 ($re2, $im2) = @{$z1->cartesian};
1155 ($re1, $im1) = @{$z1->cartesian};
1156 ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
1159 return cplx(atan2($re1, $re2), 0) if $im1 == 0;
1160 return cplx(($im1<=>0) * pip2, 0) if $re2 == 0;
1162 my $w = atan($z1/$z2);
1163 my ($u, $v) = ref $w ? @{$w->cartesian} : ($w, 0);
1164 $u += pi if $re2 < 0;
1165 $u -= pit2 if $u > pi;
1166 return cplx($u, $v);
1173 # Set (fetch if no argument) display format for all complex numbers that
1174 # don't happen to have overridden it via ->display_format
1176 # When called as a method, this actually sets the display format for
1177 # the current object.
1179 # Valid object formats are 'c' and 'p' for cartesian and polar. The first
1180 # letter is used actually, so the type can be fully spelled out for clarity.
1182 sub display_format {
1186 if (ref $self) { # Called as a method
1188 } else { # Regular procedure call
1193 if (defined $self) {
1194 return defined $self->{display} ? $self->{display} : $display
1195 unless defined $format;
1196 return $self->{display} = $format;
1199 return $display unless defined $format;
1200 return $display = $format;
1206 # Show nicely formatted complex number under its cartesian or polar form,
1207 # depending on the current display format:
1209 # . If a specific display format has been recorded for this object, use it.
1210 # . Otherwise, use the generic current default for all complex numbers,
1211 # which is a package global variable.
1218 $format = $z->{display} if defined $z->{display};
1220 return $z->stringify_polar if $format =~ /^p/i;
1221 return $z->stringify_cartesian;
1225 # ->stringify_cartesian
1227 # Stringify as a cartesian representation 'a+bi'.
1229 sub stringify_cartesian {
1231 my ($x, $y) = @{$z->cartesian};
1234 $x = int($x + ($x < 0 ? -1 : 1) * $eps)
1235 if int(abs($x)) != int(abs($x) + $eps);
1236 $y = int($y + ($y < 0 ? -1 : 1) * $eps)
1237 if int(abs($y)) != int(abs($y) + $eps);
1239 $re = "$x" if abs($x) >= $eps;
1240 if ($y == 1) { $im = 'i' }
1241 elsif ($y == -1) { $im = '-i' }
1242 elsif (abs($y) >= $eps) { $im = $y . "i" }
1245 $str = $re if defined $re;
1246 $str .= "+$im" if defined $im;
1249 $str =~ s/([-+])1i/$1i/; # Not redundant with the above 1/-1 tests.
1250 $str = '0' unless $str;
1256 # Helper for stringify_polar, a Greatest Common Divisor with a memory.
1263 # Loops forever if given negative inputs.
1265 if ($b and $a > $b) { return gcd($a % $b, $b) }
1266 elsif ($a and $b > $a) { return gcd($b % $a, $a) }
1267 else { return $a ? $a : $b }
1277 unless (exists $gcd{$id}) {
1278 $gcd{$id} = _gcd($a, $b);
1279 $gcd{"$b $a"} = $gcd{$id};
1288 # Stringify as a polar representation '[r,t]'.
1290 sub stringify_polar {
1292 my ($r, $t) = @{$z->polar};
1295 return '[0,0]' if $r <= $eps;
1298 $nt = ($nt - int($nt)) * pit2;
1299 $nt += pit2 if $nt < 0; # Range [0, 2pi]
1301 if (abs($nt) <= $eps) { $theta = 0 }
1302 elsif (abs(pi-$nt) <= $eps) { $theta = 'pi' }
1304 if (defined $theta) {
1305 $r = int($r + ($r < 0 ? -1 : 1) * $eps)
1306 if int(abs($r)) != int(abs($r) + $eps);
1307 $theta = int($theta + ($theta < 0 ? -1 : 1) * $eps)
1308 if ($theta ne 'pi' and
1309 int(abs($theta)) != int(abs($theta) + $eps));
1310 return "\[$r,$theta\]";
1314 # Okay, number is not a real. Try to identify pi/n and friends...
1317 $nt -= pit2 if $nt > pi;
1319 if (abs($nt) >= deg1) {
1322 for ($k = 1, $kpi = pi; $k < 10; $k++, $kpi += pi) {
1323 $n = int($kpi / $nt + ($nt > 0 ? 1 : -1) * 0.5);
1324 if (abs($kpi/$n - $nt) <= $eps) {
1326 my $gcd = gcd($k, $n);
1332 $theta = ($nt < 0 ? '-':'').
1333 ($k == 1 ? 'pi':"${k}pi");
1334 $theta .= '/'.$n if $n > 1;
1340 $theta = $nt unless defined $theta;
1342 $r = int($r + ($r < 0 ? -1 : 1) * $eps)
1343 if int(abs($r)) != int(abs($r) + $eps);
1344 $theta = int($theta + ($theta < 0 ? -1 : 1) * $eps)
1345 if ($theta !~ m(^-?\d*pi/\d+$) and
1346 int(abs($theta)) != int(abs($theta) + $eps));
1348 return "\[$r,$theta\]";
1356 Math::Complex - complex numbers and associated mathematical functions
1362 $z = Math::Complex->make(5, 6);
1364 $j = cplxe(1, 2*pi/3);
1368 This package lets you create and manipulate complex numbers. By default,
1369 I<Perl> limits itself to real numbers, but an extra C<use> statement brings
1370 full complex support, along with a full set of mathematical functions
1371 typically associated with and/or extended to complex numbers.
1373 If you wonder what complex numbers are, they were invented to be able to solve
1374 the following equation:
1378 and by definition, the solution is noted I<i> (engineers use I<j> instead since
1379 I<i> usually denotes an intensity, but the name does not matter). The number
1380 I<i> is a pure I<imaginary> number.
1382 The arithmetics with pure imaginary numbers works just like you would expect
1383 it with real numbers... you just have to remember that
1389 5i + 7i = i * (5 + 7) = 12i
1390 4i - 3i = i * (4 - 3) = i
1395 Complex numbers are numbers that have both a real part and an imaginary
1396 part, and are usually noted:
1400 where C<a> is the I<real> part and C<b> is the I<imaginary> part. The
1401 arithmetic with complex numbers is straightforward. You have to
1402 keep track of the real and the imaginary parts, but otherwise the
1403 rules used for real numbers just apply:
1405 (4 + 3i) + (5 - 2i) = (4 + 5) + i(3 - 2) = 9 + i
1406 (2 + i) * (4 - i) = 2*4 + 4i -2i -i*i = 8 + 2i + 1 = 9 + 2i
1408 A graphical representation of complex numbers is possible in a plane
1409 (also called the I<complex plane>, but it's really a 2D plane).
1414 is the point whose coordinates are (a, b). Actually, it would
1415 be the vector originating from (0, 0) to (a, b). It follows that the addition
1416 of two complex numbers is a vectorial addition.
1418 Since there is a bijection between a point in the 2D plane and a complex
1419 number (i.e. the mapping is unique and reciprocal), a complex number
1420 can also be uniquely identified with polar coordinates:
1424 where C<rho> is the distance to the origin, and C<theta> the angle between
1425 the vector and the I<x> axis. There is a notation for this using the
1426 exponential form, which is:
1428 rho * exp(i * theta)
1430 where I<i> is the famous imaginary number introduced above. Conversion
1431 between this form and the cartesian form C<a + bi> is immediate:
1433 a = rho * cos(theta)
1434 b = rho * sin(theta)
1436 which is also expressed by this formula:
1438 z = rho * exp(i * theta) = rho * (cos theta + i * sin theta)
1440 In other words, it's the projection of the vector onto the I<x> and I<y>
1441 axes. Mathematicians call I<rho> the I<norm> or I<modulus> and I<theta>
1442 the I<argument> of the complex number. The I<norm> of C<z> will be
1445 The polar notation (also known as the trigonometric
1446 representation) is much more handy for performing multiplications and
1447 divisions of complex numbers, whilst the cartesian notation is better
1448 suited for additions and subtractions. Real numbers are on the I<x>
1449 axis, and therefore I<theta> is zero or I<pi>.
1451 All the common operations that can be performed on a real number have
1452 been defined to work on complex numbers as well, and are merely
1453 I<extensions> of the operations defined on real numbers. This means
1454 they keep their natural meaning when there is no imaginary part, provided
1455 the number is within their definition set.
1457 For instance, the C<sqrt> routine which computes the square root of
1458 its argument is only defined for non-negative real numbers and yields a
1459 non-negative real number (it is an application from B<R+> to B<R+>).
1460 If we allow it to return a complex number, then it can be extended to
1461 negative real numbers to become an application from B<R> to B<C> (the
1462 set of complex numbers):
1464 sqrt(x) = x >= 0 ? sqrt(x) : sqrt(-x)*i
1466 It can also be extended to be an application from B<C> to B<C>,
1467 whilst its restriction to B<R> behaves as defined above by using
1468 the following definition:
1470 sqrt(z = [r,t]) = sqrt(r) * exp(i * t/2)
1472 Indeed, a negative real number can be noted C<[x,pi]> (the modulus
1473 I<x> is always non-negative, so C<[x,pi]> is really C<-x>, a negative
1474 number) and the above definition states that
1476 sqrt([x,pi]) = sqrt(x) * exp(i*pi/2) = [sqrt(x),pi/2] = sqrt(x)*i
1478 which is exactly what we had defined for negative real numbers above.
1479 The C<sqrt> returns only one of the solutions: if you want the both,
1480 use the C<root> function.
1482 All the common mathematical functions defined on real numbers that
1483 are extended to complex numbers share that same property of working
1484 I<as usual> when the imaginary part is zero (otherwise, it would not
1485 be called an extension, would it?).
1487 A I<new> operation possible on a complex number that is
1488 the identity for real numbers is called the I<conjugate>, and is noted
1489 with an horizontal bar above the number, or C<~z> here.
1496 z * ~z = (a + bi) * (a - bi) = a*a + b*b
1498 We saw that the norm of C<z> was noted C<abs(z)> and was defined as the
1499 distance to the origin, also known as:
1501 rho = abs(z) = sqrt(a*a + b*b)
1505 z * ~z = abs(z) ** 2
1507 If z is a pure real number (i.e. C<b == 0>), then the above yields:
1511 which is true (C<abs> has the regular meaning for real number, i.e. stands
1512 for the absolute value). This example explains why the norm of C<z> is
1513 noted C<abs(z)>: it extends the C<abs> function to complex numbers, yet
1514 is the regular C<abs> we know when the complex number actually has no
1515 imaginary part... This justifies I<a posteriori> our use of the C<abs>
1516 notation for the norm.
1520 Given the following notations:
1522 z1 = a + bi = r1 * exp(i * t1)
1523 z2 = c + di = r2 * exp(i * t2)
1524 z = <any complex or real number>
1526 the following (overloaded) operations are supported on complex numbers:
1528 z1 + z2 = (a + c) + i(b + d)
1529 z1 - z2 = (a - c) + i(b - d)
1530 z1 * z2 = (r1 * r2) * exp(i * (t1 + t2))
1531 z1 / z2 = (r1 / r2) * exp(i * (t1 - t2))
1532 z1 ** z2 = exp(z2 * log z1)
1534 abs(z) = r1 = sqrt(a*a + b*b)
1535 sqrt(z) = sqrt(r1) * exp(i * t/2)
1536 exp(z) = exp(a) * exp(i * b)
1537 log(z) = log(r1) + i*t
1538 sin(z) = 1/2i (exp(i * z1) - exp(-i * z))
1539 cos(z) = 1/2 (exp(i * z1) + exp(-i * z))
1540 atan2(z1, z2) = atan(z1/z2)
1542 The following extra operations are supported on both real and complex
1550 cbrt(z) = z ** (1/3)
1551 log10(z) = log(z) / log(10)
1552 logn(z, n) = log(z) / log(n)
1554 tan(z) = sin(z) / cos(z)
1560 asin(z) = -i * log(i*z + sqrt(1-z*z))
1561 acos(z) = -i * log(z + i*sqrt(1-z*z))
1562 atan(z) = i/2 * log((i+z) / (i-z))
1564 acsc(z) = asin(1 / z)
1565 asec(z) = acos(1 / z)
1566 acot(z) = atan(1 / z) = -i/2 * log((i+z) / (z-i))
1568 sinh(z) = 1/2 (exp(z) - exp(-z))
1569 cosh(z) = 1/2 (exp(z) + exp(-z))
1570 tanh(z) = sinh(z) / cosh(z) = (exp(z) - exp(-z)) / (exp(z) + exp(-z))
1572 csch(z) = 1 / sinh(z)
1573 sech(z) = 1 / cosh(z)
1574 coth(z) = 1 / tanh(z)
1576 asinh(z) = log(z + sqrt(z*z+1))
1577 acosh(z) = log(z + sqrt(z*z-1))
1578 atanh(z) = 1/2 * log((1+z) / (1-z))
1580 acsch(z) = asinh(1 / z)
1581 asech(z) = acosh(1 / z)
1582 acoth(z) = atanh(1 / z) = 1/2 * log((1+z) / (z-1))
1584 I<arg>, I<abs>, I<log>, I<csc>, I<cot>, I<acsc>, I<acot>, I<csch>,
1585 I<coth>, I<acosech>, I<acotanh>, have aliases I<rho>, I<theta>, I<ln>,
1586 I<cosec>, I<cotan>, I<acosec>, I<acotan>, I<cosech>, I<cotanh>,
1587 I<acosech>, I<acotanh>, respectively. C<Re>, C<Im>, C<arg>, C<abs>,
1588 C<rho>, and C<theta> can be used also also mutators. The C<cbrt>
1589 returns only one of the solutions: if you want all three, use the
1592 The I<root> function is available to compute all the I<n>
1593 roots of some complex, where I<n> is a strictly positive integer.
1594 There are exactly I<n> such roots, returned as a list. Getting the
1595 number mathematicians call C<j> such that:
1599 is a simple matter of writing:
1601 $j = ((root(1, 3))[1];
1603 The I<k>th root for C<z = [r,t]> is given by:
1605 (root(z, n))[k] = r**(1/n) * exp(i * (t + 2*k*pi)/n)
1607 The I<spaceship> comparison operator, E<lt>=E<gt>, is also defined. In
1608 order to ensure its restriction to real numbers is conform to what you
1609 would expect, the comparison is run on the real part of the complex
1610 number first, and imaginary parts are compared only when the real
1615 To create a complex number, use either:
1617 $z = Math::Complex->make(3, 4);
1620 if you know the cartesian form of the number, or
1624 if you like. To create a number using the polar form, use either:
1626 $z = Math::Complex->emake(5, pi/3);
1627 $x = cplxe(5, pi/3);
1629 instead. The first argument is the modulus, the second is the angle
1630 (in radians, the full circle is 2*pi). (Mnemonic: C<e> is used as a
1631 notation for complex numbers in the polar form).
1633 It is possible to write:
1635 $x = cplxe(-3, pi/4);
1637 but that will be silently converted into C<[3,-3pi/4]>, since the modulus
1638 must be non-negative (it represents the distance to the origin in the complex
1641 It is also possible to have a complex number as either argument of
1642 either the C<make> or C<emake>: the appropriate component of
1643 the argument will be used.
1648 =head1 STRINGIFICATION
1650 When printed, a complex number is usually shown under its cartesian
1651 form I<a+bi>, but there are legitimate cases where the polar format
1652 I<[r,t]> is more appropriate.
1654 By calling the routine C<Math::Complex::display_format> and supplying either
1655 C<"polar"> or C<"cartesian">, you override the default display format,
1656 which is C<"cartesian">. Not supplying any argument returns the current
1659 This default can be overridden on a per-number basis by calling the
1660 C<display_format> method instead. As before, not supplying any argument
1661 returns the current display format for this number. Otherwise whatever you
1662 specify will be the new display format for I<this> particular number.
1668 Math::Complex::display_format('polar');
1669 $j = ((root(1, 3))[1];
1670 print "j = $j\n"; # Prints "j = [1,2pi/3]
1671 $j->display_format('cartesian');
1672 print "j = $j\n"; # Prints "j = -0.5+0.866025403784439i"
1674 The polar format attempts to emphasize arguments like I<k*pi/n>
1675 (where I<n> is a positive integer and I<k> an integer within [-9,+9]).
1679 Thanks to overloading, the handling of arithmetics with complex numbers
1680 is simple and almost transparent.
1682 Here are some examples:
1686 $j = cplxe(1, 2*pi/3); # $j ** 3 == 1
1687 print "j = $j, j**3 = ", $j ** 3, "\n";
1688 print "1 + j + j**2 = ", 1 + $j + $j**2, "\n";
1690 $z = -16 + 0*i; # Force it to be a complex
1691 print "sqrt($z) = ", sqrt($z), "\n";
1693 $k = exp(i * 2*pi/3);
1694 print "$j - $k = ", $j - $k, "\n";
1696 $z->Re(3); # Re, Im, arg, abs,
1697 $j->arg(2); # (the last two aka rho, theta)
1698 # can be used also as mutators.
1700 =head1 ERRORS DUE TO DIVISION BY ZERO OR LOGARITHM OF ZERO
1702 The division (/) and the following functions
1708 atanh asech acsch acoth
1710 cannot be computed for all arguments because that would mean dividing
1711 by zero or taking logarithm of zero. These situations cause fatal
1712 runtime errors looking like this
1714 cot(0): Division by zero.
1715 (Because in the definition of cot(0), the divisor sin(0) is 0)
1720 atanh(-1): Logarithm of zero.
1723 For the C<csc>, C<cot>, C<asec>, C<acsc>, C<acot>, C<csch>, C<coth>,
1724 C<asech>, C<acsch>, the argument cannot be C<0> (zero). For the the
1725 logarithmic functions and the C<atanh>, C<acoth>, the argument cannot
1726 be C<1> (one). For the C<atanh>, C<acoth>, the argument cannot be
1727 C<-1> (minus one). For the C<atan>, C<acot>, the argument cannot be
1728 C<i> (the imaginary unit). For the C<atan>, C<acoth>, the argument
1729 cannot be C<-i> (the negative imaginary unit). For the C<tan>,
1730 C<sec>, C<tanh>, the argument cannot be I<pi/2 + k * pi>, where I<k>
1733 Note that because we are operating on approximations of real numbers,
1734 these errors can happen when merely `too close' to the singularities
1735 listed above. For example C<tan(2*atan2(1,1)+1e-15)> will die of
1738 =head1 ERRORS DUE TO INDIGESTIBLE ARGUMENTS
1740 The C<make> and C<emake> accept both real and complex arguments.
1741 When they cannot recognize the arguments they will die with error
1742 messages like the following
1744 Math::Complex::make: Cannot take real part of ...
1745 Math::Complex::make: Cannot take real part of ...
1746 Math::Complex::emake: Cannot take rho of ...
1747 Math::Complex::emake: Cannot take theta of ...
1751 Saying C<use Math::Complex;> exports many mathematical routines in the
1752 caller environment and even overrides some (C<sqrt>, C<log>).
1753 This is construed as a feature by the Authors, actually... ;-)
1755 All routines expect to be given real or complex numbers. Don't attempt to
1756 use BigFloat, since Perl has currently no rule to disambiguate a '+'
1757 operation (for instance) between two overloaded entities.
1759 In Cray UNICOS there is some strange numerical instability that results
1760 in root(), cos(), sin(), cosh(), sinh(), losing accuracy fast. Beware.
1761 The bug may be in UNICOS math libs, in UNICOS C compiler, in Math::Complex.
1762 Whatever it is, it does not manifest itself anywhere else where Perl runs.
1766 Raphael Manfredi <F<Raphael_Manfredi@grenoble.hp.com>> and
1767 Jarkko Hietaniemi <F<jhi@iki.fi>>.
1769 Extensive patches by Daniel S. Lewart <F<d-lewart@uiuc.edu>>.