2 # Complex numbers and associated mathematical functions
3 # -- Raphael Manfredi Since Sep 1996
4 # -- Jarkko Hietaniemi Since Mar 1997
5 # -- Daniel S. Lewart Since Sep 1997
14 our($VERSION, @ISA, @EXPORT, %EXPORT_TAGS);
16 my ( $i, $ip2, %logn );
18 $VERSION = sprintf("%s", q$Id: Complex.pm,v 1.26 1998/11/01 00:00:00 dsl Exp $ =~ /(\d+\.\d+)/);
25 csc cosec sec cot cotan
27 acsc acosec asec acot acotan
29 csch cosech sech coth cotanh
31 acsch acosech asech acoth acotanh
69 my $package = 'Math::Complex'; # Package name
70 my $display = 'cartesian'; # Default display format
71 my $eps = 1e-14; # Epsilon
74 # Object attributes (internal):
75 # cartesian [real, imaginary] -- cartesian form
76 # polar [rho, theta] -- polar form
77 # c_dirty cartesian form not up-to-date
78 # p_dirty polar form not up-to-date
79 # display display format (package's global when not set)
82 # Die on bad *make() arguments.
85 die "@{[(caller(1))[3]]}: Cannot take $_[0] of $_[1].\n";
91 # Create a new complex number (cartesian form)
94 my $self = bless {}, shift;
98 if ( $rre eq ref $self ) {
101 _cannot_make("real part", $rre);
106 if ( $rim eq ref $self ) {
109 _cannot_make("imaginary part", $rim);
112 $self->{'cartesian'} = [ $re, $im ];
113 $self->{c_dirty} = 0;
114 $self->{p_dirty} = 1;
115 $self->display_format('cartesian');
122 # Create a new complex number (exponential form)
125 my $self = bless {}, shift;
126 my ($rho, $theta) = @_;
129 if ( $rrh eq ref $self ) {
132 _cannot_make("rho", $rrh);
135 my $rth = ref $theta;
137 if ( $rth eq ref $self ) {
138 $theta = theta($theta);
140 _cannot_make("theta", $rth);
145 $theta = ($theta <= 0) ? $theta + pi() : $theta - pi();
147 $self->{'polar'} = [$rho, $theta];
148 $self->{p_dirty} = 0;
149 $self->{c_dirty} = 1;
150 $self->display_format('polar');
154 sub new { &make } # For backward compatibility only.
159 # Creates a complex number from a (re, im) tuple.
160 # This avoids the burden of writing Math::Complex->make(re, im).
164 return $package->make($re, defined $im ? $im : 0);
170 # Creates a complex number from a (rho, theta) tuple.
171 # This avoids the burden of writing Math::Complex->emake(rho, theta).
174 my ($rho, $theta) = @_;
175 return $package->emake($rho, defined $theta ? $theta : 0);
181 # The number defined as pi = 180 degrees
183 sub pi () { 4 * CORE::atan2(1, 1) }
190 sub pit2 () { 2 * pi }
197 sub pip2 () { pi / 2 }
202 # One degree in radians, used in stringify_polar.
205 sub deg1 () { pi / 180 }
212 sub uplog10 () { 1 / CORE::log(10) }
217 # The number defined as i*i = -1;
222 $i->{'cartesian'} = [0, 1];
223 $i->{'polar'} = [1, pip2];
230 # Attribute access/set routines
233 sub cartesian {$_[0]->{c_dirty} ?
234 $_[0]->update_cartesian : $_[0]->{'cartesian'}}
235 sub polar {$_[0]->{p_dirty} ?
236 $_[0]->update_polar : $_[0]->{'polar'}}
238 sub set_cartesian { $_[0]->{p_dirty}++; $_[0]->{'cartesian'} = $_[1] }
239 sub set_polar { $_[0]->{c_dirty}++; $_[0]->{'polar'} = $_[1] }
244 # Recompute and return the cartesian form, given accurate polar form.
246 sub update_cartesian {
248 my ($r, $t) = @{$self->{'polar'}};
249 $self->{c_dirty} = 0;
250 return $self->{'cartesian'} = [$r * CORE::cos($t), $r * CORE::sin($t)];
257 # Recompute and return the polar form, given accurate cartesian form.
261 my ($x, $y) = @{$self->{'cartesian'}};
262 $self->{p_dirty} = 0;
263 return $self->{'polar'} = [0, 0] if $x == 0 && $y == 0;
264 return $self->{'polar'} = [CORE::sqrt($x*$x + $y*$y), CORE::atan2($y, $x)];
273 my ($z1, $z2, $regular) = @_;
274 my ($re1, $im1) = @{$z1->cartesian};
275 $z2 = cplx($z2) unless ref $z2;
276 my ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
277 unless (defined $regular) {
278 $z1->set_cartesian([$re1 + $re2, $im1 + $im2]);
281 return (ref $z1)->make($re1 + $re2, $im1 + $im2);
290 my ($z1, $z2, $inverted) = @_;
291 my ($re1, $im1) = @{$z1->cartesian};
292 $z2 = cplx($z2) unless ref $z2;
293 my ($re2, $im2) = @{$z2->cartesian};
294 unless (defined $inverted) {
295 $z1->set_cartesian([$re1 - $re2, $im1 - $im2]);
299 (ref $z1)->make($re2 - $re1, $im2 - $im1) :
300 (ref $z1)->make($re1 - $re2, $im1 - $im2);
310 my ($z1, $z2, $regular) = @_;
311 if ($z1->{p_dirty} == 0 and ref $z2 and $z2->{p_dirty} == 0) {
312 # if both polar better use polar to avoid rounding errors
313 my ($r1, $t1) = @{$z1->polar};
314 my ($r2, $t2) = @{$z2->polar};
316 if ($t > pi()) { $t -= pit2 }
317 elsif ($t <= -pi()) { $t += pit2 }
318 unless (defined $regular) {
319 $z1->set_polar([$r1 * $r2, $t]);
322 return (ref $z1)->emake($r1 * $r2, $t);
324 my ($x1, $y1) = @{$z1->cartesian};
326 my ($x2, $y2) = @{$z2->cartesian};
327 return (ref $z1)->make($x1*$x2-$y1*$y2, $x1*$y2+$y1*$x2);
329 return (ref $z1)->make($x1*$z2, $y1*$z2);
337 # Die on division by zero.
340 my $mess = "$_[0]: Division by zero.\n";
343 $mess .= "(Because in the definition of $_[0], the divisor ";
344 $mess .= "$_[1] " unless ($_[1] eq '0');
350 $mess .= "Died at $up[1] line $up[2].\n";
361 my ($z1, $z2, $inverted) = @_;
362 if ($z1->{p_dirty} == 0 and ref $z2 and $z2->{p_dirty} == 0) {
363 # if both polar better use polar to avoid rounding errors
364 my ($r1, $t1) = @{$z1->polar};
365 my ($r2, $t2) = @{$z2->polar};
368 _divbyzero "$z2/0" if ($r1 == 0);
370 if ($t > pi()) { $t -= pit2 }
371 elsif ($t <= -pi()) { $t += pit2 }
372 return (ref $z1)->emake($r2 / $r1, $t);
374 _divbyzero "$z1/0" if ($r2 == 0);
376 if ($t > pi()) { $t -= pit2 }
377 elsif ($t <= -pi()) { $t += pit2 }
378 return (ref $z1)->emake($r1 / $r2, $t);
383 ($x2, $y2) = @{$z1->cartesian};
384 $d = $x2*$x2 + $y2*$y2;
385 _divbyzero "$z2/0" if $d == 0;
386 return (ref $z1)->make(($x2*$z2)/$d, -($y2*$z2)/$d);
388 my ($x1, $y1) = @{$z1->cartesian};
390 ($x2, $y2) = @{$z2->cartesian};
391 $d = $x2*$x2 + $y2*$y2;
392 _divbyzero "$z1/0" if $d == 0;
393 my $u = ($x1*$x2 + $y1*$y2)/$d;
394 my $v = ($y1*$x2 - $x1*$y2)/$d;
395 return (ref $z1)->make($u, $v);
397 _divbyzero "$z1/0" if $z2 == 0;
398 return (ref $z1)->make($x1/$z2, $y1/$z2);
407 # Computes z1**z2 = exp(z2 * log z1)).
410 my ($z1, $z2, $inverted) = @_;
412 return 1 if $z1 == 0 || $z2 == 1;
413 return 0 if $z2 == 0 && Re($z1) > 0;
415 return 1 if $z2 == 0 || $z1 == 1;
416 return 0 if $z1 == 0 && Re($z2) > 0;
418 my $w = $inverted ? CORE::exp($z1 * CORE::log($z2))
419 : CORE::exp($z2 * CORE::log($z1));
420 # If both arguments cartesian, return cartesian, else polar.
421 return $z1->{c_dirty} == 0 &&
422 (not ref $z2 or $z2->{c_dirty} == 0) ?
423 cplx(@{$w->cartesian}) : $w;
429 # Computes z1 <=> z2.
430 # Sorts on the real part first, then on the imaginary part. Thus 2-4i < 3+8i.
433 my ($z1, $z2, $inverted) = @_;
434 my ($re1, $im1) = ref $z1 ? @{$z1->cartesian} : ($z1, 0);
435 my ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
436 my $sgn = $inverted ? -1 : 1;
437 return $sgn * ($re1 <=> $re2) if $re1 != $re2;
438 return $sgn * ($im1 <=> $im2);
449 my ($r, $t) = @{$z->polar};
450 $t = ($t <= 0) ? $t + pi : $t - pi;
451 return (ref $z)->emake($r, $t);
453 my ($re, $im) = @{$z->cartesian};
454 return (ref $z)->make(-$re, -$im);
460 # Compute complex's conjugate.
465 my ($r, $t) = @{$z->polar};
466 return (ref $z)->emake($r, -$t);
468 my ($re, $im) = @{$z->cartesian};
469 return (ref $z)->make($re, -$im);
475 # Compute or set complex's norm (rho).
479 return $z unless ref $z;
481 $z->{'polar'} = [ $rho, ${$z->polar}[1] ];
486 return ${$z->polar}[0];
493 if ($$theta > pi()) { $$theta -= pit2 }
494 elsif ($$theta <= -pi()) { $$theta += pit2 }
500 # Compute or set complex's argument (theta).
503 my ($z, $theta) = @_;
504 return $z unless ref $z;
505 if (defined $theta) {
507 $z->{'polar'} = [ ${$z->polar}[0], $theta ];
511 $theta = ${$z->polar}[1];
522 # It is quite tempting to use wantarray here so that in list context
523 # sqrt() would return the two solutions. This, however, would
526 # print "sqrt(z) = ", sqrt($z), "\n";
528 # The two values would be printed side by side without no intervening
529 # whitespace, quite confusing.
530 # Therefore if you want the two solutions use the root().
534 my ($re, $im) = ref $z ? @{$z->cartesian} : ($z, 0);
535 return $re < 0 ? cplx(0, CORE::sqrt(-$re)) : CORE::sqrt($re) if $im == 0;
536 my ($r, $t) = @{$z->polar};
537 return (ref $z)->emake(CORE::sqrt($r), $t/2);
543 # Compute cbrt(z) (cubic root).
545 # Why are we not returning three values? The same answer as for sqrt().
549 return $z < 0 ? -CORE::exp(CORE::log(-$z)/3) : ($z > 0 ? CORE::exp(CORE::log($z)/3): 0)
551 my ($r, $t) = @{$z->polar};
552 return (ref $z)->emake(CORE::exp(CORE::log($r)/3), $t/3);
561 my $mess = "Root $_[0] not defined, root must be positive integer.\n";
565 $mess .= "Died at $up[1] line $up[2].\n";
573 # Computes all nth root for z, returning an array whose size is n.
574 # `n' must be a positive integer.
576 # The roots are given by (for k = 0..n-1):
578 # z^(1/n) = r^(1/n) (cos ((t+2 k pi)/n) + i sin ((t+2 k pi)/n))
582 _rootbad($n) if ($n < 1 or int($n) != $n);
583 my ($r, $t) = ref $z ? @{$z->polar} : (CORE::abs($z), $z >= 0 ? 0 : pi);
586 my $theta_inc = pit2 / $n;
587 my $rho = $r ** (1/$n);
589 my $cartesian = ref $z && $z->{c_dirty} == 0;
590 for ($k = 0, $theta = $t / $n; $k < $n; $k++, $theta += $theta_inc) {
591 my $w = cplxe($rho, $theta);
592 # Yes, $cartesian is loop invariant.
593 push @root, $cartesian ? cplx(@{$w->cartesian}) : $w;
601 # Return or set Re(z).
605 return $z unless ref $z;
607 $z->{'cartesian'} = [ $Re, ${$z->cartesian}[1] ];
611 return ${$z->cartesian}[0];
618 # Return or set Im(z).
622 return $z unless ref $z;
624 $z->{'cartesian'} = [ ${$z->cartesian}[0], $Im ];
628 return ${$z->cartesian}[1];
635 # Return or set rho(w).
638 Math::Complex::abs(@_);
644 # Return or set theta(w).
647 Math::Complex::arg(@_);
657 my ($x, $y) = @{$z->cartesian};
658 return (ref $z)->emake(CORE::exp($x), $y);
664 # Die on logarithm of zero.
667 my $mess = "$_[0]: Logarithm of zero.\n";
670 $mess .= "(Because in the definition of $_[0], the argument ";
671 $mess .= "$_[1] " unless ($_[1] eq '0');
677 $mess .= "Died at $up[1] line $up[2].\n";
690 _logofzero("log") if $z == 0;
691 return $z > 0 ? CORE::log($z) : cplx(CORE::log(-$z), pi);
693 my ($r, $t) = @{$z->polar};
694 _logofzero("log") if $r == 0;
695 if ($t > pi()) { $t -= pit2 }
696 elsif ($t <= -pi()) { $t += pit2 }
697 return (ref $z)->make(CORE::log($r), $t);
705 sub ln { Math::Complex::log(@_) }
714 return Math::Complex::log($_[0]) * uplog10;
720 # Compute logn(z,n) = log(z) / log(n)
724 $z = cplx($z, 0) unless ref $z;
725 my $logn = $logn{$n};
726 $logn = $logn{$n} = CORE::log($n) unless defined $logn; # Cache log(n)
727 return CORE::log($z) / $logn;
733 # Compute cos(z) = (exp(iz) + exp(-iz))/2.
737 my ($x, $y) = @{$z->cartesian};
738 my $ey = CORE::exp($y);
740 return (ref $z)->make(CORE::cos($x) * ($ey + $ey_1)/2,
741 CORE::sin($x) * ($ey_1 - $ey)/2);
747 # Compute sin(z) = (exp(iz) - exp(-iz))/2.
751 my ($x, $y) = @{$z->cartesian};
752 my $ey = CORE::exp($y);
754 return (ref $z)->make(CORE::sin($x) * ($ey + $ey_1)/2,
755 CORE::cos($x) * ($ey - $ey_1)/2);
761 # Compute tan(z) = sin(z) / cos(z).
765 my $cz = CORE::cos($z);
766 _divbyzero "tan($z)", "cos($z)" if (CORE::abs($cz) < $eps);
767 return CORE::sin($z) / $cz;
773 # Computes the secant sec(z) = 1 / cos(z).
777 my $cz = CORE::cos($z);
778 _divbyzero "sec($z)", "cos($z)" if ($cz == 0);
785 # Computes the cosecant csc(z) = 1 / sin(z).
789 my $sz = CORE::sin($z);
790 _divbyzero "csc($z)", "sin($z)" if ($sz == 0);
799 sub cosec { Math::Complex::csc(@_) }
804 # Computes cot(z) = cos(z) / sin(z).
808 my $sz = CORE::sin($z);
809 _divbyzero "cot($z)", "sin($z)" if ($sz == 0);
810 return CORE::cos($z) / $sz;
818 sub cotan { Math::Complex::cot(@_) }
823 # Computes the arc cosine acos(z) = -i log(z + sqrt(z*z-1)).
827 return CORE::atan2(CORE::sqrt(1-$z*$z), $z) if (! ref $z) && CORE::abs($z) <= 1;
828 my ($x, $y) = ref $z ? @{$z->cartesian} : ($z, 0);
829 my $t1 = CORE::sqrt(($x+1)*($x+1) + $y*$y);
830 my $t2 = CORE::sqrt(($x-1)*($x-1) + $y*$y);
831 my $alpha = ($t1 + $t2)/2;
832 my $beta = ($t1 - $t2)/2;
833 $alpha = 1 if $alpha < 1;
834 if ($beta > 1) { $beta = 1 }
835 elsif ($beta < -1) { $beta = -1 }
836 my $u = CORE::atan2(CORE::sqrt(1-$beta*$beta), $beta);
837 my $v = CORE::log($alpha + CORE::sqrt($alpha*$alpha-1));
838 $v = -$v if $y > 0 || ($y == 0 && $x < -1);
839 return $package->make($u, $v);
845 # Computes the arc sine asin(z) = -i log(iz + sqrt(1-z*z)).
849 return CORE::atan2($z, CORE::sqrt(1-$z*$z)) if (! ref $z) && CORE::abs($z) <= 1;
850 my ($x, $y) = ref $z ? @{$z->cartesian} : ($z, 0);
851 my $t1 = CORE::sqrt(($x+1)*($x+1) + $y*$y);
852 my $t2 = CORE::sqrt(($x-1)*($x-1) + $y*$y);
853 my $alpha = ($t1 + $t2)/2;
854 my $beta = ($t1 - $t2)/2;
855 $alpha = 1 if $alpha < 1;
856 if ($beta > 1) { $beta = 1 }
857 elsif ($beta < -1) { $beta = -1 }
858 my $u = CORE::atan2($beta, CORE::sqrt(1-$beta*$beta));
859 my $v = -CORE::log($alpha + CORE::sqrt($alpha*$alpha-1));
860 $v = -$v if $y > 0 || ($y == 0 && $x < -1);
861 return $package->make($u, $v);
867 # Computes the arc tangent atan(z) = i/2 log((i+z) / (i-z)).
871 return CORE::atan2($z, 1) unless ref $z;
872 _divbyzero "atan(i)" if ( $z == i);
873 _divbyzero "atan(-i)" if (-$z == i);
874 my $log = CORE::log((i + $z) / (i - $z));
875 $ip2 = 0.5 * i unless defined $ip2;
882 # Computes the arc secant asec(z) = acos(1 / z).
886 _divbyzero "asec($z)", $z if ($z == 0);
893 # Computes the arc cosecant acsc(z) = asin(1 / z).
897 _divbyzero "acsc($z)", $z if ($z == 0);
906 sub acosec { Math::Complex::acsc(@_) }
911 # Computes the arc cotangent acot(z) = atan(1 / z)
915 _divbyzero "acot(0)" if (CORE::abs($z) < $eps);
916 return ($z >= 0) ? CORE::atan2(1, $z) : CORE::atan2(-1, -$z) unless ref $z;
917 _divbyzero "acot(i)" if (CORE::abs($z - i) < $eps);
918 _logofzero "acot(-i)" if (CORE::abs($z + i) < $eps);
927 sub acotan { Math::Complex::acot(@_) }
932 # Computes the hyperbolic cosine cosh(z) = (exp(z) + exp(-z))/2.
939 return ($ex + 1/$ex)/2;
941 my ($x, $y) = @{$z->cartesian};
944 return (ref $z)->make(CORE::cos($y) * ($ex + $ex_1)/2,
945 CORE::sin($y) * ($ex - $ex_1)/2);
951 # Computes the hyperbolic sine sinh(z) = (exp(z) - exp(-z))/2.
958 return ($ex - 1/$ex)/2;
960 my ($x, $y) = @{$z->cartesian};
963 return (ref $z)->make(CORE::cos($y) * ($ex - $ex_1)/2,
964 CORE::sin($y) * ($ex + $ex_1)/2);
970 # Computes the hyperbolic tangent tanh(z) = sinh(z) / cosh(z).
975 _divbyzero "tanh($z)", "cosh($z)" if ($cz == 0);
976 return sinh($z) / $cz;
982 # Computes the hyperbolic secant sech(z) = 1 / cosh(z).
987 _divbyzero "sech($z)", "cosh($z)" if ($cz == 0);
994 # Computes the hyperbolic cosecant csch(z) = 1 / sinh(z).
999 _divbyzero "csch($z)", "sinh($z)" if ($sz == 0);
1008 sub cosech { Math::Complex::csch(@_) }
1013 # Computes the hyperbolic cotangent coth(z) = cosh(z) / sinh(z).
1018 _divbyzero "coth($z)", "sinh($z)" if ($sz == 0);
1019 return cosh($z) / $sz;
1027 sub cotanh { Math::Complex::coth(@_) }
1032 # Computes the arc hyperbolic cosine acosh(z) = log(z + sqrt(z*z-1)).
1037 return CORE::log($z + CORE::sqrt($z*$z-1)) if $z >= 1;
1040 my ($re, $im) = @{$z->cartesian};
1042 return cplx(CORE::log($re + CORE::sqrt($re*$re - 1)), 0) if $re >= 1;
1043 return cplx(0, CORE::atan2(CORE::sqrt(1-$re*$re), $re)) if CORE::abs($re) <= 1;
1045 return CORE::log($z + CORE::sqrt($z*$z - 1));
1051 # Computes the arc hyperbolic sine asinh(z) = log(z + sqrt(z*z-1))
1055 return CORE::log($z + CORE::sqrt($z*$z + 1));
1061 # Computes the arc hyperbolic tangent atanh(z) = 1/2 log((1+z) / (1-z)).
1066 return CORE::log((1 + $z)/(1 - $z))/2 if CORE::abs($z) < 1;
1069 _divbyzero 'atanh(1)', "1 - $z" if ($z == 1);
1070 _logofzero 'atanh(-1)' if ($z == -1);
1071 return 0.5 * CORE::log((1 + $z) / (1 - $z));
1077 # Computes the hyperbolic arc secant asech(z) = acosh(1 / z).
1081 _divbyzero 'asech(0)', $z if ($z == 0);
1082 return acosh(1 / $z);
1088 # Computes the hyperbolic arc cosecant acsch(z) = asinh(1 / z).
1092 _divbyzero 'acsch(0)', $z if ($z == 0);
1093 return asinh(1 / $z);
1099 # Alias for acosh().
1101 sub acosech { Math::Complex::acsch(@_) }
1106 # Computes the arc hyperbolic cotangent acoth(z) = 1/2 log((1+z) / (z-1)).
1110 _divbyzero 'acoth(0)' if (CORE::abs($z) < $eps);
1112 return CORE::log(($z + 1)/($z - 1))/2 if CORE::abs($z) > 1;
1115 _divbyzero 'acoth(1)', "$z - 1" if (CORE::abs($z - 1) < $eps);
1116 _logofzero 'acoth(-1)', "1 / $z" if (CORE::abs($z + 1) < $eps);
1117 return CORE::log((1 + $z) / ($z - 1)) / 2;
1125 sub acotanh { Math::Complex::acoth(@_) }
1130 # Compute atan(z1/z2).
1133 my ($z1, $z2, $inverted) = @_;
1134 my ($re1, $im1, $re2, $im2);
1136 ($re1, $im1) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
1137 ($re2, $im2) = @{$z1->cartesian};
1139 ($re1, $im1) = @{$z1->cartesian};
1140 ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
1143 return cplx(CORE::atan2($re1, $re2), 0) if $im1 == 0;
1144 return cplx(($im1<=>0) * pip2, 0) if $re2 == 0;
1146 my $w = atan($z1/$z2);
1147 my ($u, $v) = ref $w ? @{$w->cartesian} : ($w, 0);
1148 $u += pi if $re2 < 0;
1149 $u -= pit2 if $u > pi;
1150 return cplx($u, $v);
1157 # Set (fetch if no argument) display format for all complex numbers that
1158 # don't happen to have overridden it via ->display_format
1160 # When called as a method, this actually sets the display format for
1161 # the current object.
1163 # Valid object formats are 'c' and 'p' for cartesian and polar. The first
1164 # letter is used actually, so the type can be fully spelled out for clarity.
1166 sub display_format {
1170 if (ref $self) { # Called as a method
1172 } else { # Regular procedure call
1177 if (defined $self) {
1178 return defined $self->{display} ? $self->{display} : $display
1179 unless defined $format;
1180 return $self->{display} = $format;
1183 return $display unless defined $format;
1184 return $display = $format;
1190 # Show nicely formatted complex number under its cartesian or polar form,
1191 # depending on the current display format:
1193 # . If a specific display format has been recorded for this object, use it.
1194 # . Otherwise, use the generic current default for all complex numbers,
1195 # which is a package global variable.
1202 $format = $z->{display} if defined $z->{display};
1204 return $z->stringify_polar if $format =~ /^p/i;
1205 return $z->stringify_cartesian;
1209 # ->stringify_cartesian
1211 # Stringify as a cartesian representation 'a+bi'.
1213 sub stringify_cartesian {
1215 my ($x, $y) = @{$z->cartesian};
1218 $x = int($x + ($x < 0 ? -1 : 1) * $eps)
1219 if int(CORE::abs($x)) != int(CORE::abs($x) + $eps);
1220 $y = int($y + ($y < 0 ? -1 : 1) * $eps)
1221 if int(CORE::abs($y)) != int(CORE::abs($y) + $eps);
1223 $re = "$x" if CORE::abs($x) >= $eps;
1224 if ($y == 1) { $im = 'i' }
1225 elsif ($y == -1) { $im = '-i' }
1226 elsif (CORE::abs($y) >= $eps) { $im = $y . "i" }
1229 $str = $re if defined $re;
1230 $str .= "+$im" if defined $im;
1233 $str =~ s/([-+])1i/$1i/; # Not redundant with the above 1/-1 tests.
1234 $str = '0' unless $str;
1240 # Helper for stringify_polar, a Greatest Common Divisor with a memory.
1247 # Loops forever if given negative inputs.
1249 if ($b and $a > $b) { return gcd($a % $b, $b) }
1250 elsif ($a and $b > $a) { return gcd($b % $a, $a) }
1251 else { return $a ? $a : $b }
1261 unless (exists $gcd{$id}) {
1262 $gcd{$id} = _gcd($a, $b);
1263 $gcd{"$b $a"} = $gcd{$id};
1272 # Stringify as a polar representation '[r,t]'.
1274 sub stringify_polar {
1276 my ($r, $t) = @{$z->polar};
1279 return '[0,0]' if $r <= $eps;
1282 $nt = ($nt - int($nt)) * pit2;
1283 $nt += pit2 if $nt < 0; # Range [0, 2pi]
1285 if (CORE::abs($nt) <= $eps) { $theta = 0 }
1286 elsif (CORE::abs(pi-$nt) <= $eps) { $theta = 'pi' }
1288 if (defined $theta) {
1289 $r = int($r + ($r < 0 ? -1 : 1) * $eps)
1290 if int(CORE::abs($r)) != int(CORE::abs($r) + $eps);
1291 $theta = int($theta + ($theta < 0 ? -1 : 1) * $eps)
1292 if ($theta ne 'pi' and
1293 int(CORE::abs($theta)) != int(CORE::abs($theta) + $eps));
1294 return "\[$r,$theta\]";
1298 # Okay, number is not a real. Try to identify pi/n and friends...
1301 $nt -= pit2 if $nt > pi;
1303 if (CORE::abs($nt) >= deg1) {
1306 for ($k = 1, $kpi = pi; $k < 10; $k++, $kpi += pi) {
1307 $n = int($kpi / $nt + ($nt > 0 ? 1 : -1) * 0.5);
1308 if (CORE::abs($kpi/$n - $nt) <= $eps) {
1310 my $gcd = gcd($k, $n);
1316 $theta = ($nt < 0 ? '-':'').
1317 ($k == 1 ? 'pi':"${k}pi");
1318 $theta .= '/'.$n if $n > 1;
1324 $theta = $nt unless defined $theta;
1326 $r = int($r + ($r < 0 ? -1 : 1) * $eps)
1327 if int(CORE::abs($r)) != int(CORE::abs($r) + $eps);
1328 $theta = int($theta + ($theta < 0 ? -1 : 1) * $eps)
1329 if ($theta !~ m(^-?\d*pi/\d+$) and
1330 int(CORE::abs($theta)) != int(CORE::abs($theta) + $eps));
1332 return "\[$r,$theta\]";
1340 Math::Complex - complex numbers and associated mathematical functions
1346 $z = Math::Complex->make(5, 6);
1348 $j = cplxe(1, 2*pi/3);
1352 This package lets you create and manipulate complex numbers. By default,
1353 I<Perl> limits itself to real numbers, but an extra C<use> statement brings
1354 full complex support, along with a full set of mathematical functions
1355 typically associated with and/or extended to complex numbers.
1357 If you wonder what complex numbers are, they were invented to be able to solve
1358 the following equation:
1362 and by definition, the solution is noted I<i> (engineers use I<j> instead since
1363 I<i> usually denotes an intensity, but the name does not matter). The number
1364 I<i> is a pure I<imaginary> number.
1366 The arithmetics with pure imaginary numbers works just like you would expect
1367 it with real numbers... you just have to remember that
1373 5i + 7i = i * (5 + 7) = 12i
1374 4i - 3i = i * (4 - 3) = i
1379 Complex numbers are numbers that have both a real part and an imaginary
1380 part, and are usually noted:
1384 where C<a> is the I<real> part and C<b> is the I<imaginary> part. The
1385 arithmetic with complex numbers is straightforward. You have to
1386 keep track of the real and the imaginary parts, but otherwise the
1387 rules used for real numbers just apply:
1389 (4 + 3i) + (5 - 2i) = (4 + 5) + i(3 - 2) = 9 + i
1390 (2 + i) * (4 - i) = 2*4 + 4i -2i -i*i = 8 + 2i + 1 = 9 + 2i
1392 A graphical representation of complex numbers is possible in a plane
1393 (also called the I<complex plane>, but it's really a 2D plane).
1398 is the point whose coordinates are (a, b). Actually, it would
1399 be the vector originating from (0, 0) to (a, b). It follows that the addition
1400 of two complex numbers is a vectorial addition.
1402 Since there is a bijection between a point in the 2D plane and a complex
1403 number (i.e. the mapping is unique and reciprocal), a complex number
1404 can also be uniquely identified with polar coordinates:
1408 where C<rho> is the distance to the origin, and C<theta> the angle between
1409 the vector and the I<x> axis. There is a notation for this using the
1410 exponential form, which is:
1412 rho * exp(i * theta)
1414 where I<i> is the famous imaginary number introduced above. Conversion
1415 between this form and the cartesian form C<a + bi> is immediate:
1417 a = rho * cos(theta)
1418 b = rho * sin(theta)
1420 which is also expressed by this formula:
1422 z = rho * exp(i * theta) = rho * (cos theta + i * sin theta)
1424 In other words, it's the projection of the vector onto the I<x> and I<y>
1425 axes. Mathematicians call I<rho> the I<norm> or I<modulus> and I<theta>
1426 the I<argument> of the complex number. The I<norm> of C<z> will be
1429 The polar notation (also known as the trigonometric
1430 representation) is much more handy for performing multiplications and
1431 divisions of complex numbers, whilst the cartesian notation is better
1432 suited for additions and subtractions. Real numbers are on the I<x>
1433 axis, and therefore I<theta> is zero or I<pi>.
1435 All the common operations that can be performed on a real number have
1436 been defined to work on complex numbers as well, and are merely
1437 I<extensions> of the operations defined on real numbers. This means
1438 they keep their natural meaning when there is no imaginary part, provided
1439 the number is within their definition set.
1441 For instance, the C<sqrt> routine which computes the square root of
1442 its argument is only defined for non-negative real numbers and yields a
1443 non-negative real number (it is an application from B<R+> to B<R+>).
1444 If we allow it to return a complex number, then it can be extended to
1445 negative real numbers to become an application from B<R> to B<C> (the
1446 set of complex numbers):
1448 sqrt(x) = x >= 0 ? sqrt(x) : sqrt(-x)*i
1450 It can also be extended to be an application from B<C> to B<C>,
1451 whilst its restriction to B<R> behaves as defined above by using
1452 the following definition:
1454 sqrt(z = [r,t]) = sqrt(r) * exp(i * t/2)
1456 Indeed, a negative real number can be noted C<[x,pi]> (the modulus
1457 I<x> is always non-negative, so C<[x,pi]> is really C<-x>, a negative
1458 number) and the above definition states that
1460 sqrt([x,pi]) = sqrt(x) * exp(i*pi/2) = [sqrt(x),pi/2] = sqrt(x)*i
1462 which is exactly what we had defined for negative real numbers above.
1463 The C<sqrt> returns only one of the solutions: if you want the both,
1464 use the C<root> function.
1466 All the common mathematical functions defined on real numbers that
1467 are extended to complex numbers share that same property of working
1468 I<as usual> when the imaginary part is zero (otherwise, it would not
1469 be called an extension, would it?).
1471 A I<new> operation possible on a complex number that is
1472 the identity for real numbers is called the I<conjugate>, and is noted
1473 with an horizontal bar above the number, or C<~z> here.
1480 z * ~z = (a + bi) * (a - bi) = a*a + b*b
1482 We saw that the norm of C<z> was noted C<abs(z)> and was defined as the
1483 distance to the origin, also known as:
1485 rho = abs(z) = sqrt(a*a + b*b)
1489 z * ~z = abs(z) ** 2
1491 If z is a pure real number (i.e. C<b == 0>), then the above yields:
1495 which is true (C<abs> has the regular meaning for real number, i.e. stands
1496 for the absolute value). This example explains why the norm of C<z> is
1497 noted C<abs(z)>: it extends the C<abs> function to complex numbers, yet
1498 is the regular C<abs> we know when the complex number actually has no
1499 imaginary part... This justifies I<a posteriori> our use of the C<abs>
1500 notation for the norm.
1504 Given the following notations:
1506 z1 = a + bi = r1 * exp(i * t1)
1507 z2 = c + di = r2 * exp(i * t2)
1508 z = <any complex or real number>
1510 the following (overloaded) operations are supported on complex numbers:
1512 z1 + z2 = (a + c) + i(b + d)
1513 z1 - z2 = (a - c) + i(b - d)
1514 z1 * z2 = (r1 * r2) * exp(i * (t1 + t2))
1515 z1 / z2 = (r1 / r2) * exp(i * (t1 - t2))
1516 z1 ** z2 = exp(z2 * log z1)
1518 abs(z) = r1 = sqrt(a*a + b*b)
1519 sqrt(z) = sqrt(r1) * exp(i * t/2)
1520 exp(z) = exp(a) * exp(i * b)
1521 log(z) = log(r1) + i*t
1522 sin(z) = 1/2i (exp(i * z1) - exp(-i * z))
1523 cos(z) = 1/2 (exp(i * z1) + exp(-i * z))
1524 atan2(z1, z2) = atan(z1/z2)
1526 The following extra operations are supported on both real and complex
1534 cbrt(z) = z ** (1/3)
1535 log10(z) = log(z) / log(10)
1536 logn(z, n) = log(z) / log(n)
1538 tan(z) = sin(z) / cos(z)
1544 asin(z) = -i * log(i*z + sqrt(1-z*z))
1545 acos(z) = -i * log(z + i*sqrt(1-z*z))
1546 atan(z) = i/2 * log((i+z) / (i-z))
1548 acsc(z) = asin(1 / z)
1549 asec(z) = acos(1 / z)
1550 acot(z) = atan(1 / z) = -i/2 * log((i+z) / (z-i))
1552 sinh(z) = 1/2 (exp(z) - exp(-z))
1553 cosh(z) = 1/2 (exp(z) + exp(-z))
1554 tanh(z) = sinh(z) / cosh(z) = (exp(z) - exp(-z)) / (exp(z) + exp(-z))
1556 csch(z) = 1 / sinh(z)
1557 sech(z) = 1 / cosh(z)
1558 coth(z) = 1 / tanh(z)
1560 asinh(z) = log(z + sqrt(z*z+1))
1561 acosh(z) = log(z + sqrt(z*z-1))
1562 atanh(z) = 1/2 * log((1+z) / (1-z))
1564 acsch(z) = asinh(1 / z)
1565 asech(z) = acosh(1 / z)
1566 acoth(z) = atanh(1 / z) = 1/2 * log((1+z) / (z-1))
1568 I<arg>, I<abs>, I<log>, I<csc>, I<cot>, I<acsc>, I<acot>, I<csch>,
1569 I<coth>, I<acosech>, I<acotanh>, have aliases I<rho>, I<theta>, I<ln>,
1570 I<cosec>, I<cotan>, I<acosec>, I<acotan>, I<cosech>, I<cotanh>,
1571 I<acosech>, I<acotanh>, respectively. C<Re>, C<Im>, C<arg>, C<abs>,
1572 C<rho>, and C<theta> can be used also also mutators. The C<cbrt>
1573 returns only one of the solutions: if you want all three, use the
1576 The I<root> function is available to compute all the I<n>
1577 roots of some complex, where I<n> is a strictly positive integer.
1578 There are exactly I<n> such roots, returned as a list. Getting the
1579 number mathematicians call C<j> such that:
1583 is a simple matter of writing:
1585 $j = ((root(1, 3))[1];
1587 The I<k>th root for C<z = [r,t]> is given by:
1589 (root(z, n))[k] = r**(1/n) * exp(i * (t + 2*k*pi)/n)
1591 The I<spaceship> comparison operator, E<lt>=E<gt>, is also defined. In
1592 order to ensure its restriction to real numbers is conform to what you
1593 would expect, the comparison is run on the real part of the complex
1594 number first, and imaginary parts are compared only when the real
1599 To create a complex number, use either:
1601 $z = Math::Complex->make(3, 4);
1604 if you know the cartesian form of the number, or
1608 if you like. To create a number using the polar form, use either:
1610 $z = Math::Complex->emake(5, pi/3);
1611 $x = cplxe(5, pi/3);
1613 instead. The first argument is the modulus, the second is the angle
1614 (in radians, the full circle is 2*pi). (Mnemonic: C<e> is used as a
1615 notation for complex numbers in the polar form).
1617 It is possible to write:
1619 $x = cplxe(-3, pi/4);
1621 but that will be silently converted into C<[3,-3pi/4]>, since the modulus
1622 must be non-negative (it represents the distance to the origin in the complex
1625 It is also possible to have a complex number as either argument of
1626 either the C<make> or C<emake>: the appropriate component of
1627 the argument will be used.
1632 =head1 STRINGIFICATION
1634 When printed, a complex number is usually shown under its cartesian
1635 form I<a+bi>, but there are legitimate cases where the polar format
1636 I<[r,t]> is more appropriate.
1638 By calling the routine C<Math::Complex::display_format> and supplying either
1639 C<"polar"> or C<"cartesian">, you override the default display format,
1640 which is C<"cartesian">. Not supplying any argument returns the current
1643 This default can be overridden on a per-number basis by calling the
1644 C<display_format> method instead. As before, not supplying any argument
1645 returns the current display format for this number. Otherwise whatever you
1646 specify will be the new display format for I<this> particular number.
1652 Math::Complex::display_format('polar');
1653 $j = ((root(1, 3))[1];
1654 print "j = $j\n"; # Prints "j = [1,2pi/3]
1655 $j->display_format('cartesian');
1656 print "j = $j\n"; # Prints "j = -0.5+0.866025403784439i"
1658 The polar format attempts to emphasize arguments like I<k*pi/n>
1659 (where I<n> is a positive integer and I<k> an integer within [-9,+9]).
1663 Thanks to overloading, the handling of arithmetics with complex numbers
1664 is simple and almost transparent.
1666 Here are some examples:
1670 $j = cplxe(1, 2*pi/3); # $j ** 3 == 1
1671 print "j = $j, j**3 = ", $j ** 3, "\n";
1672 print "1 + j + j**2 = ", 1 + $j + $j**2, "\n";
1674 $z = -16 + 0*i; # Force it to be a complex
1675 print "sqrt($z) = ", sqrt($z), "\n";
1677 $k = exp(i * 2*pi/3);
1678 print "$j - $k = ", $j - $k, "\n";
1680 $z->Re(3); # Re, Im, arg, abs,
1681 $j->arg(2); # (the last two aka rho, theta)
1682 # can be used also as mutators.
1684 =head1 ERRORS DUE TO DIVISION BY ZERO OR LOGARITHM OF ZERO
1686 The division (/) and the following functions
1692 atanh asech acsch acoth
1694 cannot be computed for all arguments because that would mean dividing
1695 by zero or taking logarithm of zero. These situations cause fatal
1696 runtime errors looking like this
1698 cot(0): Division by zero.
1699 (Because in the definition of cot(0), the divisor sin(0) is 0)
1704 atanh(-1): Logarithm of zero.
1707 For the C<csc>, C<cot>, C<asec>, C<acsc>, C<acot>, C<csch>, C<coth>,
1708 C<asech>, C<acsch>, the argument cannot be C<0> (zero). For the the
1709 logarithmic functions and the C<atanh>, C<acoth>, the argument cannot
1710 be C<1> (one). For the C<atanh>, C<acoth>, the argument cannot be
1711 C<-1> (minus one). For the C<atan>, C<acot>, the argument cannot be
1712 C<i> (the imaginary unit). For the C<atan>, C<acoth>, the argument
1713 cannot be C<-i> (the negative imaginary unit). For the C<tan>,
1714 C<sec>, C<tanh>, the argument cannot be I<pi/2 + k * pi>, where I<k>
1717 Note that because we are operating on approximations of real numbers,
1718 these errors can happen when merely `too close' to the singularities
1719 listed above. For example C<tan(2*atan2(1,1)+1e-15)> will die of
1722 =head1 ERRORS DUE TO INDIGESTIBLE ARGUMENTS
1724 The C<make> and C<emake> accept both real and complex arguments.
1725 When they cannot recognize the arguments they will die with error
1726 messages like the following
1728 Math::Complex::make: Cannot take real part of ...
1729 Math::Complex::make: Cannot take real part of ...
1730 Math::Complex::emake: Cannot take rho of ...
1731 Math::Complex::emake: Cannot take theta of ...
1735 Saying C<use Math::Complex;> exports many mathematical routines in the
1736 caller environment and even overrides some (C<sqrt>, C<log>).
1737 This is construed as a feature by the Authors, actually... ;-)
1739 All routines expect to be given real or complex numbers. Don't attempt to
1740 use BigFloat, since Perl has currently no rule to disambiguate a '+'
1741 operation (for instance) between two overloaded entities.
1743 In Cray UNICOS there is some strange numerical instability that results
1744 in root(), cos(), sin(), cosh(), sinh(), losing accuracy fast. Beware.
1745 The bug may be in UNICOS math libs, in UNICOS C compiler, in Math::Complex.
1746 Whatever it is, it does not manifest itself anywhere else where Perl runs.
1750 Raphael Manfredi <F<Raphael_Manfredi@pobox.com>> and
1751 Jarkko Hietaniemi <F<jhi@iki.fi>>.
1753 Extensive patches by Daniel S. Lewart <F<d-lewart@uiuc.edu>>.