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_FORMAT = ('style' => 'cartesian',
71 'polar_pretty_print' => 1);
72 my $eps = 1e-14; # Epsilon
75 # Object attributes (internal):
76 # cartesian [real, imaginary] -- cartesian form
77 # polar [rho, theta] -- polar form
78 # c_dirty cartesian form not up-to-date
79 # p_dirty polar form not up-to-date
80 # display display format (package's global when not set)
83 # Die on bad *make() arguments.
86 die "@{[(caller(1))[3]]}: Cannot take $_[0] of $_[1].\n";
92 # Create a new complex number (cartesian form)
95 my $self = bless {}, shift;
99 if ( $rre eq ref $self ) {
102 _cannot_make("real part", $rre);
107 if ( $rim eq ref $self ) {
110 _cannot_make("imaginary part", $rim);
113 $self->{'cartesian'} = [ $re, $im ];
114 $self->{c_dirty} = 0;
115 $self->{p_dirty} = 1;
116 $self->display_format('cartesian');
123 # Create a new complex number (exponential form)
126 my $self = bless {}, shift;
127 my ($rho, $theta) = @_;
130 if ( $rrh eq ref $self ) {
133 _cannot_make("rho", $rrh);
136 my $rth = ref $theta;
138 if ( $rth eq ref $self ) {
139 $theta = theta($theta);
141 _cannot_make("theta", $rth);
146 $theta = ($theta <= 0) ? $theta + pi() : $theta - pi();
148 $self->{'polar'} = [$rho, $theta];
149 $self->{p_dirty} = 0;
150 $self->{c_dirty} = 1;
151 $self->display_format('polar');
155 sub new { &make } # For backward compatibility only.
160 # Creates a complex number from a (re, im) tuple.
161 # This avoids the burden of writing Math::Complex->make(re, im).
165 return __PACKAGE__->make($re, defined $im ? $im : 0);
171 # Creates a complex number from a (rho, theta) tuple.
172 # This avoids the burden of writing Math::Complex->emake(rho, theta).
175 my ($rho, $theta) = @_;
176 return __PACKAGE__->emake($rho, defined $theta ? $theta : 0);
182 # The number defined as pi = 180 degrees
184 sub pi () { 4 * CORE::atan2(1, 1) }
191 sub pit2 () { 2 * pi }
198 sub pip2 () { pi / 2 }
203 # One degree in radians, used in stringify_polar.
206 sub deg1 () { pi / 180 }
213 sub uplog10 () { 1 / CORE::log(10) }
218 # The number defined as i*i = -1;
223 $i->{'cartesian'} = [0, 1];
224 $i->{'polar'} = [1, pip2];
231 # Attribute access/set routines
234 sub cartesian {$_[0]->{c_dirty} ?
235 $_[0]->update_cartesian : $_[0]->{'cartesian'}}
236 sub polar {$_[0]->{p_dirty} ?
237 $_[0]->update_polar : $_[0]->{'polar'}}
239 sub set_cartesian { $_[0]->{p_dirty}++; $_[0]->{'cartesian'} = $_[1] }
240 sub set_polar { $_[0]->{c_dirty}++; $_[0]->{'polar'} = $_[1] }
245 # Recompute and return the cartesian form, given accurate polar form.
247 sub update_cartesian {
249 my ($r, $t) = @{$self->{'polar'}};
250 $self->{c_dirty} = 0;
251 return $self->{'cartesian'} = [$r * CORE::cos($t), $r * CORE::sin($t)];
258 # Recompute and return the polar form, given accurate cartesian form.
262 my ($x, $y) = @{$self->{'cartesian'}};
263 $self->{p_dirty} = 0;
264 return $self->{'polar'} = [0, 0] if $x == 0 && $y == 0;
265 return $self->{'polar'} = [CORE::sqrt($x*$x + $y*$y), CORE::atan2($y, $x)];
274 my ($z1, $z2, $regular) = @_;
275 my ($re1, $im1) = @{$z1->cartesian};
276 $z2 = cplx($z2) unless ref $z2;
277 my ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
278 unless (defined $regular) {
279 $z1->set_cartesian([$re1 + $re2, $im1 + $im2]);
282 return (ref $z1)->make($re1 + $re2, $im1 + $im2);
291 my ($z1, $z2, $inverted) = @_;
292 my ($re1, $im1) = @{$z1->cartesian};
293 $z2 = cplx($z2) unless ref $z2;
294 my ($re2, $im2) = @{$z2->cartesian};
295 unless (defined $inverted) {
296 $z1->set_cartesian([$re1 - $re2, $im1 - $im2]);
300 (ref $z1)->make($re2 - $re1, $im2 - $im1) :
301 (ref $z1)->make($re1 - $re2, $im1 - $im2);
311 my ($z1, $z2, $regular) = @_;
312 if ($z1->{p_dirty} == 0 and ref $z2 and $z2->{p_dirty} == 0) {
313 # if both polar better use polar to avoid rounding errors
314 my ($r1, $t1) = @{$z1->polar};
315 my ($r2, $t2) = @{$z2->polar};
317 if ($t > pi()) { $t -= pit2 }
318 elsif ($t <= -pi()) { $t += pit2 }
319 unless (defined $regular) {
320 $z1->set_polar([$r1 * $r2, $t]);
323 return (ref $z1)->emake($r1 * $r2, $t);
325 my ($x1, $y1) = @{$z1->cartesian};
327 my ($x2, $y2) = @{$z2->cartesian};
328 return (ref $z1)->make($x1*$x2-$y1*$y2, $x1*$y2+$y1*$x2);
330 return (ref $z1)->make($x1*$z2, $y1*$z2);
338 # Die on division by zero.
341 my $mess = "$_[0]: Division by zero.\n";
344 $mess .= "(Because in the definition of $_[0], the divisor ";
345 $mess .= "$_[1] " unless ($_[1] eq '0');
351 $mess .= "Died at $up[1] line $up[2].\n";
362 my ($z1, $z2, $inverted) = @_;
363 if ($z1->{p_dirty} == 0 and ref $z2 and $z2->{p_dirty} == 0) {
364 # if both polar better use polar to avoid rounding errors
365 my ($r1, $t1) = @{$z1->polar};
366 my ($r2, $t2) = @{$z2->polar};
369 _divbyzero "$z2/0" if ($r1 == 0);
371 if ($t > pi()) { $t -= pit2 }
372 elsif ($t <= -pi()) { $t += pit2 }
373 return (ref $z1)->emake($r2 / $r1, $t);
375 _divbyzero "$z1/0" if ($r2 == 0);
377 if ($t > pi()) { $t -= pit2 }
378 elsif ($t <= -pi()) { $t += pit2 }
379 return (ref $z1)->emake($r1 / $r2, $t);
384 ($x2, $y2) = @{$z1->cartesian};
385 $d = $x2*$x2 + $y2*$y2;
386 _divbyzero "$z2/0" if $d == 0;
387 return (ref $z1)->make(($x2*$z2)/$d, -($y2*$z2)/$d);
389 my ($x1, $y1) = @{$z1->cartesian};
391 ($x2, $y2) = @{$z2->cartesian};
392 $d = $x2*$x2 + $y2*$y2;
393 _divbyzero "$z1/0" if $d == 0;
394 my $u = ($x1*$x2 + $y1*$y2)/$d;
395 my $v = ($y1*$x2 - $x1*$y2)/$d;
396 return (ref $z1)->make($u, $v);
398 _divbyzero "$z1/0" if $z2 == 0;
399 return (ref $z1)->make($x1/$z2, $y1/$z2);
408 # Computes z1**z2 = exp(z2 * log z1)).
411 my ($z1, $z2, $inverted) = @_;
413 return 1 if $z1 == 0 || $z2 == 1;
414 return 0 if $z2 == 0 && Re($z1) > 0;
416 return 1 if $z2 == 0 || $z1 == 1;
417 return 0 if $z1 == 0 && Re($z2) > 0;
419 my $w = $inverted ? CORE::exp($z1 * CORE::log($z2))
420 : CORE::exp($z2 * CORE::log($z1));
421 # If both arguments cartesian, return cartesian, else polar.
422 return $z1->{c_dirty} == 0 &&
423 (not ref $z2 or $z2->{c_dirty} == 0) ?
424 cplx(@{$w->cartesian}) : $w;
430 # Computes z1 <=> z2.
431 # Sorts on the real part first, then on the imaginary part. Thus 2-4i < 3+8i.
434 my ($z1, $z2, $inverted) = @_;
435 my ($re1, $im1) = ref $z1 ? @{$z1->cartesian} : ($z1, 0);
436 my ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
437 my $sgn = $inverted ? -1 : 1;
438 return $sgn * ($re1 <=> $re2) if $re1 != $re2;
439 return $sgn * ($im1 <=> $im2);
450 my ($r, $t) = @{$z->polar};
451 $t = ($t <= 0) ? $t + pi : $t - pi;
452 return (ref $z)->emake($r, $t);
454 my ($re, $im) = @{$z->cartesian};
455 return (ref $z)->make(-$re, -$im);
461 # Compute complex's conjugate.
466 my ($r, $t) = @{$z->polar};
467 return (ref $z)->emake($r, -$t);
469 my ($re, $im) = @{$z->cartesian};
470 return (ref $z)->make($re, -$im);
476 # Compute or set complex's norm (rho).
480 return $z unless ref $z;
482 $z->{'polar'} = [ $rho, ${$z->polar}[1] ];
487 return ${$z->polar}[0];
494 if ($$theta > pi()) { $$theta -= pit2 }
495 elsif ($$theta <= -pi()) { $$theta += pit2 }
501 # Compute or set complex's argument (theta).
504 my ($z, $theta) = @_;
505 return $z unless ref $z;
506 if (defined $theta) {
508 $z->{'polar'} = [ ${$z->polar}[0], $theta ];
512 $theta = ${$z->polar}[1];
523 # It is quite tempting to use wantarray here so that in list context
524 # sqrt() would return the two solutions. This, however, would
527 # print "sqrt(z) = ", sqrt($z), "\n";
529 # The two values would be printed side by side without no intervening
530 # whitespace, quite confusing.
531 # Therefore if you want the two solutions use the root().
535 my ($re, $im) = ref $z ? @{$z->cartesian} : ($z, 0);
536 return $re < 0 ? cplx(0, CORE::sqrt(-$re)) : CORE::sqrt($re) if $im == 0;
537 my ($r, $t) = @{$z->polar};
538 return (ref $z)->emake(CORE::sqrt($r), $t/2);
544 # Compute cbrt(z) (cubic root).
546 # Why are we not returning three values? The same answer as for sqrt().
550 return $z < 0 ? -CORE::exp(CORE::log(-$z)/3) : ($z > 0 ? CORE::exp(CORE::log($z)/3): 0)
552 my ($r, $t) = @{$z->polar};
553 return (ref $z)->emake(CORE::exp(CORE::log($r)/3), $t/3);
562 my $mess = "Root $_[0] not defined, root must be positive integer.\n";
566 $mess .= "Died at $up[1] line $up[2].\n";
574 # Computes all nth root for z, returning an array whose size is n.
575 # `n' must be a positive integer.
577 # The roots are given by (for k = 0..n-1):
579 # z^(1/n) = r^(1/n) (cos ((t+2 k pi)/n) + i sin ((t+2 k pi)/n))
583 _rootbad($n) if ($n < 1 or int($n) != $n);
584 my ($r, $t) = ref $z ? @{$z->polar} : (CORE::abs($z), $z >= 0 ? 0 : pi);
587 my $theta_inc = pit2 / $n;
588 my $rho = $r ** (1/$n);
590 my $cartesian = ref $z && $z->{c_dirty} == 0;
591 for ($k = 0, $theta = $t / $n; $k < $n; $k++, $theta += $theta_inc) {
592 my $w = cplxe($rho, $theta);
593 # Yes, $cartesian is loop invariant.
594 push @root, $cartesian ? cplx(@{$w->cartesian}) : $w;
602 # Return or set Re(z).
606 return $z unless ref $z;
608 $z->{'cartesian'} = [ $Re, ${$z->cartesian}[1] ];
612 return ${$z->cartesian}[0];
619 # Return or set Im(z).
623 return $z unless ref $z;
625 $z->{'cartesian'} = [ ${$z->cartesian}[0], $Im ];
629 return ${$z->cartesian}[1];
636 # Return or set rho(w).
639 Math::Complex::abs(@_);
645 # Return or set theta(w).
648 Math::Complex::arg(@_);
658 my ($x, $y) = @{$z->cartesian};
659 return (ref $z)->emake(CORE::exp($x), $y);
665 # Die on logarithm of zero.
668 my $mess = "$_[0]: Logarithm of zero.\n";
671 $mess .= "(Because in the definition of $_[0], the argument ";
672 $mess .= "$_[1] " unless ($_[1] eq '0');
678 $mess .= "Died at $up[1] line $up[2].\n";
691 _logofzero("log") if $z == 0;
692 return $z > 0 ? CORE::log($z) : cplx(CORE::log(-$z), pi);
694 my ($r, $t) = @{$z->polar};
695 _logofzero("log") if $r == 0;
696 if ($t > pi()) { $t -= pit2 }
697 elsif ($t <= -pi()) { $t += pit2 }
698 return (ref $z)->make(CORE::log($r), $t);
706 sub ln { Math::Complex::log(@_) }
715 return Math::Complex::log($_[0]) * uplog10;
721 # Compute logn(z,n) = log(z) / log(n)
725 $z = cplx($z, 0) unless ref $z;
726 my $logn = $logn{$n};
727 $logn = $logn{$n} = CORE::log($n) unless defined $logn; # Cache log(n)
728 return CORE::log($z) / $logn;
734 # Compute cos(z) = (exp(iz) + exp(-iz))/2.
738 my ($x, $y) = @{$z->cartesian};
739 my $ey = CORE::exp($y);
741 return (ref $z)->make(CORE::cos($x) * ($ey + $ey_1)/2,
742 CORE::sin($x) * ($ey_1 - $ey)/2);
748 # Compute sin(z) = (exp(iz) - exp(-iz))/2.
752 my ($x, $y) = @{$z->cartesian};
753 my $ey = CORE::exp($y);
755 return (ref $z)->make(CORE::sin($x) * ($ey + $ey_1)/2,
756 CORE::cos($x) * ($ey - $ey_1)/2);
762 # Compute tan(z) = sin(z) / cos(z).
766 my $cz = CORE::cos($z);
767 _divbyzero "tan($z)", "cos($z)" if (CORE::abs($cz) < $eps);
768 return CORE::sin($z) / $cz;
774 # Computes the secant sec(z) = 1 / cos(z).
778 my $cz = CORE::cos($z);
779 _divbyzero "sec($z)", "cos($z)" if ($cz == 0);
786 # Computes the cosecant csc(z) = 1 / sin(z).
790 my $sz = CORE::sin($z);
791 _divbyzero "csc($z)", "sin($z)" if ($sz == 0);
800 sub cosec { Math::Complex::csc(@_) }
805 # Computes cot(z) = cos(z) / sin(z).
809 my $sz = CORE::sin($z);
810 _divbyzero "cot($z)", "sin($z)" if ($sz == 0);
811 return CORE::cos($z) / $sz;
819 sub cotan { Math::Complex::cot(@_) }
824 # Computes the arc cosine acos(z) = -i log(z + sqrt(z*z-1)).
828 return CORE::atan2(CORE::sqrt(1-$z*$z), $z) if (! ref $z) && CORE::abs($z) <= 1;
829 my ($x, $y) = ref $z ? @{$z->cartesian} : ($z, 0);
830 my $t1 = CORE::sqrt(($x+1)*($x+1) + $y*$y);
831 my $t2 = CORE::sqrt(($x-1)*($x-1) + $y*$y);
832 my $alpha = ($t1 + $t2)/2;
833 my $beta = ($t1 - $t2)/2;
834 $alpha = 1 if $alpha < 1;
835 if ($beta > 1) { $beta = 1 }
836 elsif ($beta < -1) { $beta = -1 }
837 my $u = CORE::atan2(CORE::sqrt(1-$beta*$beta), $beta);
838 my $v = CORE::log($alpha + CORE::sqrt($alpha*$alpha-1));
839 $v = -$v if $y > 0 || ($y == 0 && $x < -1);
840 return __PACKAGE__->make($u, $v);
846 # Computes the arc sine asin(z) = -i log(iz + sqrt(1-z*z)).
850 return CORE::atan2($z, CORE::sqrt(1-$z*$z)) if (! ref $z) && CORE::abs($z) <= 1;
851 my ($x, $y) = ref $z ? @{$z->cartesian} : ($z, 0);
852 my $t1 = CORE::sqrt(($x+1)*($x+1) + $y*$y);
853 my $t2 = CORE::sqrt(($x-1)*($x-1) + $y*$y);
854 my $alpha = ($t1 + $t2)/2;
855 my $beta = ($t1 - $t2)/2;
856 $alpha = 1 if $alpha < 1;
857 if ($beta > 1) { $beta = 1 }
858 elsif ($beta < -1) { $beta = -1 }
859 my $u = CORE::atan2($beta, CORE::sqrt(1-$beta*$beta));
860 my $v = -CORE::log($alpha + CORE::sqrt($alpha*$alpha-1));
861 $v = -$v if $y > 0 || ($y == 0 && $x < -1);
862 return __PACKAGE__->make($u, $v);
868 # Computes the arc tangent atan(z) = i/2 log((i+z) / (i-z)).
872 return CORE::atan2($z, 1) unless ref $z;
873 _divbyzero "atan(i)" if ( $z == i);
874 _divbyzero "atan(-i)" if (-$z == i);
875 my $log = CORE::log((i + $z) / (i - $z));
876 $ip2 = 0.5 * i unless defined $ip2;
883 # Computes the arc secant asec(z) = acos(1 / z).
887 _divbyzero "asec($z)", $z if ($z == 0);
894 # Computes the arc cosecant acsc(z) = asin(1 / z).
898 _divbyzero "acsc($z)", $z if ($z == 0);
907 sub acosec { Math::Complex::acsc(@_) }
912 # Computes the arc cotangent acot(z) = atan(1 / z)
916 _divbyzero "acot(0)" if (CORE::abs($z) < $eps);
917 return ($z >= 0) ? CORE::atan2(1, $z) : CORE::atan2(-1, -$z) unless ref $z;
918 _divbyzero "acot(i)" if (CORE::abs($z - i) < $eps);
919 _logofzero "acot(-i)" if (CORE::abs($z + i) < $eps);
928 sub acotan { Math::Complex::acot(@_) }
933 # Computes the hyperbolic cosine cosh(z) = (exp(z) + exp(-z))/2.
940 return ($ex + 1/$ex)/2;
942 my ($x, $y) = @{$z->cartesian};
945 return (ref $z)->make(CORE::cos($y) * ($ex + $ex_1)/2,
946 CORE::sin($y) * ($ex - $ex_1)/2);
952 # Computes the hyperbolic sine sinh(z) = (exp(z) - exp(-z))/2.
959 return ($ex - 1/$ex)/2;
961 my ($x, $y) = @{$z->cartesian};
964 return (ref $z)->make(CORE::cos($y) * ($ex - $ex_1)/2,
965 CORE::sin($y) * ($ex + $ex_1)/2);
971 # Computes the hyperbolic tangent tanh(z) = sinh(z) / cosh(z).
976 _divbyzero "tanh($z)", "cosh($z)" if ($cz == 0);
977 return sinh($z) / $cz;
983 # Computes the hyperbolic secant sech(z) = 1 / cosh(z).
988 _divbyzero "sech($z)", "cosh($z)" if ($cz == 0);
995 # Computes the hyperbolic cosecant csch(z) = 1 / sinh(z).
1000 _divbyzero "csch($z)", "sinh($z)" if ($sz == 0);
1009 sub cosech { Math::Complex::csch(@_) }
1014 # Computes the hyperbolic cotangent coth(z) = cosh(z) / sinh(z).
1019 _divbyzero "coth($z)", "sinh($z)" if ($sz == 0);
1020 return cosh($z) / $sz;
1028 sub cotanh { Math::Complex::coth(@_) }
1033 # Computes the arc hyperbolic cosine acosh(z) = log(z + sqrt(z*z-1)).
1038 return CORE::log($z + CORE::sqrt($z*$z-1)) if $z >= 1;
1041 my ($re, $im) = @{$z->cartesian};
1043 return cplx(CORE::log($re + CORE::sqrt($re*$re - 1)), 0) if $re >= 1;
1044 return cplx(0, CORE::atan2(CORE::sqrt(1-$re*$re), $re)) if CORE::abs($re) <= 1;
1046 return CORE::log($z + CORE::sqrt($z*$z - 1));
1052 # Computes the arc hyperbolic sine asinh(z) = log(z + sqrt(z*z-1))
1056 return CORE::log($z + CORE::sqrt($z*$z + 1));
1062 # Computes the arc hyperbolic tangent atanh(z) = 1/2 log((1+z) / (1-z)).
1067 return CORE::log((1 + $z)/(1 - $z))/2 if CORE::abs($z) < 1;
1070 _divbyzero 'atanh(1)', "1 - $z" if ($z == 1);
1071 _logofzero 'atanh(-1)' if ($z == -1);
1072 return 0.5 * CORE::log((1 + $z) / (1 - $z));
1078 # Computes the hyperbolic arc secant asech(z) = acosh(1 / z).
1082 _divbyzero 'asech(0)', $z if ($z == 0);
1083 return acosh(1 / $z);
1089 # Computes the hyperbolic arc cosecant acsch(z) = asinh(1 / z).
1093 _divbyzero 'acsch(0)', $z if ($z == 0);
1094 return asinh(1 / $z);
1100 # Alias for acosh().
1102 sub acosech { Math::Complex::acsch(@_) }
1107 # Computes the arc hyperbolic cotangent acoth(z) = 1/2 log((1+z) / (z-1)).
1111 _divbyzero 'acoth(0)' if (CORE::abs($z) < $eps);
1113 return CORE::log(($z + 1)/($z - 1))/2 if CORE::abs($z) > 1;
1116 _divbyzero 'acoth(1)', "$z - 1" if (CORE::abs($z - 1) < $eps);
1117 _logofzero 'acoth(-1)', "1 / $z" if (CORE::abs($z + 1) < $eps);
1118 return CORE::log((1 + $z) / ($z - 1)) / 2;
1126 sub acotanh { Math::Complex::acoth(@_) }
1131 # Compute atan(z1/z2).
1134 my ($z1, $z2, $inverted) = @_;
1135 my ($re1, $im1, $re2, $im2);
1137 ($re1, $im1) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
1138 ($re2, $im2) = @{$z1->cartesian};
1140 ($re1, $im1) = @{$z1->cartesian};
1141 ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
1144 return cplx(CORE::atan2($re1, $re2), 0) if $im1 == 0;
1145 return cplx(($im1<=>0) * pip2, 0) if $re2 == 0;
1147 my $w = atan($z1/$z2);
1148 my ($u, $v) = ref $w ? @{$w->cartesian} : ($w, 0);
1149 $u += pi if $re2 < 0;
1150 $u -= pit2 if $u > pi;
1151 return cplx($u, $v);
1158 # Set (get if no argument) the display format for all complex numbers that
1159 # don't happen to have overridden it via ->display_format
1161 # When called as an object method, this actually sets the display format for
1162 # the current object.
1164 # Valid object formats are 'c' and 'p' for cartesian and polar. The first
1165 # letter is used actually, so the type can be fully spelled out for clarity.
1167 sub display_format {
1169 my %display_format = %DISPLAY_FORMAT;
1171 if (ref $self) { # Called as an object method
1172 if (exists $self->{display_format}) {
1173 my %obj = %{$self->{display_format}};
1174 @display_format{keys %obj} = values %obj;
1177 $display_format{style} = shift;
1180 @display_format{keys %new} = values %new;
1182 } else { # Called as a class method
1184 $display_format{style} = $self;
1187 @display_format{keys %new} = values %new;
1192 if (defined $self) {
1193 $self->{display_format} = { %display_format };
1196 %{$self->{display_format}} :
1197 $self->{display_format}->{style};
1200 %DISPLAY_FORMAT = %display_format;
1204 $DISPLAY_FORMAT{style};
1210 # Show nicely formatted complex number under its cartesian or polar form,
1211 # depending on the current display format:
1213 # . If a specific display format has been recorded for this object, use it.
1214 # . Otherwise, use the generic current default for all complex numbers,
1215 # which is a package global variable.
1220 my $style = $z->display_format;
1222 $style = $DISPLAY_FORMAT{style} unless defined $style;
1224 return $z->stringify_polar if $style =~ /^p/i;
1225 return $z->stringify_cartesian;
1229 # ->stringify_cartesian
1231 # Stringify as a cartesian representation 'a+bi'.
1233 sub stringify_cartesian {
1235 my ($x, $y) = @{$z->cartesian};
1238 $x = int($x + ($x < 0 ? -1 : 1) * $eps)
1239 if int(CORE::abs($x)) != int(CORE::abs($x) + $eps);
1240 $y = int($y + ($y < 0 ? -1 : 1) * $eps)
1241 if int(CORE::abs($y)) != int(CORE::abs($y) + $eps);
1243 $re = "$x" if CORE::abs($x) >= $eps;
1245 my %format = $z->display_format;
1246 my $format = $format{format};
1248 if ($y == 1) { $im = 'i' }
1249 elsif ($y == -1) { $im = '-i' }
1250 elsif (CORE::abs($y) >= $eps) {
1251 $im = (defined $format ? sprintf($format, $y) : $y) . "i";
1255 $str = defined $format ? sprintf($format, $re) : $re
1261 $str .= "+" if defined $re;
1270 # Helper for stringify_polar, a Greatest Common Divisor with a memory.
1277 # Loops forever if given negative inputs.
1279 if ($b and $a > $b) { return gcd($a % $b, $b) }
1280 elsif ($a and $b > $a) { return gcd($b % $a, $a) }
1281 else { return $a ? $a : $b }
1291 unless (exists $gcd{$id}) {
1292 $gcd{$id} = _gcd($a, $b);
1293 $gcd{"$b $a"} = $gcd{$id};
1302 # Stringify as a polar representation '[r,t]'.
1304 sub stringify_polar {
1306 my ($r, $t) = @{$z->polar};
1309 return '[0,0]' if $r <= $eps;
1311 my %format = $z->display_format;
1314 $nt = ($nt - int($nt)) * pit2;
1315 $nt += pit2 if $nt < 0; # Range [0, 2pi]
1317 if (CORE::abs($nt) <= $eps) { $theta = 0 }
1318 elsif (CORE::abs(pi-$nt) <= $eps) { $theta = 'pi' }
1320 if (defined $theta) {
1321 $r = int($r + ($r < 0 ? -1 : 1) * $eps)
1322 if int(CORE::abs($r)) != int(CORE::abs($r) + $eps);
1323 $theta = int($theta + ($theta < 0 ? -1 : 1) * $eps)
1324 if ($theta ne 'pi' and
1325 int(CORE::abs($theta)) != int(CORE::abs($theta) + $eps));
1326 return "\[$r,$theta\]";
1330 # Okay, number is not a real. Try to identify pi/n and friends...
1333 $nt -= pit2 if $nt > pi;
1335 if ($format{polar_pretty_print} && CORE::abs($nt) >= deg1) {
1338 for ($k = 1, $kpi = pi; $k < 10; $k++, $kpi += pi) {
1339 $n = int($kpi / $nt + ($nt > 0 ? 1 : -1) * 0.5);
1340 if (CORE::abs($kpi/$n - $nt) <= $eps) {
1342 my $gcd = gcd($k, $n);
1348 $theta = ($nt < 0 ? '-':'').
1349 ($k == 1 ? 'pi':"${k}pi");
1350 $theta .= '/'.$n if $n > 1;
1356 $theta = $nt unless defined $theta;
1358 $r = int($r + ($r < 0 ? -1 : 1) * $eps)
1359 if int(CORE::abs($r)) != int(CORE::abs($r) + $eps);
1360 $theta = int($theta + ($theta < 0 ? -1 : 1) * $eps)
1361 if ($theta !~ m(^-?\d*pi/\d+$) and
1362 int(CORE::abs($theta)) != int(CORE::abs($theta) + $eps));
1364 my $format = $format{format};
1365 if (defined $format) {
1366 $r = sprintf($format, $r);
1367 $theta = sprintf($format, $theta);
1370 return "\[$r,$theta\]";
1378 Math::Complex - complex numbers and associated mathematical functions
1384 $z = Math::Complex->make(5, 6);
1386 $j = cplxe(1, 2*pi/3);
1390 This package lets you create and manipulate complex numbers. By default,
1391 I<Perl> limits itself to real numbers, but an extra C<use> statement brings
1392 full complex support, along with a full set of mathematical functions
1393 typically associated with and/or extended to complex numbers.
1395 If you wonder what complex numbers are, they were invented to be able to solve
1396 the following equation:
1400 and by definition, the solution is noted I<i> (engineers use I<j> instead since
1401 I<i> usually denotes an intensity, but the name does not matter). The number
1402 I<i> is a pure I<imaginary> number.
1404 The arithmetics with pure imaginary numbers works just like you would expect
1405 it with real numbers... you just have to remember that
1411 5i + 7i = i * (5 + 7) = 12i
1412 4i - 3i = i * (4 - 3) = i
1417 Complex numbers are numbers that have both a real part and an imaginary
1418 part, and are usually noted:
1422 where C<a> is the I<real> part and C<b> is the I<imaginary> part. The
1423 arithmetic with complex numbers is straightforward. You have to
1424 keep track of the real and the imaginary parts, but otherwise the
1425 rules used for real numbers just apply:
1427 (4 + 3i) + (5 - 2i) = (4 + 5) + i(3 - 2) = 9 + i
1428 (2 + i) * (4 - i) = 2*4 + 4i -2i -i*i = 8 + 2i + 1 = 9 + 2i
1430 A graphical representation of complex numbers is possible in a plane
1431 (also called the I<complex plane>, but it's really a 2D plane).
1436 is the point whose coordinates are (a, b). Actually, it would
1437 be the vector originating from (0, 0) to (a, b). It follows that the addition
1438 of two complex numbers is a vectorial addition.
1440 Since there is a bijection between a point in the 2D plane and a complex
1441 number (i.e. the mapping is unique and reciprocal), a complex number
1442 can also be uniquely identified with polar coordinates:
1446 where C<rho> is the distance to the origin, and C<theta> the angle between
1447 the vector and the I<x> axis. There is a notation for this using the
1448 exponential form, which is:
1450 rho * exp(i * theta)
1452 where I<i> is the famous imaginary number introduced above. Conversion
1453 between this form and the cartesian form C<a + bi> is immediate:
1455 a = rho * cos(theta)
1456 b = rho * sin(theta)
1458 which is also expressed by this formula:
1460 z = rho * exp(i * theta) = rho * (cos theta + i * sin theta)
1462 In other words, it's the projection of the vector onto the I<x> and I<y>
1463 axes. Mathematicians call I<rho> the I<norm> or I<modulus> and I<theta>
1464 the I<argument> of the complex number. The I<norm> of C<z> will be
1467 The polar notation (also known as the trigonometric
1468 representation) is much more handy for performing multiplications and
1469 divisions of complex numbers, whilst the cartesian notation is better
1470 suited for additions and subtractions. Real numbers are on the I<x>
1471 axis, and therefore I<theta> is zero or I<pi>.
1473 All the common operations that can be performed on a real number have
1474 been defined to work on complex numbers as well, and are merely
1475 I<extensions> of the operations defined on real numbers. This means
1476 they keep their natural meaning when there is no imaginary part, provided
1477 the number is within their definition set.
1479 For instance, the C<sqrt> routine which computes the square root of
1480 its argument is only defined for non-negative real numbers and yields a
1481 non-negative real number (it is an application from B<R+> to B<R+>).
1482 If we allow it to return a complex number, then it can be extended to
1483 negative real numbers to become an application from B<R> to B<C> (the
1484 set of complex numbers):
1486 sqrt(x) = x >= 0 ? sqrt(x) : sqrt(-x)*i
1488 It can also be extended to be an application from B<C> to B<C>,
1489 whilst its restriction to B<R> behaves as defined above by using
1490 the following definition:
1492 sqrt(z = [r,t]) = sqrt(r) * exp(i * t/2)
1494 Indeed, a negative real number can be noted C<[x,pi]> (the modulus
1495 I<x> is always non-negative, so C<[x,pi]> is really C<-x>, a negative
1496 number) and the above definition states that
1498 sqrt([x,pi]) = sqrt(x) * exp(i*pi/2) = [sqrt(x),pi/2] = sqrt(x)*i
1500 which is exactly what we had defined for negative real numbers above.
1501 The C<sqrt> returns only one of the solutions: if you want the both,
1502 use the C<root> function.
1504 All the common mathematical functions defined on real numbers that
1505 are extended to complex numbers share that same property of working
1506 I<as usual> when the imaginary part is zero (otherwise, it would not
1507 be called an extension, would it?).
1509 A I<new> operation possible on a complex number that is
1510 the identity for real numbers is called the I<conjugate>, and is noted
1511 with an horizontal bar above the number, or C<~z> here.
1518 z * ~z = (a + bi) * (a - bi) = a*a + b*b
1520 We saw that the norm of C<z> was noted C<abs(z)> and was defined as the
1521 distance to the origin, also known as:
1523 rho = abs(z) = sqrt(a*a + b*b)
1527 z * ~z = abs(z) ** 2
1529 If z is a pure real number (i.e. C<b == 0>), then the above yields:
1533 which is true (C<abs> has the regular meaning for real number, i.e. stands
1534 for the absolute value). This example explains why the norm of C<z> is
1535 noted C<abs(z)>: it extends the C<abs> function to complex numbers, yet
1536 is the regular C<abs> we know when the complex number actually has no
1537 imaginary part... This justifies I<a posteriori> our use of the C<abs>
1538 notation for the norm.
1542 Given the following notations:
1544 z1 = a + bi = r1 * exp(i * t1)
1545 z2 = c + di = r2 * exp(i * t2)
1546 z = <any complex or real number>
1548 the following (overloaded) operations are supported on complex numbers:
1550 z1 + z2 = (a + c) + i(b + d)
1551 z1 - z2 = (a - c) + i(b - d)
1552 z1 * z2 = (r1 * r2) * exp(i * (t1 + t2))
1553 z1 / z2 = (r1 / r2) * exp(i * (t1 - t2))
1554 z1 ** z2 = exp(z2 * log z1)
1556 abs(z) = r1 = sqrt(a*a + b*b)
1557 sqrt(z) = sqrt(r1) * exp(i * t/2)
1558 exp(z) = exp(a) * exp(i * b)
1559 log(z) = log(r1) + i*t
1560 sin(z) = 1/2i (exp(i * z1) - exp(-i * z))
1561 cos(z) = 1/2 (exp(i * z1) + exp(-i * z))
1562 atan2(z1, z2) = atan(z1/z2)
1564 The following extra operations are supported on both real and complex
1572 cbrt(z) = z ** (1/3)
1573 log10(z) = log(z) / log(10)
1574 logn(z, n) = log(z) / log(n)
1576 tan(z) = sin(z) / cos(z)
1582 asin(z) = -i * log(i*z + sqrt(1-z*z))
1583 acos(z) = -i * log(z + i*sqrt(1-z*z))
1584 atan(z) = i/2 * log((i+z) / (i-z))
1586 acsc(z) = asin(1 / z)
1587 asec(z) = acos(1 / z)
1588 acot(z) = atan(1 / z) = -i/2 * log((i+z) / (z-i))
1590 sinh(z) = 1/2 (exp(z) - exp(-z))
1591 cosh(z) = 1/2 (exp(z) + exp(-z))
1592 tanh(z) = sinh(z) / cosh(z) = (exp(z) - exp(-z)) / (exp(z) + exp(-z))
1594 csch(z) = 1 / sinh(z)
1595 sech(z) = 1 / cosh(z)
1596 coth(z) = 1 / tanh(z)
1598 asinh(z) = log(z + sqrt(z*z+1))
1599 acosh(z) = log(z + sqrt(z*z-1))
1600 atanh(z) = 1/2 * log((1+z) / (1-z))
1602 acsch(z) = asinh(1 / z)
1603 asech(z) = acosh(1 / z)
1604 acoth(z) = atanh(1 / z) = 1/2 * log((1+z) / (z-1))
1606 I<arg>, I<abs>, I<log>, I<csc>, I<cot>, I<acsc>, I<acot>, I<csch>,
1607 I<coth>, I<acosech>, I<acotanh>, have aliases I<rho>, I<theta>, I<ln>,
1608 I<cosec>, I<cotan>, I<acosec>, I<acotan>, I<cosech>, I<cotanh>,
1609 I<acosech>, I<acotanh>, respectively. C<Re>, C<Im>, C<arg>, C<abs>,
1610 C<rho>, and C<theta> can be used also also mutators. The C<cbrt>
1611 returns only one of the solutions: if you want all three, use the
1614 The I<root> function is available to compute all the I<n>
1615 roots of some complex, where I<n> is a strictly positive integer.
1616 There are exactly I<n> such roots, returned as a list. Getting the
1617 number mathematicians call C<j> such that:
1621 is a simple matter of writing:
1623 $j = ((root(1, 3))[1];
1625 The I<k>th root for C<z = [r,t]> is given by:
1627 (root(z, n))[k] = r**(1/n) * exp(i * (t + 2*k*pi)/n)
1629 The I<spaceship> comparison operator, E<lt>=E<gt>, is also defined. In
1630 order to ensure its restriction to real numbers is conform to what you
1631 would expect, the comparison is run on the real part of the complex
1632 number first, and imaginary parts are compared only when the real
1637 To create a complex number, use either:
1639 $z = Math::Complex->make(3, 4);
1642 if you know the cartesian form of the number, or
1646 if you like. To create a number using the polar form, use either:
1648 $z = Math::Complex->emake(5, pi/3);
1649 $x = cplxe(5, pi/3);
1651 instead. The first argument is the modulus, the second is the angle
1652 (in radians, the full circle is 2*pi). (Mnemonic: C<e> is used as a
1653 notation for complex numbers in the polar form).
1655 It is possible to write:
1657 $x = cplxe(-3, pi/4);
1659 but that will be silently converted into C<[3,-3pi/4]>, since the
1660 modulus must be non-negative (it represents the distance to the origin
1661 in the complex plane).
1663 It is also possible to have a complex number as either argument of
1664 either the C<make> or C<emake>: the appropriate component of
1665 the argument will be used.
1670 =head1 STRINGIFICATION
1672 When printed, a complex number is usually shown under its cartesian
1673 style I<a+bi>, but there are legitimate cases where the polar style
1674 I<[r,t]> is more appropriate.
1676 In the polar style Math::Complex will try to recognize certain common
1677 numbers such as multiples or small rationals of pi (2pi, pi/2) and
1678 prettyprint those numbers.
1680 By calling the class method C<Math::Complex::display_format> and
1681 supplying either C<"polar"> or C<"cartesian"> as an argument, you
1682 override the default display format, which is C<"cartesian">. Not
1683 supplying any argument returns the current settings.
1685 This default can be overridden on a per-number basis by calling the
1686 C<display_format> method instead. As before, not supplying any argument
1687 returns the current display format for this number. Otherwise whatever you
1688 specify will be the new display format for I<this> particular number.
1694 Math::Complex::display_format('polar');
1695 $j = (root(1, 3))[1];
1696 print "j = $j\n"; # Prints "j = [1,2pi/3]"
1697 $j->display_format('cartesian');
1698 print "j = $j\n"; # Prints "j = -0.5+0.866025403784439i"
1700 The polar format attempts to emphasize arguments like I<k*pi/n>
1701 (where I<n> is a positive integer and I<k> an integer within [-9,+9]).
1703 =head2 CHANGED IN PERL 5.6
1705 The C<display_format> class method and the corresponding
1706 C<display_format> object method can now be called using
1707 a parameter hash instead of just a one parameter.
1709 The old display format style, which can have values C<"cartesian"> or
1710 C<"polar">, can be changed using the C<"style"> parameter. (The one
1711 parameter calling convention also still works.)
1713 There are two new display parameters.
1715 The first one is C<"format">, which is a sprintf()-style format
1716 string to be used for both parts of the complex number(s). The
1717 default is C<undef>, which corresponds usually (this is somewhat
1718 system-dependent) to C<"%.15g">. You can revert to the default by
1719 setting the format string to C<undef>.
1721 # the $j from the above example
1723 $j->display_format('format' => '%.5f');
1724 print "j = $j\n"; # Prints "j = -0.50000+0.86603i"
1725 $j->display_format('format' => '%.6f');
1726 print "j = $j\n"; # Prints "j = -0.5+0.86603i"
1728 Notice that this affects also the return values of the
1729 C<display_format> methods: in list context the whole parameter hash
1730 will be returned, as opposed to only the style parameter value. If
1731 you want to know the whole truth for a complex number, you must call
1732 both the class method and the object method:
1734 The second new display parameter is C<"polar_pretty_print">, which can be
1735 set to true or false, the default being true. See above for what this
1740 Thanks to overloading, the handling of arithmetics with complex numbers
1741 is simple and almost transparent.
1743 Here are some examples:
1747 $j = cplxe(1, 2*pi/3); # $j ** 3 == 1
1748 print "j = $j, j**3 = ", $j ** 3, "\n";
1749 print "1 + j + j**2 = ", 1 + $j + $j**2, "\n";
1751 $z = -16 + 0*i; # Force it to be a complex
1752 print "sqrt($z) = ", sqrt($z), "\n";
1754 $k = exp(i * 2*pi/3);
1755 print "$j - $k = ", $j - $k, "\n";
1757 $z->Re(3); # Re, Im, arg, abs,
1758 $j->arg(2); # (the last two aka rho, theta)
1759 # can be used also as mutators.
1761 =head1 ERRORS DUE TO DIVISION BY ZERO OR LOGARITHM OF ZERO
1763 The division (/) and the following functions
1769 atanh asech acsch acoth
1771 cannot be computed for all arguments because that would mean dividing
1772 by zero or taking logarithm of zero. These situations cause fatal
1773 runtime errors looking like this
1775 cot(0): Division by zero.
1776 (Because in the definition of cot(0), the divisor sin(0) is 0)
1781 atanh(-1): Logarithm of zero.
1784 For the C<csc>, C<cot>, C<asec>, C<acsc>, C<acot>, C<csch>, C<coth>,
1785 C<asech>, C<acsch>, the argument cannot be C<0> (zero). For the the
1786 logarithmic functions and the C<atanh>, C<acoth>, the argument cannot
1787 be C<1> (one). For the C<atanh>, C<acoth>, the argument cannot be
1788 C<-1> (minus one). For the C<atan>, C<acot>, the argument cannot be
1789 C<i> (the imaginary unit). For the C<atan>, C<acoth>, the argument
1790 cannot be C<-i> (the negative imaginary unit). For the C<tan>,
1791 C<sec>, C<tanh>, the argument cannot be I<pi/2 + k * pi>, where I<k>
1794 Note that because we are operating on approximations of real numbers,
1795 these errors can happen when merely `too close' to the singularities
1796 listed above. For example C<tan(2*atan2(1,1)+1e-15)> will die of
1799 =head1 ERRORS DUE TO INDIGESTIBLE ARGUMENTS
1801 The C<make> and C<emake> accept both real and complex arguments.
1802 When they cannot recognize the arguments they will die with error
1803 messages like the following
1805 Math::Complex::make: Cannot take real part of ...
1806 Math::Complex::make: Cannot take real part of ...
1807 Math::Complex::emake: Cannot take rho of ...
1808 Math::Complex::emake: Cannot take theta of ...
1812 Saying C<use Math::Complex;> exports many mathematical routines in the
1813 caller environment and even overrides some (C<sqrt>, C<log>).
1814 This is construed as a feature by the Authors, actually... ;-)
1816 All routines expect to be given real or complex numbers. Don't attempt to
1817 use BigFloat, since Perl has currently no rule to disambiguate a '+'
1818 operation (for instance) between two overloaded entities.
1820 In Cray UNICOS there is some strange numerical instability that results
1821 in root(), cos(), sin(), cosh(), sinh(), losing accuracy fast. Beware.
1822 The bug may be in UNICOS math libs, in UNICOS C compiler, in Math::Complex.
1823 Whatever it is, it does not manifest itself anywhere else where Perl runs.
1827 Raphael Manfredi <F<Raphael_Manfredi@pobox.com>> and
1828 Jarkko Hietaniemi <F<jhi@iki.fi>>.
1830 Extensive patches by Daniel S. Lewart <F<d-lewart@uiuc.edu>>.