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\]";
1379 Math::Complex - complex numbers and associated mathematical functions
1385 $z = Math::Complex->make(5, 6);
1387 $j = cplxe(1, 2*pi/3);
1391 This package lets you create and manipulate complex numbers. By default,
1392 I<Perl> limits itself to real numbers, but an extra C<use> statement brings
1393 full complex support, along with a full set of mathematical functions
1394 typically associated with and/or extended to complex numbers.
1396 If you wonder what complex numbers are, they were invented to be able to solve
1397 the following equation:
1401 and by definition, the solution is noted I<i> (engineers use I<j> instead since
1402 I<i> usually denotes an intensity, but the name does not matter). The number
1403 I<i> is a pure I<imaginary> number.
1405 The arithmetics with pure imaginary numbers works just like you would expect
1406 it with real numbers... you just have to remember that
1412 5i + 7i = i * (5 + 7) = 12i
1413 4i - 3i = i * (4 - 3) = i
1418 Complex numbers are numbers that have both a real part and an imaginary
1419 part, and are usually noted:
1423 where C<a> is the I<real> part and C<b> is the I<imaginary> part. The
1424 arithmetic with complex numbers is straightforward. You have to
1425 keep track of the real and the imaginary parts, but otherwise the
1426 rules used for real numbers just apply:
1428 (4 + 3i) + (5 - 2i) = (4 + 5) + i(3 - 2) = 9 + i
1429 (2 + i) * (4 - i) = 2*4 + 4i -2i -i*i = 8 + 2i + 1 = 9 + 2i
1431 A graphical representation of complex numbers is possible in a plane
1432 (also called the I<complex plane>, but it's really a 2D plane).
1437 is the point whose coordinates are (a, b). Actually, it would
1438 be the vector originating from (0, 0) to (a, b). It follows that the addition
1439 of two complex numbers is a vectorial addition.
1441 Since there is a bijection between a point in the 2D plane and a complex
1442 number (i.e. the mapping is unique and reciprocal), a complex number
1443 can also be uniquely identified with polar coordinates:
1447 where C<rho> is the distance to the origin, and C<theta> the angle between
1448 the vector and the I<x> axis. There is a notation for this using the
1449 exponential form, which is:
1451 rho * exp(i * theta)
1453 where I<i> is the famous imaginary number introduced above. Conversion
1454 between this form and the cartesian form C<a + bi> is immediate:
1456 a = rho * cos(theta)
1457 b = rho * sin(theta)
1459 which is also expressed by this formula:
1461 z = rho * exp(i * theta) = rho * (cos theta + i * sin theta)
1463 In other words, it's the projection of the vector onto the I<x> and I<y>
1464 axes. Mathematicians call I<rho> the I<norm> or I<modulus> and I<theta>
1465 the I<argument> of the complex number. The I<norm> of C<z> will be
1468 The polar notation (also known as the trigonometric
1469 representation) is much more handy for performing multiplications and
1470 divisions of complex numbers, whilst the cartesian notation is better
1471 suited for additions and subtractions. Real numbers are on the I<x>
1472 axis, and therefore I<theta> is zero or I<pi>.
1474 All the common operations that can be performed on a real number have
1475 been defined to work on complex numbers as well, and are merely
1476 I<extensions> of the operations defined on real numbers. This means
1477 they keep their natural meaning when there is no imaginary part, provided
1478 the number is within their definition set.
1480 For instance, the C<sqrt> routine which computes the square root of
1481 its argument is only defined for non-negative real numbers and yields a
1482 non-negative real number (it is an application from B<R+> to B<R+>).
1483 If we allow it to return a complex number, then it can be extended to
1484 negative real numbers to become an application from B<R> to B<C> (the
1485 set of complex numbers):
1487 sqrt(x) = x >= 0 ? sqrt(x) : sqrt(-x)*i
1489 It can also be extended to be an application from B<C> to B<C>,
1490 whilst its restriction to B<R> behaves as defined above by using
1491 the following definition:
1493 sqrt(z = [r,t]) = sqrt(r) * exp(i * t/2)
1495 Indeed, a negative real number can be noted C<[x,pi]> (the modulus
1496 I<x> is always non-negative, so C<[x,pi]> is really C<-x>, a negative
1497 number) and the above definition states that
1499 sqrt([x,pi]) = sqrt(x) * exp(i*pi/2) = [sqrt(x),pi/2] = sqrt(x)*i
1501 which is exactly what we had defined for negative real numbers above.
1502 The C<sqrt> returns only one of the solutions: if you want the both,
1503 use the C<root> function.
1505 All the common mathematical functions defined on real numbers that
1506 are extended to complex numbers share that same property of working
1507 I<as usual> when the imaginary part is zero (otherwise, it would not
1508 be called an extension, would it?).
1510 A I<new> operation possible on a complex number that is
1511 the identity for real numbers is called the I<conjugate>, and is noted
1512 with an horizontal bar above the number, or C<~z> here.
1519 z * ~z = (a + bi) * (a - bi) = a*a + b*b
1521 We saw that the norm of C<z> was noted C<abs(z)> and was defined as the
1522 distance to the origin, also known as:
1524 rho = abs(z) = sqrt(a*a + b*b)
1528 z * ~z = abs(z) ** 2
1530 If z is a pure real number (i.e. C<b == 0>), then the above yields:
1534 which is true (C<abs> has the regular meaning for real number, i.e. stands
1535 for the absolute value). This example explains why the norm of C<z> is
1536 noted C<abs(z)>: it extends the C<abs> function to complex numbers, yet
1537 is the regular C<abs> we know when the complex number actually has no
1538 imaginary part... This justifies I<a posteriori> our use of the C<abs>
1539 notation for the norm.
1543 Given the following notations:
1545 z1 = a + bi = r1 * exp(i * t1)
1546 z2 = c + di = r2 * exp(i * t2)
1547 z = <any complex or real number>
1549 the following (overloaded) operations are supported on complex numbers:
1551 z1 + z2 = (a + c) + i(b + d)
1552 z1 - z2 = (a - c) + i(b - d)
1553 z1 * z2 = (r1 * r2) * exp(i * (t1 + t2))
1554 z1 / z2 = (r1 / r2) * exp(i * (t1 - t2))
1555 z1 ** z2 = exp(z2 * log z1)
1557 abs(z) = r1 = sqrt(a*a + b*b)
1558 sqrt(z) = sqrt(r1) * exp(i * t/2)
1559 exp(z) = exp(a) * exp(i * b)
1560 log(z) = log(r1) + i*t
1561 sin(z) = 1/2i (exp(i * z1) - exp(-i * z))
1562 cos(z) = 1/2 (exp(i * z1) + exp(-i * z))
1563 atan2(z1, z2) = atan(z1/z2)
1565 The following extra operations are supported on both real and complex
1573 cbrt(z) = z ** (1/3)
1574 log10(z) = log(z) / log(10)
1575 logn(z, n) = log(z) / log(n)
1577 tan(z) = sin(z) / cos(z)
1583 asin(z) = -i * log(i*z + sqrt(1-z*z))
1584 acos(z) = -i * log(z + i*sqrt(1-z*z))
1585 atan(z) = i/2 * log((i+z) / (i-z))
1587 acsc(z) = asin(1 / z)
1588 asec(z) = acos(1 / z)
1589 acot(z) = atan(1 / z) = -i/2 * log((i+z) / (z-i))
1591 sinh(z) = 1/2 (exp(z) - exp(-z))
1592 cosh(z) = 1/2 (exp(z) + exp(-z))
1593 tanh(z) = sinh(z) / cosh(z) = (exp(z) - exp(-z)) / (exp(z) + exp(-z))
1595 csch(z) = 1 / sinh(z)
1596 sech(z) = 1 / cosh(z)
1597 coth(z) = 1 / tanh(z)
1599 asinh(z) = log(z + sqrt(z*z+1))
1600 acosh(z) = log(z + sqrt(z*z-1))
1601 atanh(z) = 1/2 * log((1+z) / (1-z))
1603 acsch(z) = asinh(1 / z)
1604 asech(z) = acosh(1 / z)
1605 acoth(z) = atanh(1 / z) = 1/2 * log((1+z) / (z-1))
1607 I<arg>, I<abs>, I<log>, I<csc>, I<cot>, I<acsc>, I<acot>, I<csch>,
1608 I<coth>, I<acosech>, I<acotanh>, have aliases I<rho>, I<theta>, I<ln>,
1609 I<cosec>, I<cotan>, I<acosec>, I<acotan>, I<cosech>, I<cotanh>,
1610 I<acosech>, I<acotanh>, respectively. C<Re>, C<Im>, C<arg>, C<abs>,
1611 C<rho>, and C<theta> can be used also also mutators. The C<cbrt>
1612 returns only one of the solutions: if you want all three, use the
1615 The I<root> function is available to compute all the I<n>
1616 roots of some complex, where I<n> is a strictly positive integer.
1617 There are exactly I<n> such roots, returned as a list. Getting the
1618 number mathematicians call C<j> such that:
1622 is a simple matter of writing:
1624 $j = ((root(1, 3))[1];
1626 The I<k>th root for C<z = [r,t]> is given by:
1628 (root(z, n))[k] = r**(1/n) * exp(i * (t + 2*k*pi)/n)
1630 The I<spaceship> comparison operator, E<lt>=E<gt>, is also defined. In
1631 order to ensure its restriction to real numbers is conform to what you
1632 would expect, the comparison is run on the real part of the complex
1633 number first, and imaginary parts are compared only when the real
1638 To create a complex number, use either:
1640 $z = Math::Complex->make(3, 4);
1643 if you know the cartesian form of the number, or
1647 if you like. To create a number using the polar form, use either:
1649 $z = Math::Complex->emake(5, pi/3);
1650 $x = cplxe(5, pi/3);
1652 instead. The first argument is the modulus, the second is the angle
1653 (in radians, the full circle is 2*pi). (Mnemonic: C<e> is used as a
1654 notation for complex numbers in the polar form).
1656 It is possible to write:
1658 $x = cplxe(-3, pi/4);
1660 but that will be silently converted into C<[3,-3pi/4]>, since the
1661 modulus must be non-negative (it represents the distance to the origin
1662 in the complex plane).
1664 It is also possible to have a complex number as either argument of
1665 either the C<make> or C<emake>: the appropriate component of
1666 the argument will be used.
1671 =head1 STRINGIFICATION
1673 When printed, a complex number is usually shown under its cartesian
1674 style I<a+bi>, but there are legitimate cases where the polar style
1675 I<[r,t]> is more appropriate.
1677 By calling the class method C<Math::Complex::display_format> and
1678 supplying either C<"polar"> or C<"cartesian"> as an argument, you
1679 override the default display style, which is C<"cartesian">. Not
1680 supplying any argument returns the current settings.
1682 This default can be overridden on a per-number basis by calling the
1683 C<display_format> method instead. As before, not supplying any argument
1684 returns the current display style for this number. Otherwise whatever you
1685 specify will be the new display style for I<this> particular number.
1691 Math::Complex::display_format('polar');
1692 $j = (root(1, 3))[1];
1693 print "j = $j\n"; # Prints "j = [1,2pi/3]"
1694 $j->display_format('cartesian');
1695 print "j = $j\n"; # Prints "j = -0.5+0.866025403784439i"
1697 The polar style attempts to emphasize arguments like I<k*pi/n>
1698 (where I<n> is a positive integer and I<k> an integer within [-9,+9]),
1699 this is called I<polar pretty-printing>.
1701 =head2 CHANGED IN PERL 5.6
1703 The C<display_format> class method and the corresponding
1704 C<display_format> object method can now be called using
1705 a parameter hash instead of just a one parameter.
1707 The old display format style, which can have values C<"cartesian"> or
1708 C<"polar">, can be changed using the C<"style"> parameter. (The one
1709 parameter calling convention also still works.)
1711 There are two new display parameters.
1713 The first one is C<"format">, which is a sprintf()-style format
1714 string to be used for both parts of the complex number(s). The
1715 default is C<undef>, which corresponds usually (this is somewhat
1716 system-dependent) to C<"%.15g">. You can revert to the default by
1717 setting the format string to C<undef>.
1719 # the $j from the above example
1721 $j->display_format('format' => '%.5f');
1722 print "j = $j\n"; # Prints "j = -0.50000+0.86603i"
1723 $j->display_format('format' => '%.6f');
1724 print "j = $j\n"; # Prints "j = -0.5+0.86603i"
1726 Notice that this affects also the return values of the
1727 C<display_format> methods: in list context the whole parameter hash
1728 will be returned, as opposed to only the style parameter value. If
1729 you want to know the whole truth for a complex number, you must call
1730 both the class method and the object method:
1732 The second new display parameter is C<"polar_pretty_print">, which can
1733 be set to true or false, the default being true. See the previous
1734 section for what this means.
1738 Thanks to overloading, the handling of arithmetics with complex numbers
1739 is simple and almost transparent.
1741 Here are some examples:
1745 $j = cplxe(1, 2*pi/3); # $j ** 3 == 1
1746 print "j = $j, j**3 = ", $j ** 3, "\n";
1747 print "1 + j + j**2 = ", 1 + $j + $j**2, "\n";
1749 $z = -16 + 0*i; # Force it to be a complex
1750 print "sqrt($z) = ", sqrt($z), "\n";
1752 $k = exp(i * 2*pi/3);
1753 print "$j - $k = ", $j - $k, "\n";
1755 $z->Re(3); # Re, Im, arg, abs,
1756 $j->arg(2); # (the last two aka rho, theta)
1757 # can be used also as mutators.
1759 =head1 ERRORS DUE TO DIVISION BY ZERO OR LOGARITHM OF ZERO
1761 The division (/) and the following functions
1767 atanh asech acsch acoth
1769 cannot be computed for all arguments because that would mean dividing
1770 by zero or taking logarithm of zero. These situations cause fatal
1771 runtime errors looking like this
1773 cot(0): Division by zero.
1774 (Because in the definition of cot(0), the divisor sin(0) is 0)
1779 atanh(-1): Logarithm of zero.
1782 For the C<csc>, C<cot>, C<asec>, C<acsc>, C<acot>, C<csch>, C<coth>,
1783 C<asech>, C<acsch>, the argument cannot be C<0> (zero). For the the
1784 logarithmic functions and the C<atanh>, C<acoth>, the argument cannot
1785 be C<1> (one). For the C<atanh>, C<acoth>, the argument cannot be
1786 C<-1> (minus one). For the C<atan>, C<acot>, the argument cannot be
1787 C<i> (the imaginary unit). For the C<atan>, C<acoth>, the argument
1788 cannot be C<-i> (the negative imaginary unit). For the C<tan>,
1789 C<sec>, C<tanh>, the argument cannot be I<pi/2 + k * pi>, where I<k>
1792 Note that because we are operating on approximations of real numbers,
1793 these errors can happen when merely `too close' to the singularities
1794 listed above. For example C<tan(2*atan2(1,1)+1e-15)> will die of
1797 =head1 ERRORS DUE TO INDIGESTIBLE ARGUMENTS
1799 The C<make> and C<emake> accept both real and complex arguments.
1800 When they cannot recognize the arguments they will die with error
1801 messages like the following
1803 Math::Complex::make: Cannot take real part of ...
1804 Math::Complex::make: Cannot take real part of ...
1805 Math::Complex::emake: Cannot take rho of ...
1806 Math::Complex::emake: Cannot take theta of ...
1810 Saying C<use Math::Complex;> exports many mathematical routines in the
1811 caller environment and even overrides some (C<sqrt>, C<log>).
1812 This is construed as a feature by the Authors, actually... ;-)
1814 All routines expect to be given real or complex numbers. Don't attempt to
1815 use BigFloat, since Perl has currently no rule to disambiguate a '+'
1816 operation (for instance) between two overloaded entities.
1818 In Cray UNICOS there is some strange numerical instability that results
1819 in root(), cos(), sin(), cosh(), sinh(), losing accuracy fast. Beware.
1820 The bug may be in UNICOS math libs, in UNICOS C compiler, in Math::Complex.
1821 Whatever it is, it does not manifest itself anywhere else where Perl runs.
1825 Raphael Manfredi <F<Raphael_Manfredi@pobox.com>> and
1826 Jarkko Hietaniemi <F<jhi@iki.fi>>.
1828 Extensive patches by Daniel S. Lewart <F<d-lewart@uiuc.edu>>.