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.26 1998/11/01 00:00:00 dsl 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 * CORE::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 / CORE::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 * CORE::cos($t), $r * CORE::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'} = [CORE::sqrt($x*$x + $y*$y), CORE::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 # Computes z1**z2 = exp(z2 * log z1)).
409 my ($z1, $z2, $inverted) = @_;
411 return 1 if $z1 == 0 || $z2 == 1;
412 return 0 if $z2 == 0 && Re($z1) > 0;
414 return 1 if $z2 == 0 || $z1 == 1;
415 return 0 if $z1 == 0 && Re($z2) > 0;
417 my $w = $inverted ? CORE::exp($z1 * CORE::log($z2))
418 : CORE::exp($z2 * CORE::log($z1));
419 # If both arguments cartesian, return cartesian, else polar.
420 return $z1->{c_dirty} == 0 &&
421 (not ref $z2 or $z2->{c_dirty} == 0) ?
422 cplx(@{$w->cartesian}) : $w;
428 # Computes z1 <=> z2.
429 # Sorts on the real part first, then on the imaginary part. Thus 2-4i < 3+8i.
432 my ($z1, $z2, $inverted) = @_;
433 my ($re1, $im1) = ref $z1 ? @{$z1->cartesian} : ($z1, 0);
434 my ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
435 my $sgn = $inverted ? -1 : 1;
436 return $sgn * ($re1 <=> $re2) if $re1 != $re2;
437 return $sgn * ($im1 <=> $im2);
448 my ($r, $t) = @{$z->polar};
449 $t = ($t <= 0) ? $t + pi : $t - pi;
450 return (ref $z)->emake($r, $t);
452 my ($re, $im) = @{$z->cartesian};
453 return (ref $z)->make(-$re, -$im);
459 # Compute complex's conjugate.
464 my ($r, $t) = @{$z->polar};
465 return (ref $z)->emake($r, -$t);
467 my ($re, $im) = @{$z->cartesian};
468 return (ref $z)->make($re, -$im);
474 # Compute or set complex's norm (rho).
478 return $z unless ref $z;
480 $z->{'polar'} = [ $rho, ${$z->polar}[1] ];
485 return ${$z->polar}[0];
492 if ($$theta > pi()) { $$theta -= pit2 }
493 elsif ($$theta <= -pi()) { $$theta += pit2 }
499 # Compute or set complex's argument (theta).
502 my ($z, $theta) = @_;
503 return $z unless ref $z;
504 if (defined $theta) {
506 $z->{'polar'} = [ ${$z->polar}[0], $theta ];
510 $theta = ${$z->polar}[1];
521 # It is quite tempting to use wantarray here so that in list context
522 # sqrt() would return the two solutions. This, however, would
525 # print "sqrt(z) = ", sqrt($z), "\n";
527 # The two values would be printed side by side without no intervening
528 # whitespace, quite confusing.
529 # Therefore if you want the two solutions use the root().
533 my ($re, $im) = ref $z ? @{$z->cartesian} : ($z, 0);
534 return $re < 0 ? cplx(0, CORE::sqrt(-$re)) : CORE::sqrt($re) if $im == 0;
535 my ($r, $t) = @{$z->polar};
536 return (ref $z)->emake(CORE::sqrt($r), $t/2);
542 # Compute cbrt(z) (cubic root).
544 # Why are we not returning three values? The same answer as for sqrt().
548 return $z < 0 ? -CORE::exp(CORE::log(-$z)/3) : ($z > 0 ? CORE::exp(CORE::log($z)/3): 0)
550 my ($r, $t) = @{$z->polar};
551 return (ref $z)->emake(CORE::exp(CORE::log($r)/3), $t/3);
560 my $mess = "Root $_[0] not defined, root must be positive integer.\n";
564 $mess .= "Died at $up[1] line $up[2].\n";
572 # Computes all nth root for z, returning an array whose size is n.
573 # `n' must be a positive integer.
575 # The roots are given by (for k = 0..n-1):
577 # z^(1/n) = r^(1/n) (cos ((t+2 k pi)/n) + i sin ((t+2 k pi)/n))
581 _rootbad($n) if ($n < 1 or int($n) != $n);
582 my ($r, $t) = ref $z ? @{$z->polar} : (CORE::abs($z), $z >= 0 ? 0 : pi);
585 my $theta_inc = pit2 / $n;
586 my $rho = $r ** (1/$n);
588 my $cartesian = ref $z && $z->{c_dirty} == 0;
589 for ($k = 0, $theta = $t / $n; $k < $n; $k++, $theta += $theta_inc) {
590 my $w = cplxe($rho, $theta);
591 # Yes, $cartesian is loop invariant.
592 push @root, $cartesian ? cplx(@{$w->cartesian}) : $w;
600 # Return or set Re(z).
604 return $z unless ref $z;
606 $z->{'cartesian'} = [ $Re, ${$z->cartesian}[1] ];
610 return ${$z->cartesian}[0];
617 # Return or set Im(z).
621 return $z unless ref $z;
623 $z->{'cartesian'} = [ ${$z->cartesian}[0], $Im ];
627 return ${$z->cartesian}[1];
634 # Return or set rho(w).
637 Math::Complex::abs(@_);
643 # Return or set theta(w).
646 Math::Complex::arg(@_);
656 my ($x, $y) = @{$z->cartesian};
657 return (ref $z)->emake(CORE::exp($x), $y);
663 # Die on logarithm of zero.
666 my $mess = "$_[0]: Logarithm of zero.\n";
669 $mess .= "(Because in the definition of $_[0], the argument ";
670 $mess .= "$_[1] " unless ($_[1] eq '0');
676 $mess .= "Died at $up[1] line $up[2].\n";
689 _logofzero("log") if $z == 0;
690 return $z > 0 ? CORE::log($z) : cplx(CORE::log(-$z), pi);
692 my ($r, $t) = @{$z->polar};
693 _logofzero("log") if $r == 0;
694 if ($t > pi()) { $t -= pit2 }
695 elsif ($t <= -pi()) { $t += pit2 }
696 return (ref $z)->make(CORE::log($r), $t);
704 sub ln { Math::Complex::log(@_) }
713 return Math::Complex::log($_[0]) * uplog10;
719 # Compute logn(z,n) = log(z) / log(n)
723 $z = cplx($z, 0) unless ref $z;
724 my $logn = $logn{$n};
725 $logn = $logn{$n} = CORE::log($n) unless defined $logn; # Cache log(n)
726 return CORE::log($z) / $logn;
732 # Compute cos(z) = (exp(iz) + exp(-iz))/2.
736 my ($x, $y) = @{$z->cartesian};
737 my $ey = CORE::exp($y);
739 return (ref $z)->make(CORE::cos($x) * ($ey + $ey_1)/2,
740 CORE::sin($x) * ($ey_1 - $ey)/2);
746 # Compute sin(z) = (exp(iz) - exp(-iz))/2.
750 my ($x, $y) = @{$z->cartesian};
751 my $ey = CORE::exp($y);
753 return (ref $z)->make(CORE::sin($x) * ($ey + $ey_1)/2,
754 CORE::cos($x) * ($ey - $ey_1)/2);
760 # Compute tan(z) = sin(z) / cos(z).
764 my $cz = CORE::cos($z);
765 _divbyzero "tan($z)", "cos($z)" if (CORE::abs($cz) < $eps);
766 return CORE::sin($z) / $cz;
772 # Computes the secant sec(z) = 1 / cos(z).
776 my $cz = CORE::cos($z);
777 _divbyzero "sec($z)", "cos($z)" if ($cz == 0);
784 # Computes the cosecant csc(z) = 1 / sin(z).
788 my $sz = CORE::sin($z);
789 _divbyzero "csc($z)", "sin($z)" if ($sz == 0);
798 sub cosec { Math::Complex::csc(@_) }
803 # Computes cot(z) = cos(z) / sin(z).
807 my $sz = CORE::sin($z);
808 _divbyzero "cot($z)", "sin($z)" if ($sz == 0);
809 return CORE::cos($z) / $sz;
817 sub cotan { Math::Complex::cot(@_) }
822 # Computes the arc cosine acos(z) = -i log(z + sqrt(z*z-1)).
826 return CORE::atan2(CORE::sqrt(1-$z*$z), $z) if (! ref $z) && CORE::abs($z) <= 1;
827 my ($x, $y) = ref $z ? @{$z->cartesian} : ($z, 0);
828 my $t1 = CORE::sqrt(($x+1)*($x+1) + $y*$y);
829 my $t2 = CORE::sqrt(($x-1)*($x-1) + $y*$y);
830 my $alpha = ($t1 + $t2)/2;
831 my $beta = ($t1 - $t2)/2;
832 $alpha = 1 if $alpha < 1;
833 if ($beta > 1) { $beta = 1 }
834 elsif ($beta < -1) { $beta = -1 }
835 my $u = CORE::atan2(CORE::sqrt(1-$beta*$beta), $beta);
836 my $v = CORE::log($alpha + CORE::sqrt($alpha*$alpha-1));
837 $v = -$v if $y > 0 || ($y == 0 && $x < -1);
838 return $package->make($u, $v);
844 # Computes the arc sine asin(z) = -i log(iz + sqrt(1-z*z)).
848 return CORE::atan2($z, CORE::sqrt(1-$z*$z)) if (! ref $z) && CORE::abs($z) <= 1;
849 my ($x, $y) = ref $z ? @{$z->cartesian} : ($z, 0);
850 my $t1 = CORE::sqrt(($x+1)*($x+1) + $y*$y);
851 my $t2 = CORE::sqrt(($x-1)*($x-1) + $y*$y);
852 my $alpha = ($t1 + $t2)/2;
853 my $beta = ($t1 - $t2)/2;
854 $alpha = 1 if $alpha < 1;
855 if ($beta > 1) { $beta = 1 }
856 elsif ($beta < -1) { $beta = -1 }
857 my $u = CORE::atan2($beta, CORE::sqrt(1-$beta*$beta));
858 my $v = -CORE::log($alpha + CORE::sqrt($alpha*$alpha-1));
859 $v = -$v if $y > 0 || ($y == 0 && $x < -1);
860 return $package->make($u, $v);
866 # Computes the arc tangent atan(z) = i/2 log((i+z) / (i-z)).
870 return CORE::atan2($z, 1) unless ref $z;
871 _divbyzero "atan(i)" if ( $z == i);
872 _divbyzero "atan(-i)" if (-$z == i);
873 my $log = CORE::log((i + $z) / (i - $z));
874 $ip2 = 0.5 * i unless defined $ip2;
881 # Computes the arc secant asec(z) = acos(1 / z).
885 _divbyzero "asec($z)", $z if ($z == 0);
892 # Computes the arc cosecant acsc(z) = asin(1 / z).
896 _divbyzero "acsc($z)", $z if ($z == 0);
905 sub acosec { Math::Complex::acsc(@_) }
910 # Computes the arc cotangent acot(z) = atan(1 / z)
914 _divbyzero "acot(0)" if (CORE::abs($z) < $eps);
915 return ($z >= 0) ? CORE::atan2(1, $z) : CORE::atan2(-1, -$z) unless ref $z;
916 _divbyzero "acot(i)" if (CORE::abs($z - i) < $eps);
917 _logofzero "acot(-i)" if (CORE::abs($z + i) < $eps);
926 sub acotan { Math::Complex::acot(@_) }
931 # Computes the hyperbolic cosine cosh(z) = (exp(z) + exp(-z))/2.
938 return ($ex + 1/$ex)/2;
940 my ($x, $y) = @{$z->cartesian};
943 return (ref $z)->make(CORE::cos($y) * ($ex + $ex_1)/2,
944 CORE::sin($y) * ($ex - $ex_1)/2);
950 # Computes the hyperbolic sine sinh(z) = (exp(z) - exp(-z))/2.
957 return ($ex - 1/$ex)/2;
959 my ($x, $y) = @{$z->cartesian};
962 return (ref $z)->make(CORE::cos($y) * ($ex - $ex_1)/2,
963 CORE::sin($y) * ($ex + $ex_1)/2);
969 # Computes the hyperbolic tangent tanh(z) = sinh(z) / cosh(z).
974 _divbyzero "tanh($z)", "cosh($z)" if ($cz == 0);
975 return sinh($z) / $cz;
981 # Computes the hyperbolic secant sech(z) = 1 / cosh(z).
986 _divbyzero "sech($z)", "cosh($z)" if ($cz == 0);
993 # Computes the hyperbolic cosecant csch(z) = 1 / sinh(z).
998 _divbyzero "csch($z)", "sinh($z)" if ($sz == 0);
1007 sub cosech { Math::Complex::csch(@_) }
1012 # Computes the hyperbolic cotangent coth(z) = cosh(z) / sinh(z).
1017 _divbyzero "coth($z)", "sinh($z)" if ($sz == 0);
1018 return cosh($z) / $sz;
1026 sub cotanh { Math::Complex::coth(@_) }
1031 # Computes the arc hyperbolic cosine acosh(z) = log(z + sqrt(z*z-1)).
1036 return CORE::log($z + CORE::sqrt($z*$z-1)) if $z >= 1;
1039 my ($re, $im) = @{$z->cartesian};
1041 return cplx(CORE::log($re + CORE::sqrt($re*$re - 1)), 0) if $re >= 1;
1042 return cplx(0, CORE::atan2(CORE::sqrt(1-$re*$re), $re)) if CORE::abs($re) <= 1;
1044 return CORE::log($z + CORE::sqrt($z*$z - 1));
1050 # Computes the arc hyperbolic sine asinh(z) = log(z + sqrt(z*z-1))
1054 return CORE::log($z + CORE::sqrt($z*$z + 1));
1060 # Computes the arc hyperbolic tangent atanh(z) = 1/2 log((1+z) / (1-z)).
1065 return CORE::log((1 + $z)/(1 - $z))/2 if CORE::abs($z) < 1;
1068 _divbyzero 'atanh(1)', "1 - $z" if ($z == 1);
1069 _logofzero 'atanh(-1)' if ($z == -1);
1070 return 0.5 * CORE::log((1 + $z) / (1 - $z));
1076 # Computes the hyperbolic arc secant asech(z) = acosh(1 / z).
1080 _divbyzero 'asech(0)', $z if ($z == 0);
1081 return acosh(1 / $z);
1087 # Computes the hyperbolic arc cosecant acsch(z) = asinh(1 / z).
1091 _divbyzero 'acsch(0)', $z if ($z == 0);
1092 return asinh(1 / $z);
1098 # Alias for acosh().
1100 sub acosech { Math::Complex::acsch(@_) }
1105 # Computes the arc hyperbolic cotangent acoth(z) = 1/2 log((1+z) / (z-1)).
1109 _divbyzero 'acoth(0)' if (CORE::abs($z) < $eps);
1111 return CORE::log(($z + 1)/($z - 1))/2 if CORE::abs($z) > 1;
1114 _divbyzero 'acoth(1)', "$z - 1" if (CORE::abs($z - 1) < $eps);
1115 _logofzero 'acoth(-1)', "1 / $z" if (CORE::abs($z + 1) < $eps);
1116 return CORE::log((1 + $z) / ($z - 1)) / 2;
1124 sub acotanh { Math::Complex::acoth(@_) }
1129 # Compute atan(z1/z2).
1132 my ($z1, $z2, $inverted) = @_;
1133 my ($re1, $im1, $re2, $im2);
1135 ($re1, $im1) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
1136 ($re2, $im2) = @{$z1->cartesian};
1138 ($re1, $im1) = @{$z1->cartesian};
1139 ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
1142 return cplx(CORE::atan2($re1, $re2), 0) if $im1 == 0;
1143 return cplx(($im1<=>0) * pip2, 0) if $re2 == 0;
1145 my $w = atan($z1/$z2);
1146 my ($u, $v) = ref $w ? @{$w->cartesian} : ($w, 0);
1147 $u += pi if $re2 < 0;
1148 $u -= pit2 if $u > pi;
1149 return cplx($u, $v);
1156 # Set (fetch if no argument) display format for all complex numbers that
1157 # don't happen to have overridden it via ->display_format
1159 # When called as a method, this actually sets the display format for
1160 # the current object.
1162 # Valid object formats are 'c' and 'p' for cartesian and polar. The first
1163 # letter is used actually, so the type can be fully spelled out for clarity.
1165 sub display_format {
1169 if (ref $self) { # Called as a method
1171 } else { # Regular procedure call
1176 if (defined $self) {
1177 return defined $self->{display} ? $self->{display} : $display
1178 unless defined $format;
1179 return $self->{display} = $format;
1182 return $display unless defined $format;
1183 return $display = $format;
1189 # Show nicely formatted complex number under its cartesian or polar form,
1190 # depending on the current display format:
1192 # . If a specific display format has been recorded for this object, use it.
1193 # . Otherwise, use the generic current default for all complex numbers,
1194 # which is a package global variable.
1201 $format = $z->{display} if defined $z->{display};
1203 return $z->stringify_polar if $format =~ /^p/i;
1204 return $z->stringify_cartesian;
1208 # ->stringify_cartesian
1210 # Stringify as a cartesian representation 'a+bi'.
1212 sub stringify_cartesian {
1214 my ($x, $y) = @{$z->cartesian};
1217 $x = int($x + ($x < 0 ? -1 : 1) * $eps)
1218 if int(CORE::abs($x)) != int(CORE::abs($x) + $eps);
1219 $y = int($y + ($y < 0 ? -1 : 1) * $eps)
1220 if int(CORE::abs($y)) != int(CORE::abs($y) + $eps);
1222 $re = "$x" if CORE::abs($x) >= $eps;
1223 if ($y == 1) { $im = 'i' }
1224 elsif ($y == -1) { $im = '-i' }
1225 elsif (CORE::abs($y) >= $eps) { $im = $y . "i" }
1228 $str = $re if defined $re;
1229 $str .= "+$im" if defined $im;
1232 $str =~ s/([-+])1i/$1i/; # Not redundant with the above 1/-1 tests.
1233 $str = '0' unless $str;
1239 # Helper for stringify_polar, a Greatest Common Divisor with a memory.
1246 # Loops forever if given negative inputs.
1248 if ($b and $a > $b) { return gcd($a % $b, $b) }
1249 elsif ($a and $b > $a) { return gcd($b % $a, $a) }
1250 else { return $a ? $a : $b }
1260 unless (exists $gcd{$id}) {
1261 $gcd{$id} = _gcd($a, $b);
1262 $gcd{"$b $a"} = $gcd{$id};
1271 # Stringify as a polar representation '[r,t]'.
1273 sub stringify_polar {
1275 my ($r, $t) = @{$z->polar};
1278 return '[0,0]' if $r <= $eps;
1281 $nt = ($nt - int($nt)) * pit2;
1282 $nt += pit2 if $nt < 0; # Range [0, 2pi]
1284 if (CORE::abs($nt) <= $eps) { $theta = 0 }
1285 elsif (CORE::abs(pi-$nt) <= $eps) { $theta = 'pi' }
1287 if (defined $theta) {
1288 $r = int($r + ($r < 0 ? -1 : 1) * $eps)
1289 if int(CORE::abs($r)) != int(CORE::abs($r) + $eps);
1290 $theta = int($theta + ($theta < 0 ? -1 : 1) * $eps)
1291 if ($theta ne 'pi' and
1292 int(CORE::abs($theta)) != int(CORE::abs($theta) + $eps));
1293 return "\[$r,$theta\]";
1297 # Okay, number is not a real. Try to identify pi/n and friends...
1300 $nt -= pit2 if $nt > pi;
1302 if (CORE::abs($nt) >= deg1) {
1305 for ($k = 1, $kpi = pi; $k < 10; $k++, $kpi += pi) {
1306 $n = int($kpi / $nt + ($nt > 0 ? 1 : -1) * 0.5);
1307 if (CORE::abs($kpi/$n - $nt) <= $eps) {
1309 my $gcd = gcd($k, $n);
1315 $theta = ($nt < 0 ? '-':'').
1316 ($k == 1 ? 'pi':"${k}pi");
1317 $theta .= '/'.$n if $n > 1;
1323 $theta = $nt unless defined $theta;
1325 $r = int($r + ($r < 0 ? -1 : 1) * $eps)
1326 if int(CORE::abs($r)) != int(CORE::abs($r) + $eps);
1327 $theta = int($theta + ($theta < 0 ? -1 : 1) * $eps)
1328 if ($theta !~ m(^-?\d*pi/\d+$) and
1329 int(CORE::abs($theta)) != int(CORE::abs($theta) + $eps));
1331 return "\[$r,$theta\]";
1339 Math::Complex - complex numbers and associated mathematical functions
1345 $z = Math::Complex->make(5, 6);
1347 $j = cplxe(1, 2*pi/3);
1351 This package lets you create and manipulate complex numbers. By default,
1352 I<Perl> limits itself to real numbers, but an extra C<use> statement brings
1353 full complex support, along with a full set of mathematical functions
1354 typically associated with and/or extended to complex numbers.
1356 If you wonder what complex numbers are, they were invented to be able to solve
1357 the following equation:
1361 and by definition, the solution is noted I<i> (engineers use I<j> instead since
1362 I<i> usually denotes an intensity, but the name does not matter). The number
1363 I<i> is a pure I<imaginary> number.
1365 The arithmetics with pure imaginary numbers works just like you would expect
1366 it with real numbers... you just have to remember that
1372 5i + 7i = i * (5 + 7) = 12i
1373 4i - 3i = i * (4 - 3) = i
1378 Complex numbers are numbers that have both a real part and an imaginary
1379 part, and are usually noted:
1383 where C<a> is the I<real> part and C<b> is the I<imaginary> part. The
1384 arithmetic with complex numbers is straightforward. You have to
1385 keep track of the real and the imaginary parts, but otherwise the
1386 rules used for real numbers just apply:
1388 (4 + 3i) + (5 - 2i) = (4 + 5) + i(3 - 2) = 9 + i
1389 (2 + i) * (4 - i) = 2*4 + 4i -2i -i*i = 8 + 2i + 1 = 9 + 2i
1391 A graphical representation of complex numbers is possible in a plane
1392 (also called the I<complex plane>, but it's really a 2D plane).
1397 is the point whose coordinates are (a, b). Actually, it would
1398 be the vector originating from (0, 0) to (a, b). It follows that the addition
1399 of two complex numbers is a vectorial addition.
1401 Since there is a bijection between a point in the 2D plane and a complex
1402 number (i.e. the mapping is unique and reciprocal), a complex number
1403 can also be uniquely identified with polar coordinates:
1407 where C<rho> is the distance to the origin, and C<theta> the angle between
1408 the vector and the I<x> axis. There is a notation for this using the
1409 exponential form, which is:
1411 rho * exp(i * theta)
1413 where I<i> is the famous imaginary number introduced above. Conversion
1414 between this form and the cartesian form C<a + bi> is immediate:
1416 a = rho * cos(theta)
1417 b = rho * sin(theta)
1419 which is also expressed by this formula:
1421 z = rho * exp(i * theta) = rho * (cos theta + i * sin theta)
1423 In other words, it's the projection of the vector onto the I<x> and I<y>
1424 axes. Mathematicians call I<rho> the I<norm> or I<modulus> and I<theta>
1425 the I<argument> of the complex number. The I<norm> of C<z> will be
1428 The polar notation (also known as the trigonometric
1429 representation) is much more handy for performing multiplications and
1430 divisions of complex numbers, whilst the cartesian notation is better
1431 suited for additions and subtractions. Real numbers are on the I<x>
1432 axis, and therefore I<theta> is zero or I<pi>.
1434 All the common operations that can be performed on a real number have
1435 been defined to work on complex numbers as well, and are merely
1436 I<extensions> of the operations defined on real numbers. This means
1437 they keep their natural meaning when there is no imaginary part, provided
1438 the number is within their definition set.
1440 For instance, the C<sqrt> routine which computes the square root of
1441 its argument is only defined for non-negative real numbers and yields a
1442 non-negative real number (it is an application from B<R+> to B<R+>).
1443 If we allow it to return a complex number, then it can be extended to
1444 negative real numbers to become an application from B<R> to B<C> (the
1445 set of complex numbers):
1447 sqrt(x) = x >= 0 ? sqrt(x) : sqrt(-x)*i
1449 It can also be extended to be an application from B<C> to B<C>,
1450 whilst its restriction to B<R> behaves as defined above by using
1451 the following definition:
1453 sqrt(z = [r,t]) = sqrt(r) * exp(i * t/2)
1455 Indeed, a negative real number can be noted C<[x,pi]> (the modulus
1456 I<x> is always non-negative, so C<[x,pi]> is really C<-x>, a negative
1457 number) and the above definition states that
1459 sqrt([x,pi]) = sqrt(x) * exp(i*pi/2) = [sqrt(x),pi/2] = sqrt(x)*i
1461 which is exactly what we had defined for negative real numbers above.
1462 The C<sqrt> returns only one of the solutions: if you want the both,
1463 use the C<root> function.
1465 All the common mathematical functions defined on real numbers that
1466 are extended to complex numbers share that same property of working
1467 I<as usual> when the imaginary part is zero (otherwise, it would not
1468 be called an extension, would it?).
1470 A I<new> operation possible on a complex number that is
1471 the identity for real numbers is called the I<conjugate>, and is noted
1472 with an horizontal bar above the number, or C<~z> here.
1479 z * ~z = (a + bi) * (a - bi) = a*a + b*b
1481 We saw that the norm of C<z> was noted C<abs(z)> and was defined as the
1482 distance to the origin, also known as:
1484 rho = abs(z) = sqrt(a*a + b*b)
1488 z * ~z = abs(z) ** 2
1490 If z is a pure real number (i.e. C<b == 0>), then the above yields:
1494 which is true (C<abs> has the regular meaning for real number, i.e. stands
1495 for the absolute value). This example explains why the norm of C<z> is
1496 noted C<abs(z)>: it extends the C<abs> function to complex numbers, yet
1497 is the regular C<abs> we know when the complex number actually has no
1498 imaginary part... This justifies I<a posteriori> our use of the C<abs>
1499 notation for the norm.
1503 Given the following notations:
1505 z1 = a + bi = r1 * exp(i * t1)
1506 z2 = c + di = r2 * exp(i * t2)
1507 z = <any complex or real number>
1509 the following (overloaded) operations are supported on complex numbers:
1511 z1 + z2 = (a + c) + i(b + d)
1512 z1 - z2 = (a - c) + i(b - d)
1513 z1 * z2 = (r1 * r2) * exp(i * (t1 + t2))
1514 z1 / z2 = (r1 / r2) * exp(i * (t1 - t2))
1515 z1 ** z2 = exp(z2 * log z1)
1517 abs(z) = r1 = sqrt(a*a + b*b)
1518 sqrt(z) = sqrt(r1) * exp(i * t/2)
1519 exp(z) = exp(a) * exp(i * b)
1520 log(z) = log(r1) + i*t
1521 sin(z) = 1/2i (exp(i * z1) - exp(-i * z))
1522 cos(z) = 1/2 (exp(i * z1) + exp(-i * z))
1523 atan2(z1, z2) = atan(z1/z2)
1525 The following extra operations are supported on both real and complex
1533 cbrt(z) = z ** (1/3)
1534 log10(z) = log(z) / log(10)
1535 logn(z, n) = log(z) / log(n)
1537 tan(z) = sin(z) / cos(z)
1543 asin(z) = -i * log(i*z + sqrt(1-z*z))
1544 acos(z) = -i * log(z + i*sqrt(1-z*z))
1545 atan(z) = i/2 * log((i+z) / (i-z))
1547 acsc(z) = asin(1 / z)
1548 asec(z) = acos(1 / z)
1549 acot(z) = atan(1 / z) = -i/2 * log((i+z) / (z-i))
1551 sinh(z) = 1/2 (exp(z) - exp(-z))
1552 cosh(z) = 1/2 (exp(z) + exp(-z))
1553 tanh(z) = sinh(z) / cosh(z) = (exp(z) - exp(-z)) / (exp(z) + exp(-z))
1555 csch(z) = 1 / sinh(z)
1556 sech(z) = 1 / cosh(z)
1557 coth(z) = 1 / tanh(z)
1559 asinh(z) = log(z + sqrt(z*z+1))
1560 acosh(z) = log(z + sqrt(z*z-1))
1561 atanh(z) = 1/2 * log((1+z) / (1-z))
1563 acsch(z) = asinh(1 / z)
1564 asech(z) = acosh(1 / z)
1565 acoth(z) = atanh(1 / z) = 1/2 * log((1+z) / (z-1))
1567 I<arg>, I<abs>, I<log>, I<csc>, I<cot>, I<acsc>, I<acot>, I<csch>,
1568 I<coth>, I<acosech>, I<acotanh>, have aliases I<rho>, I<theta>, I<ln>,
1569 I<cosec>, I<cotan>, I<acosec>, I<acotan>, I<cosech>, I<cotanh>,
1570 I<acosech>, I<acotanh>, respectively. C<Re>, C<Im>, C<arg>, C<abs>,
1571 C<rho>, and C<theta> can be used also also mutators. The C<cbrt>
1572 returns only one of the solutions: if you want all three, use the
1575 The I<root> function is available to compute all the I<n>
1576 roots of some complex, where I<n> is a strictly positive integer.
1577 There are exactly I<n> such roots, returned as a list. Getting the
1578 number mathematicians call C<j> such that:
1582 is a simple matter of writing:
1584 $j = ((root(1, 3))[1];
1586 The I<k>th root for C<z = [r,t]> is given by:
1588 (root(z, n))[k] = r**(1/n) * exp(i * (t + 2*k*pi)/n)
1590 The I<spaceship> comparison operator, E<lt>=E<gt>, is also defined. In
1591 order to ensure its restriction to real numbers is conform to what you
1592 would expect, the comparison is run on the real part of the complex
1593 number first, and imaginary parts are compared only when the real
1598 To create a complex number, use either:
1600 $z = Math::Complex->make(3, 4);
1603 if you know the cartesian form of the number, or
1607 if you like. To create a number using the polar form, use either:
1609 $z = Math::Complex->emake(5, pi/3);
1610 $x = cplxe(5, pi/3);
1612 instead. The first argument is the modulus, the second is the angle
1613 (in radians, the full circle is 2*pi). (Mnemonic: C<e> is used as a
1614 notation for complex numbers in the polar form).
1616 It is possible to write:
1618 $x = cplxe(-3, pi/4);
1620 but that will be silently converted into C<[3,-3pi/4]>, since the modulus
1621 must be non-negative (it represents the distance to the origin in the complex
1624 It is also possible to have a complex number as either argument of
1625 either the C<make> or C<emake>: the appropriate component of
1626 the argument will be used.
1631 =head1 STRINGIFICATION
1633 When printed, a complex number is usually shown under its cartesian
1634 form I<a+bi>, but there are legitimate cases where the polar format
1635 I<[r,t]> is more appropriate.
1637 By calling the routine C<Math::Complex::display_format> and supplying either
1638 C<"polar"> or C<"cartesian">, you override the default display format,
1639 which is C<"cartesian">. Not supplying any argument returns the current
1642 This default can be overridden on a per-number basis by calling the
1643 C<display_format> method instead. As before, not supplying any argument
1644 returns the current display format for this number. Otherwise whatever you
1645 specify will be the new display format for I<this> particular number.
1651 Math::Complex::display_format('polar');
1652 $j = ((root(1, 3))[1];
1653 print "j = $j\n"; # Prints "j = [1,2pi/3]
1654 $j->display_format('cartesian');
1655 print "j = $j\n"; # Prints "j = -0.5+0.866025403784439i"
1657 The polar format attempts to emphasize arguments like I<k*pi/n>
1658 (where I<n> is a positive integer and I<k> an integer within [-9,+9]).
1662 Thanks to overloading, the handling of arithmetics with complex numbers
1663 is simple and almost transparent.
1665 Here are some examples:
1669 $j = cplxe(1, 2*pi/3); # $j ** 3 == 1
1670 print "j = $j, j**3 = ", $j ** 3, "\n";
1671 print "1 + j + j**2 = ", 1 + $j + $j**2, "\n";
1673 $z = -16 + 0*i; # Force it to be a complex
1674 print "sqrt($z) = ", sqrt($z), "\n";
1676 $k = exp(i * 2*pi/3);
1677 print "$j - $k = ", $j - $k, "\n";
1679 $z->Re(3); # Re, Im, arg, abs,
1680 $j->arg(2); # (the last two aka rho, theta)
1681 # can be used also as mutators.
1683 =head1 ERRORS DUE TO DIVISION BY ZERO OR LOGARITHM OF ZERO
1685 The division (/) and the following functions
1691 atanh asech acsch acoth
1693 cannot be computed for all arguments because that would mean dividing
1694 by zero or taking logarithm of zero. These situations cause fatal
1695 runtime errors looking like this
1697 cot(0): Division by zero.
1698 (Because in the definition of cot(0), the divisor sin(0) is 0)
1703 atanh(-1): Logarithm of zero.
1706 For the C<csc>, C<cot>, C<asec>, C<acsc>, C<acot>, C<csch>, C<coth>,
1707 C<asech>, C<acsch>, the argument cannot be C<0> (zero). For the the
1708 logarithmic functions and the C<atanh>, C<acoth>, the argument cannot
1709 be C<1> (one). For the C<atanh>, C<acoth>, the argument cannot be
1710 C<-1> (minus one). For the C<atan>, C<acot>, the argument cannot be
1711 C<i> (the imaginary unit). For the C<atan>, C<acoth>, the argument
1712 cannot be C<-i> (the negative imaginary unit). For the C<tan>,
1713 C<sec>, C<tanh>, the argument cannot be I<pi/2 + k * pi>, where I<k>
1716 Note that because we are operating on approximations of real numbers,
1717 these errors can happen when merely `too close' to the singularities
1718 listed above. For example C<tan(2*atan2(1,1)+1e-15)> will die of
1721 =head1 ERRORS DUE TO INDIGESTIBLE ARGUMENTS
1723 The C<make> and C<emake> accept both real and complex arguments.
1724 When they cannot recognize the arguments they will die with error
1725 messages like the following
1727 Math::Complex::make: Cannot take real part of ...
1728 Math::Complex::make: Cannot take real part of ...
1729 Math::Complex::emake: Cannot take rho of ...
1730 Math::Complex::emake: Cannot take theta of ...
1734 Saying C<use Math::Complex;> exports many mathematical routines in the
1735 caller environment and even overrides some (C<sqrt>, C<log>).
1736 This is construed as a feature by the Authors, actually... ;-)
1738 All routines expect to be given real or complex numbers. Don't attempt to
1739 use BigFloat, since Perl has currently no rule to disambiguate a '+'
1740 operation (for instance) between two overloaded entities.
1742 In Cray UNICOS there is some strange numerical instability that results
1743 in root(), cos(), sin(), cosh(), sinh(), losing accuracy fast. Beware.
1744 The bug may be in UNICOS math libs, in UNICOS C compiler, in Math::Complex.
1745 Whatever it is, it does not manifest itself anywhere else where Perl runs.
1749 Raphael Manfredi <F<Raphael_Manfredi@grenoble.hp.com>> and
1750 Jarkko Hietaniemi <F<jhi@iki.fi>>.
1752 Extensive patches by Daniel S. Lewart <F<d-lewart@uiuc.edu>>.