3 # Complex numbers and associated mathematical functions
4 # -- Raphael Manfredi, September 1996
5 # -- Jarkko Hietaniemi, March-April 1997
12 use vars qw($VERSION @ISA
24 csc cosec sec cot cotan
26 acsc acosec asec acot acotan
28 csch cosech sech coth cotanh
30 acsch acosech asech acoth acotanh
68 $package = 'Math::Complex'; # Package name
69 $display = 'cartesian'; # Default display format
72 # Object attributes (internal):
73 # cartesian [real, imaginary] -- cartesian form
74 # polar [rho, theta] -- polar form
75 # c_dirty cartesian form not up-to-date
76 # p_dirty polar form not up-to-date
77 # display display format (package's global when not set)
83 # Create a new complex number (cartesian form)
86 my $self = bless {}, shift;
88 $self->{'cartesian'} = [$re, $im];
97 # Create a new complex number (exponential form)
100 my $self = bless {}, shift;
101 my ($rho, $theta) = @_;
102 $theta += pi() if $rho < 0;
103 $self->{'polar'} = [abs($rho), $theta];
104 $self->{p_dirty} = 0;
105 $self->{c_dirty} = 1;
109 sub new { &make } # For backward compatibility only.
114 # Creates a complex number from a (re, im) tuple.
115 # This avoids the burden of writing Math::Complex->make(re, im).
119 return $package->make($re, defined $im ? $im : 0);
125 # Creates a complex number from a (rho, theta) tuple.
126 # This avoids the burden of writing Math::Complex->emake(rho, theta).
129 my ($rho, $theta) = @_;
130 return $package->emake($rho, defined $theta ? $theta : 0);
136 # The number defined as 2 * pi = 360 degrees
139 use constant pi => 4 * atan2(1, 1);
147 use constant log10inv => 1 / log(10);
152 # The number defined as i*i = -1;
157 $i->{'cartesian'} = [0, 1];
158 $i->{'polar'} = [1, pi/2];
165 # Attribute access/set routines
168 sub cartesian {$_[0]->{c_dirty} ?
169 $_[0]->update_cartesian : $_[0]->{'cartesian'}}
170 sub polar {$_[0]->{p_dirty} ?
171 $_[0]->update_polar : $_[0]->{'polar'}}
173 sub set_cartesian { $_[0]->{p_dirty}++; $_[0]->{'cartesian'} = $_[1] }
174 sub set_polar { $_[0]->{c_dirty}++; $_[0]->{'polar'} = $_[1] }
179 # Recompute and return the cartesian form, given accurate polar form.
181 sub update_cartesian {
183 my ($r, $t) = @{$self->{'polar'}};
184 $self->{c_dirty} = 0;
185 return $self->{'cartesian'} = [$r * cos $t, $r * sin $t];
192 # Recompute and return the polar form, given accurate cartesian form.
196 my ($x, $y) = @{$self->{'cartesian'}};
197 $self->{p_dirty} = 0;
198 return $self->{'polar'} = [0, 0] if $x == 0 && $y == 0;
199 return $self->{'polar'} = [sqrt($x*$x + $y*$y), atan2($y, $x)];
208 my ($z1, $z2, $regular) = @_;
209 my ($re1, $im1) = @{$z1->cartesian};
210 my ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
211 unless (defined $regular) {
212 $z1->set_cartesian([$re1 + $re2, $im1 + $im2]);
215 return (ref $z1)->make($re1 + $re2, $im1 + $im2);
224 my ($z1, $z2, $inverted) = @_;
225 my ($re1, $im1) = @{$z1->cartesian};
226 my ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
227 unless (defined $inverted) {
228 $z1->set_cartesian([$re1 - $re2, $im1 - $im2]);
232 (ref $z1)->make($re2 - $re1, $im2 - $im1) :
233 (ref $z1)->make($re1 - $re2, $im1 - $im2);
242 my ($z1, $z2, $regular) = @_;
243 my ($r1, $t1) = @{$z1->polar};
244 my ($r2, $t2) = ref $z2 ?
245 @{$z2->polar} : (abs($z2), $z2 >= 0 ? 0 : pi);
246 unless (defined $regular) {
247 $z1->set_polar([$r1 * $r2, $t1 + $t2]);
250 return (ref $z1)->emake($r1 * $r2, $t1 + $t2);
256 # Die on division by zero.
259 my $mess = "$_[0]: Division by zero.\n";
262 $mess .= "(Because in the definition of $_[0], the divisor ";
263 $mess .= "$_[1] " unless ($_[1] eq '0');
269 $mess .= "Died at $up[1] line $up[2].\n";
277 # Die on zero raised to the zeroth.
280 my $mess = "The zero raised to the zeroth power is not defined.\n";
284 $mess .= "Died at $up[1] line $up[2].\n";
295 my ($z1, $z2, $inverted) = @_;
296 my ($r1, $t1) = @{$z1->polar};
297 my ($r2, $t2) = ref $z2 ?
298 @{$z2->polar} : (abs($z2), $z2 >= 0 ? 0 : pi);
299 unless (defined $inverted) {
300 divbyzero "$z1/0" if ($r2 == 0);
301 $z1->set_polar([$r1 / $r2, $t1 - $t2]);
305 divbyzero "$z2/0" if ($r1 == 0);
306 return (ref $z1)->emake($r2 / $r1, $t2 - $t1);
308 divbyzero "$z1/0" if ($r2 == 0);
309 return (ref $z1)->emake($r1 / $r2, $t1 - $t2);
316 # Computes z1**z2 = exp(z2 * log z1)).
319 my ($z1, $z2, $inverted) = @_;
320 zerotozero if ($z1 == 0 and $z2 == 0);
321 return exp($z1 * log $z2) if defined $inverted && $inverted;
322 return exp($z2 * log $z1);
328 # Computes z1 <=> z2.
329 # Sorts on the real part first, then on the imaginary part. Thus 2-4i > 3+8i.
332 my ($z1, $z2, $inverted) = @_;
333 my ($re1, $im1) = ref $z1 ? @{$z1->cartesian} : ($z1, 0);
334 my ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
335 my $sgn = $inverted ? -1 : 1;
336 return $sgn * ($re1 <=> $re2) if $re1 != $re2;
337 return $sgn * ($im1 <=> $im2);
348 my ($r, $t) = @{$z->polar};
349 return (ref $z)->emake($r, pi + $t);
351 my ($re, $im) = @{$z->cartesian};
352 return (ref $z)->make(-$re, -$im);
358 # Compute complex's conjugate.
363 my ($r, $t) = @{$z->polar};
364 return (ref $z)->emake($r, -$t);
366 my ($re, $im) = @{$z->cartesian};
367 return (ref $z)->make($re, -$im);
373 # Compute complex's norm (rho).
377 return abs($z) unless ref $z;
378 my ($r, $t) = @{$z->polar};
385 # Compute complex's argument (theta).
389 return ($z < 0 ? pi : 0) unless ref $z;
390 my ($r, $t) = @{$z->polar};
401 $z = cplx($z, 0) unless ref $z;
402 my ($r, $t) = @{$z->polar};
403 return (ref $z)->emake(sqrt($r), $t/2);
409 # Compute cbrt(z) (cubic root).
413 return cplx($z, 0) ** (1/3) unless ref $z;
414 my ($r, $t) = @{$z->polar};
415 return (ref $z)->emake($r**(1/3), $t/3);
421 # Computes all nth root for z, returning an array whose size is n.
422 # `n' must be a positive integer.
424 # The roots are given by (for k = 0..n-1):
426 # z^(1/n) = r^(1/n) (cos ((t+2 k pi)/n) + i sin ((t+2 k pi)/n))
431 return undef unless $n > 0;
432 my ($r, $t) = ref $z ? @{$z->polar} : (abs($z), $z >= 0 ? 0 : pi);
435 my $theta_inc = 2 * pi / $n;
436 my $rho = $r ** (1/$n);
438 my $complex = ref($z) || $package;
439 for ($k = 0, $theta = $t / $n; $k < $n; $k++, $theta += $theta_inc) {
440 push(@root, $complex->emake($rho, $theta));
452 return $z unless ref $z;
453 my ($re, $im) = @{$z->cartesian};
464 return 0 unless ref $z;
465 my ($re, $im) = @{$z->cartesian};
476 $z = cplx($z, 0) unless ref $z;
477 my ($x, $y) = @{$z->cartesian};
478 return (ref $z)->emake(exp($x), $y);
488 $z = cplx($z, 0) unless ref $z;
489 my ($x, $y) = @{$z->cartesian};
490 my ($r, $t) = @{$z->polar};
491 $t -= 2 * pi if ($t > pi() and $x < 0);
492 $t += 2 * pi if ($t < -pi() and $x < 0);
493 return (ref $z)->make(log($r), $t);
501 sub ln { Math::Complex::log(@_) }
512 return log(cplx($z, 0)) * log10inv unless ref $z;
513 my ($r, $t) = @{$z->polar};
514 return (ref $z)->make(log($r) * log10inv, $t * log10inv);
520 # Compute logn(z,n) = log(z) / log(n)
524 $z = cplx($z, 0) unless ref $z;
525 my $logn = $logn{$n};
526 $logn = $logn{$n} = log($n) unless defined $logn; # Cache log(n)
527 return log($z) / $logn;
533 # Compute cos(z) = (exp(iz) + exp(-iz))/2.
537 $z = cplx($z, 0) unless ref $z;
538 my ($x, $y) = @{$z->cartesian};
541 return (ref $z)->make(cos($x) * ($ey + $ey_1)/2,
542 sin($x) * ($ey_1 - $ey)/2);
548 # Compute sin(z) = (exp(iz) - exp(-iz))/2.
552 $z = cplx($z, 0) unless ref $z;
553 my ($x, $y) = @{$z->cartesian};
556 return (ref $z)->make(sin($x) * ($ey + $ey_1)/2,
557 cos($x) * ($ey - $ey_1)/2);
563 # Compute tan(z) = sin(z) / cos(z).
568 divbyzero "tan($z)", "cos($z)" if ($cz == 0);
569 return sin($z) / $cz;
575 # Computes the secant sec(z) = 1 / cos(z).
580 divbyzero "sec($z)", "cos($z)" if ($cz == 0);
587 # Computes the cosecant csc(z) = 1 / sin(z).
592 divbyzero "csc($z)", "sin($z)" if ($sz == 0);
601 sub cosec { Math::Complex::csc(@_) }
606 # Computes cot(z) = 1 / tan(z).
611 divbyzero "cot($z)", "sin($z)" if ($sz == 0);
612 return cos($z) / $sz;
620 sub cotan { Math::Complex::cot(@_) }
625 # Computes the arc cosine acos(z) = -i log(z + sqrt(z*z-1)).
629 $z = cplx($z, 0) unless ref $z;
630 return ~i * log($z + (Re($z) * Im($z) > 0 ? 1 : -1) * sqrt($z*$z - 1));
636 # Computes the arc sine asin(z) = -i log(iz + sqrt(1-z*z)).
640 $z = cplx($z, 0) unless ref $z;
641 return ~i * log(i * $z + sqrt(1 - $z*$z));
647 # Computes the arc tangent atan(z) = i/2 log((i+z) / (i-z)).
651 $z = cplx($z, 0) unless ref $z;
652 divbyzero "atan($z)", "i - $z" if ($z == i);
653 return i/2*log((i + $z) / (i - $z));
659 # Computes the arc secant asec(z) = acos(1 / z).
663 divbyzero "asec($z)", $z if ($z == 0);
670 # Computes the arc cosecant sec(z) = asin(1 / z).
674 divbyzero "acsc($z)", $z if ($z == 0);
683 sub acosec { Math::Complex::acsc(@_) }
688 # Computes the arc cotangent acot(z) = -i/2 log((i+z) / (z-i))
692 $z = cplx($z, 0) unless ref $z;
693 divbyzero "acot($z)", "$z - i" if ($z == i);
694 return i/-2 * log((i + $z) / ($z - i));
702 sub acotan { Math::Complex::acot(@_) }
707 # Computes the hyperbolic cosine cosh(z) = (exp(z) + exp(-z))/2.
711 my ($x, $y) = ref $z ? @{$z->cartesian} : ($z, 0);
714 return ($ex + $ex_1)/2 unless ref $z;
715 return (ref $z)->make(cos($y) * ($ex + $ex_1)/2,
716 sin($y) * ($ex - $ex_1)/2);
722 # Computes the hyperbolic sine sinh(z) = (exp(z) - exp(-z))/2.
726 my ($x, $y) = ref $z ? @{$z->cartesian} : ($z, 0);
729 return ($ex - $ex_1)/2 unless ref $z;
730 return (ref $z)->make(cos($y) * ($ex - $ex_1)/2,
731 sin($y) * ($ex + $ex_1)/2);
737 # Computes the hyperbolic tangent tanh(z) = sinh(z) / cosh(z).
742 divbyzero "tanh($z)", "cosh($z)" if ($cz == 0);
743 return sinh($z) / $cz;
749 # Computes the hyperbolic secant sech(z) = 1 / cosh(z).
754 divbyzero "sech($z)", "cosh($z)" if ($cz == 0);
761 # Computes the hyperbolic cosecant csch(z) = 1 / sinh(z).
766 divbyzero "csch($z)", "sinh($z)" if ($sz == 0);
775 sub cosech { Math::Complex::csch(@_) }
780 # Computes the hyperbolic cotangent coth(z) = cosh(z) / sinh(z).
785 divbyzero "coth($z)", "sinh($z)" if ($sz == 0);
786 return cosh($z) / $sz;
794 sub cotanh { Math::Complex::coth(@_) }
799 # Computes the arc hyperbolic cosine acosh(z) = log(z + sqrt(z*z-1)).
803 $z = cplx($z, 0) unless ref $z;
804 return log($z + sqrt($z*$z - 1));
810 # Computes the arc hyperbolic sine asinh(z) = log(z + sqrt(z*z-1))
814 $z = cplx($z, 0) unless ref $z;
815 return log($z + sqrt($z*$z + 1));
821 # Computes the arc hyperbolic tangent atanh(z) = 1/2 log((1+z) / (1-z)).
825 divbyzero 'atanh(1)', "1 - $z" if ($z == 1);
826 $z = cplx($z, 0) unless ref $z;
827 my $cz = (1 + $z) / (1 - $z);
834 # Computes the hyperbolic arc secant asech(z) = acosh(1 / z).
838 divbyzero 'asech(0)', $z if ($z == 0);
839 return acosh(1 / $z);
845 # Computes the hyperbolic arc cosecant acsch(z) = asinh(1 / z).
849 divbyzero 'acsch(0)', $z if ($z == 0);
850 return asinh(1 / $z);
858 sub acosech { Math::Complex::acsch(@_) }
863 # Computes the arc hyperbolic cotangent acoth(z) = 1/2 log((1+z) / (z-1)).
867 divbyzero 'acoth(1)', "$z - 1" if ($z == 1);
868 $z = cplx($z, 0) unless ref $z;
869 my $cz = (1 + $z) / ($z - 1);
878 sub acotanh { Math::Complex::acoth(@_) }
883 # Compute atan(z1/z2).
886 my ($z1, $z2, $inverted) = @_;
887 my ($re1, $im1) = ref $z1 ? @{$z1->cartesian} : ($z1, 0);
888 my ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
890 if (defined $inverted && $inverted) { # atan(z2/z1)
891 return pi * ($re2 > 0 ? 1 : -1) if $re1 == 0 && $im1 == 0;
894 return pi * ($re1 > 0 ? 1 : -1) if $re2 == 0 && $im2 == 0;
904 # Set (fetch if no argument) display format for all complex numbers that
905 # don't happen to have overrriden it via ->display_format
907 # When called as a method, this actually sets the display format for
908 # the current object.
910 # Valid object formats are 'c' and 'p' for cartesian and polar. The first
911 # letter is used actually, so the type can be fully spelled out for clarity.
917 if (ref $self) { # Called as a method
919 } else { # Regular procedure call
925 return defined $self->{display} ? $self->{display} : $display
926 unless defined $format;
927 return $self->{display} = $format;
930 return $display unless defined $format;
931 return $display = $format;
937 # Show nicely formatted complex number under its cartesian or polar form,
938 # depending on the current display format:
940 # . If a specific display format has been recorded for this object, use it.
941 # . Otherwise, use the generic current default for all complex numbers,
942 # which is a package global variable.
949 $format = $z->{display} if defined $z->{display};
951 return $z->stringify_polar if $format =~ /^p/i;
952 return $z->stringify_cartesian;
956 # ->stringify_cartesian
958 # Stringify as a cartesian representation 'a+bi'.
960 sub stringify_cartesian {
962 my ($x, $y) = @{$z->cartesian};
965 $x = int($x + ($x < 0 ? -1 : 1) * 1e-14)
966 if int(abs($x)) != int(abs($x) + 1e-14);
967 $y = int($y + ($y < 0 ? -1 : 1) * 1e-14)
968 if int(abs($y)) != int(abs($y) + 1e-14);
970 $re = "$x" if abs($x) >= 1e-14;
971 if ($y == 1) { $im = 'i' }
972 elsif ($y == -1) { $im = '-i' }
973 elsif (abs($y) >= 1e-14) { $im = $y . "i" }
976 $str = $re if defined $re;
977 $str .= "+$im" if defined $im;
980 $str = '0' unless $str;
988 # Stringify as a polar representation '[r,t]'.
990 sub stringify_polar {
992 my ($r, $t) = @{$z->polar};
996 return '[0,0]' if $r <= $eps;
1000 $nt = ($nt - int($nt)) * $tpi;
1001 $nt += $tpi if $nt < 0; # Range [0, 2pi]
1003 if (abs($nt) <= $eps) { $theta = 0 }
1004 elsif (abs(pi-$nt) <= $eps) { $theta = 'pi' }
1006 if (defined $theta) {
1007 $r = int($r + ($r < 0 ? -1 : 1) * $eps)
1008 if int(abs($r)) != int(abs($r) + $eps);
1009 $theta = int($theta + ($theta < 0 ? -1 : 1) * $eps)
1010 if ($theta ne 'pi' and
1011 int(abs($theta)) != int(abs($theta) + $eps));
1012 return "\[$r,$theta\]";
1016 # Okay, number is not a real. Try to identify pi/n and friends...
1019 $nt -= $tpi if $nt > pi;
1022 for ($k = 1, $kpi = pi; $k < 10; $k++, $kpi += pi) {
1023 $n = int($kpi / $nt + ($nt > 0 ? 1 : -1) * 0.5);
1024 if (abs($kpi/$n - $nt) <= $eps) {
1025 $theta = ($nt < 0 ? '-':'').
1026 ($k == 1 ? 'pi':"${k}pi").'/'.abs($n);
1031 $theta = $nt unless defined $theta;
1033 $r = int($r + ($r < 0 ? -1 : 1) * $eps)
1034 if int(abs($r)) != int(abs($r) + $eps);
1035 $theta = int($theta + ($theta < 0 ? -1 : 1) * $eps)
1036 if ($theta !~ m(^-?\d*pi/\d+$) and
1037 int(abs($theta)) != int(abs($theta) + $eps));
1039 return "\[$r,$theta\]";
1047 Math::Complex - complex numbers and associated mathematical functions
1053 $z = Math::Complex->make(5, 6);
1055 $j = cplxe(1, 2*pi/3);
1059 This package lets you create and manipulate complex numbers. By default,
1060 I<Perl> limits itself to real numbers, but an extra C<use> statement brings
1061 full complex support, along with a full set of mathematical functions
1062 typically associated with and/or extended to complex numbers.
1064 If you wonder what complex numbers are, they were invented to be able to solve
1065 the following equation:
1069 and by definition, the solution is noted I<i> (engineers use I<j> instead since
1070 I<i> usually denotes an intensity, but the name does not matter). The number
1071 I<i> is a pure I<imaginary> number.
1073 The arithmetics with pure imaginary numbers works just like you would expect
1074 it with real numbers... you just have to remember that
1080 5i + 7i = i * (5 + 7) = 12i
1081 4i - 3i = i * (4 - 3) = i
1086 Complex numbers are numbers that have both a real part and an imaginary
1087 part, and are usually noted:
1091 where C<a> is the I<real> part and C<b> is the I<imaginary> part. The
1092 arithmetic with complex numbers is straightforward. You have to
1093 keep track of the real and the imaginary parts, but otherwise the
1094 rules used for real numbers just apply:
1096 (4 + 3i) + (5 - 2i) = (4 + 5) + i(3 - 2) = 9 + i
1097 (2 + i) * (4 - i) = 2*4 + 4i -2i -i*i = 8 + 2i + 1 = 9 + 2i
1099 A graphical representation of complex numbers is possible in a plane
1100 (also called the I<complex plane>, but it's really a 2D plane).
1105 is the point whose coordinates are (a, b). Actually, it would
1106 be the vector originating from (0, 0) to (a, b). It follows that the addition
1107 of two complex numbers is a vectorial addition.
1109 Since there is a bijection between a point in the 2D plane and a complex
1110 number (i.e. the mapping is unique and reciprocal), a complex number
1111 can also be uniquely identified with polar coordinates:
1115 where C<rho> is the distance to the origin, and C<theta> the angle between
1116 the vector and the I<x> axis. There is a notation for this using the
1117 exponential form, which is:
1119 rho * exp(i * theta)
1121 where I<i> is the famous imaginary number introduced above. Conversion
1122 between this form and the cartesian form C<a + bi> is immediate:
1124 a = rho * cos(theta)
1125 b = rho * sin(theta)
1127 which is also expressed by this formula:
1129 z = rho * exp(i * theta) = rho * (cos theta + i * sin theta)
1131 In other words, it's the projection of the vector onto the I<x> and I<y>
1132 axes. Mathematicians call I<rho> the I<norm> or I<modulus> and I<theta>
1133 the I<argument> of the complex number. The I<norm> of C<z> will be
1136 The polar notation (also known as the trigonometric
1137 representation) is much more handy for performing multiplications and
1138 divisions of complex numbers, whilst the cartesian notation is better
1139 suited for additions and substractions. Real numbers are on the I<x>
1140 axis, and therefore I<theta> is zero.
1142 All the common operations that can be performed on a real number have
1143 been defined to work on complex numbers as well, and are merely
1144 I<extensions> of the operations defined on real numbers. This means
1145 they keep their natural meaning when there is no imaginary part, provided
1146 the number is within their definition set.
1148 For instance, the C<sqrt> routine which computes the square root of
1149 its argument is only defined for positive real numbers and yields a
1150 positive real number (it is an application from B<R+> to B<R+>).
1151 If we allow it to return a complex number, then it can be extended to
1152 negative real numbers to become an application from B<R> to B<C> (the
1153 set of complex numbers):
1155 sqrt(x) = x >= 0 ? sqrt(x) : sqrt(-x)*i
1157 It can also be extended to be an application from B<C> to B<C>,
1158 whilst its restriction to B<R> behaves as defined above by using
1159 the following definition:
1161 sqrt(z = [r,t]) = sqrt(r) * exp(i * t/2)
1163 Indeed, a negative real number can be noted C<[x,pi]>
1164 (the modulus I<x> is always positive, so C<[x,pi]> is really C<-x>, a
1166 and the above definition states that
1168 sqrt([x,pi]) = sqrt(x) * exp(i*pi/2) = [sqrt(x),pi/2] = sqrt(x)*i
1170 which is exactly what we had defined for negative real numbers above.
1172 All the common mathematical functions defined on real numbers that
1173 are extended to complex numbers share that same property of working
1174 I<as usual> when the imaginary part is zero (otherwise, it would not
1175 be called an extension, would it?).
1177 A I<new> operation possible on a complex number that is
1178 the identity for real numbers is called the I<conjugate>, and is noted
1179 with an horizontal bar above the number, or C<~z> here.
1186 z * ~z = (a + bi) * (a - bi) = a*a + b*b
1188 We saw that the norm of C<z> was noted C<abs(z)> and was defined as the
1189 distance to the origin, also known as:
1191 rho = abs(z) = sqrt(a*a + b*b)
1195 z * ~z = abs(z) ** 2
1197 If z is a pure real number (i.e. C<b == 0>), then the above yields:
1201 which is true (C<abs> has the regular meaning for real number, i.e. stands
1202 for the absolute value). This example explains why the norm of C<z> is
1203 noted C<abs(z)>: it extends the C<abs> function to complex numbers, yet
1204 is the regular C<abs> we know when the complex number actually has no
1205 imaginary part... This justifies I<a posteriori> our use of the C<abs>
1206 notation for the norm.
1210 Given the following notations:
1212 z1 = a + bi = r1 * exp(i * t1)
1213 z2 = c + di = r2 * exp(i * t2)
1214 z = <any complex or real number>
1216 the following (overloaded) operations are supported on complex numbers:
1218 z1 + z2 = (a + c) + i(b + d)
1219 z1 - z2 = (a - c) + i(b - d)
1220 z1 * z2 = (r1 * r2) * exp(i * (t1 + t2))
1221 z1 / z2 = (r1 / r2) * exp(i * (t1 - t2))
1222 z1 ** z2 = exp(z2 * log z1)
1224 abs(z1) = r1 = sqrt(a*a + b*b)
1225 sqrt(z1) = sqrt(r1) * exp(i * t1/2)
1226 exp(z1) = exp(a) * exp(i * b)
1227 log(z1) = log(r1) + i*t1
1228 sin(z1) = 1/2i (exp(i * z1) - exp(-i * z1))
1229 cos(z1) = 1/2 (exp(i * z1) + exp(-i * z1))
1231 atan2(z1, z2) = atan(z1/z2)
1233 The following extra operations are supported on both real and complex
1240 cbrt(z) = z ** (1/3)
1241 log10(z) = log(z) / log(10)
1242 logn(z, n) = log(z) / log(n)
1244 tan(z) = sin(z) / cos(z)
1250 asin(z) = -i * log(i*z + sqrt(1-z*z))
1251 acos(z) = -i * log(z + sqrt(z*z-1))
1252 atan(z) = i/2 * log((i+z) / (i-z))
1254 acsc(z) = asin(1 / z)
1255 asec(z) = acos(1 / z)
1256 acot(z) = -i/2 * log((i+z) / (z-i))
1258 sinh(z) = 1/2 (exp(z) - exp(-z))
1259 cosh(z) = 1/2 (exp(z) + exp(-z))
1260 tanh(z) = sinh(z) / cosh(z) = (exp(z) - exp(-z)) / (exp(z) + exp(-z))
1262 csch(z) = 1 / sinh(z)
1263 sech(z) = 1 / cosh(z)
1264 coth(z) = 1 / tanh(z)
1266 asinh(z) = log(z + sqrt(z*z+1))
1267 acosh(z) = log(z + sqrt(z*z-1))
1268 atanh(z) = 1/2 * log((1+z) / (1-z))
1270 acsch(z) = asinh(1 / z)
1271 asech(z) = acosh(1 / z)
1272 acoth(z) = atanh(1 / z) = 1/2 * log((1+z) / (z-1))
1274 I<log>, I<csc>, I<cot>, I<acsc>, I<acot>, I<csch>, I<coth>,
1275 I<acosech>, I<acotanh>, have aliases I<ln>, I<cosec>, I<cotan>,
1276 I<acosec>, I<acotan>, I<cosech>, I<cotanh>, I<acosech>, I<acotanh>,
1279 The I<root> function is available to compute all the I<n>
1280 roots of some complex, where I<n> is a strictly positive integer.
1281 There are exactly I<n> such roots, returned as a list. Getting the
1282 number mathematicians call C<j> such that:
1286 is a simple matter of writing:
1288 $j = ((root(1, 3))[1];
1290 The I<k>th root for C<z = [r,t]> is given by:
1292 (root(z, n))[k] = r**(1/n) * exp(i * (t + 2*k*pi)/n)
1294 The I<spaceship> comparison operator is also defined. In order to
1295 ensure its restriction to real numbers is conform to what you would
1296 expect, the comparison is run on the real part of the complex number
1297 first, and imaginary parts are compared only when the real parts
1302 To create a complex number, use either:
1304 $z = Math::Complex->make(3, 4);
1307 if you know the cartesian form of the number, or
1311 if you like. To create a number using the trigonometric form, use either:
1313 $z = Math::Complex->emake(5, pi/3);
1314 $x = cplxe(5, pi/3);
1316 instead. The first argument is the modulus, the second is the angle
1317 (in radians, the full circle is 2*pi). (Mnmemonic: C<e> is used as a
1318 notation for complex numbers in the trigonometric form).
1320 It is possible to write:
1322 $x = cplxe(-3, pi/4);
1324 but that will be silently converted into C<[3,-3pi/4]>, since the modulus
1325 must be positive (it represents the distance to the origin in the complex
1328 =head1 STRINGIFICATION
1330 When printed, a complex number is usually shown under its cartesian
1331 form I<a+bi>, but there are legitimate cases where the polar format
1332 I<[r,t]> is more appropriate.
1334 By calling the routine C<Math::Complex::display_format> and supplying either
1335 C<"polar"> or C<"cartesian">, you override the default display format,
1336 which is C<"cartesian">. Not supplying any argument returns the current
1339 This default can be overridden on a per-number basis by calling the
1340 C<display_format> method instead. As before, not supplying any argument
1341 returns the current display format for this number. Otherwise whatever you
1342 specify will be the new display format for I<this> particular number.
1348 Math::Complex::display_format('polar');
1349 $j = ((root(1, 3))[1];
1350 print "j = $j\n"; # Prints "j = [1,2pi/3]
1351 $j->display_format('cartesian');
1352 print "j = $j\n"; # Prints "j = -0.5+0.866025403784439i"
1354 The polar format attempts to emphasize arguments like I<k*pi/n>
1355 (where I<n> is a positive integer and I<k> an integer within [-9,+9]).
1359 Thanks to overloading, the handling of arithmetics with complex numbers
1360 is simple and almost transparent.
1362 Here are some examples:
1366 $j = cplxe(1, 2*pi/3); # $j ** 3 == 1
1367 print "j = $j, j**3 = ", $j ** 3, "\n";
1368 print "1 + j + j**2 = ", 1 + $j + $j**2, "\n";
1370 $z = -16 + 0*i; # Force it to be a complex
1371 print "sqrt($z) = ", sqrt($z), "\n";
1373 $k = exp(i * 2*pi/3);
1374 print "$j - $k = ", $j - $k, "\n";
1376 =head1 ERRORS DUE TO DIVISION BY ZERO
1378 The division (/) and the following functions
1397 cannot be computed for all arguments because that would mean dividing
1398 by zero. These situations cause fatal runtime errors looking like this
1400 cot(0): Division by zero.
1401 (Because in the definition of cot(0), the divisor sin(0) is 0)
1404 For the C<csc>, C<cot>, C<asec>, C<acsc>, C<csch>, C<coth>, C<asech>,
1405 C<acsch>, the argument cannot be C<0> (zero). For the C<atanh>,
1406 C<acoth>, the argument cannot be C<1> (one). For the C<atan>, C<acot>,
1407 the argument cannot be C<i> (the imaginary unit). For the C<tan>,
1408 C<sec>, C<tanh>, C<sech>, the argument cannot be I<pi/2 + k * pi>, where
1409 I<k> is any integer.
1413 Saying C<use Math::Complex;> exports many mathematical routines in the
1414 caller environment and even overrides some (C<sin>, C<cos>, C<sqrt>,
1415 C<log>, C<exp>). This is construed as a feature by the Authors,
1418 The code is not optimized for speed, although we try to use the cartesian
1419 form for addition-like operators and the trigonometric form for all
1420 multiplication-like operators.
1422 The arg() routine does not ensure the angle is within the range [-pi,+pi]
1423 (a side effect caused by multiplication and division using the trigonometric
1426 All routines expect to be given real or complex numbers. Don't attempt to
1427 use BigFloat, since Perl has currently no rule to disambiguate a '+'
1428 operation (for instance) between two overloaded entities.
1432 Raphael Manfredi <F<Raphael_Manfredi@grenoble.hp.com>>
1433 Jarkko Hietaniemi <F<jhi@iki.fi>>