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 $z2 = cplx($z2) unless ref $z2;
211 my ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
212 unless (defined $regular) {
213 $z1->set_cartesian([$re1 + $re2, $im1 + $im2]);
216 return (ref $z1)->make($re1 + $re2, $im1 + $im2);
225 my ($z1, $z2, $inverted) = @_;
226 my ($re1, $im1) = @{$z1->cartesian};
227 $z2 = cplx($z2) unless ref $z2;
228 my ($re2, $im2) = @{$z2->cartesian};
229 unless (defined $inverted) {
230 $z1->set_cartesian([$re1 - $re2, $im1 - $im2]);
234 (ref $z1)->make($re2 - $re1, $im2 - $im1) :
235 (ref $z1)->make($re1 - $re2, $im1 - $im2);
245 my ($z1, $z2, $regular) = @_;
246 my ($r1, $t1) = @{$z1->polar};
247 $z2 = cplxe(abs($z2), $z2 >= 0 ? 0 : pi) unless ref $z2;
248 my ($r2, $t2) = @{$z2->polar};
249 unless (defined $regular) {
250 $z1->set_polar([$r1 * $r2, $t1 + $t2]);
253 return (ref $z1)->emake($r1 * $r2, $t1 + $t2);
259 # Die on division by zero.
262 my $mess = "$_[0]: Division by zero.\n";
265 $mess .= "(Because in the definition of $_[0], the divisor ";
266 $mess .= "$_[1] " unless ($_[1] eq '0');
272 $mess .= "Died at $up[1] line $up[2].\n";
283 my ($z1, $z2, $inverted) = @_;
284 my ($r1, $t1) = @{$z1->polar};
285 $z2 = cplxe(abs($z2), $z2 >= 0 ? 0 : pi) unless ref $z2;
286 my ($r2, $t2) = @{$z2->polar};
287 unless (defined $inverted) {
288 _divbyzero "$z1/0" if ($r2 == 0);
289 $z1->set_polar([$r1 / $r2, $t1 - $t2]);
293 _divbyzero "$z2/0" if ($r1 == 0);
294 return (ref $z1)->emake($r2 / $r1, $t2 - $t1);
296 _divbyzero "$z1/0" if ($r2 == 0);
297 return (ref $z1)->emake($r1 / $r2, $t1 - $t2);
304 # Die on zero raised to the zeroth.
307 my $mess = "The zero raised to the zeroth power is not defined.\n";
311 $mess .= "Died at $up[1] line $up[2].\n";
319 # Computes z1**z2 = exp(z2 * log z1)).
322 my ($z1, $z2, $inverted) = @_;
325 _zerotozero if ($z1z and $z2z);
328 return 1 if ($z1z or $z2 == 1);
331 return 1 if ($z2z or $z1 == 1);
333 $z2 = cplx($z2) unless ref $z2;
334 unless (defined $inverted) {
335 my $z3 = exp($z2 * log $z1);
336 $z1->set_cartesian([@{$z3->cartesian}]);
339 return exp($z2 * log $z1) unless $inverted;
340 return exp($z1 * log $z2);
346 # Computes z1 <=> z2.
347 # Sorts on the real part first, then on the imaginary part. Thus 2-4i > 3+8i.
350 my ($z1, $z2, $inverted) = @_;
351 my ($re1, $im1) = ref $z1 ? @{$z1->cartesian} : ($z1, 0);
352 my ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
353 my $sgn = $inverted ? -1 : 1;
354 return $sgn * ($re1 <=> $re2) if $re1 != $re2;
355 return $sgn * ($im1 <=> $im2);
366 my ($r, $t) = @{$z->polar};
367 return (ref $z)->emake($r, pi + $t);
369 my ($re, $im) = @{$z->cartesian};
370 return (ref $z)->make(-$re, -$im);
376 # Compute complex's conjugate.
381 my ($r, $t) = @{$z->polar};
382 return (ref $z)->emake($r, -$t);
384 my ($re, $im) = @{$z->cartesian};
385 return (ref $z)->make($re, -$im);
391 # Compute complex's norm (rho).
395 return abs($z) unless ref $z;
396 my ($r, $t) = @{$z->polar};
403 # Compute complex's argument (theta).
407 return ($z < 0 ? pi : 0) unless ref $z;
408 my ($r, $t) = @{$z->polar};
419 $z = cplx($z, 0) unless ref $z;
420 my ($r, $t) = @{$z->polar};
421 return (ref $z)->emake(sqrt($r), $t/2);
427 # Compute cbrt(z) (cubic root).
431 return cplx($z, 0) ** (1/3) unless ref $z;
432 my ($r, $t) = @{$z->polar};
433 return (ref $z)->emake($r**(1/3), $t/3);
442 my $mess = "Root $_[0] not defined, root must be positive integer.\n";
446 $mess .= "Died at $up[1] line $up[2].\n";
454 # Computes all nth root for z, returning an array whose size is n.
455 # `n' must be a positive integer.
457 # The roots are given by (for k = 0..n-1):
459 # z^(1/n) = r^(1/n) (cos ((t+2 k pi)/n) + i sin ((t+2 k pi)/n))
463 _rootbad($n) if ($n < 1 or int($n) != $n);
464 my ($r, $t) = ref $z ? @{$z->polar} : (abs($z), $z >= 0 ? 0 : pi);
467 my $theta_inc = 2 * pi / $n;
468 my $rho = $r ** (1/$n);
470 my $complex = ref($z) || $package;
471 for ($k = 0, $theta = $t / $n; $k < $n; $k++, $theta += $theta_inc) {
472 push(@root, $complex->emake($rho, $theta));
484 return $z unless ref $z;
485 my ($re, $im) = @{$z->cartesian};
496 return 0 unless ref $z;
497 my ($re, $im) = @{$z->cartesian};
508 $z = cplx($z, 0) unless ref $z;
509 my ($x, $y) = @{$z->cartesian};
510 return (ref $z)->emake(exp($x), $y);
520 $z = cplx($z, 0) unless ref $z;
521 my ($x, $y) = @{$z->cartesian};
522 my ($r, $t) = @{$z->polar};
523 $t -= 2 * pi if ($t > pi() and $x < 0);
524 $t += 2 * pi if ($t < -pi() and $x < 0);
525 return (ref $z)->make(log($r), $t);
533 sub ln { Math::Complex::log(@_) }
544 return log(cplx($z, 0)) * log10inv unless ref $z;
545 my ($r, $t) = @{$z->polar};
546 return (ref $z)->make(log($r) * log10inv, $t * log10inv);
552 # Compute logn(z,n) = log(z) / log(n)
556 $z = cplx($z, 0) unless ref $z;
557 my $logn = $logn{$n};
558 $logn = $logn{$n} = log($n) unless defined $logn; # Cache log(n)
559 return log($z) / $logn;
565 # Compute cos(z) = (exp(iz) + exp(-iz))/2.
569 $z = cplx($z, 0) unless ref $z;
570 my ($x, $y) = @{$z->cartesian};
573 return (ref $z)->make(cos($x) * ($ey + $ey_1)/2,
574 sin($x) * ($ey_1 - $ey)/2);
580 # Compute sin(z) = (exp(iz) - exp(-iz))/2.
584 $z = cplx($z, 0) unless ref $z;
585 my ($x, $y) = @{$z->cartesian};
588 return (ref $z)->make(sin($x) * ($ey + $ey_1)/2,
589 cos($x) * ($ey - $ey_1)/2);
595 # Compute tan(z) = sin(z) / cos(z).
600 _divbyzero "tan($z)", "cos($z)" if ($cz == 0);
601 return sin($z) / $cz;
607 # Computes the secant sec(z) = 1 / cos(z).
612 _divbyzero "sec($z)", "cos($z)" if ($cz == 0);
619 # Computes the cosecant csc(z) = 1 / sin(z).
624 _divbyzero "csc($z)", "sin($z)" if ($sz == 0);
633 sub cosec { Math::Complex::csc(@_) }
638 # Computes cot(z) = 1 / tan(z).
643 _divbyzero "cot($z)", "sin($z)" if ($sz == 0);
644 return cos($z) / $sz;
652 sub cotan { Math::Complex::cot(@_) }
657 # Computes the arc cosine acos(z) = -i log(z + sqrt(z*z-1)).
661 $z = cplx($z, 0) unless ref $z;
662 return ~i * log($z + (Re($z) * Im($z) > 0 ? 1 : -1) * sqrt($z*$z - 1));
668 # Computes the arc sine asin(z) = -i log(iz + sqrt(1-z*z)).
672 $z = cplx($z, 0) unless ref $z;
673 return ~i * log(i * $z + sqrt(1 - $z*$z));
679 # Computes the arc tangent atan(z) = i/2 log((i+z) / (i-z)).
683 $z = cplx($z, 0) unless ref $z;
684 _divbyzero "atan($z)", "i - $z" if ($z == i);
685 return i/2*log((i + $z) / (i - $z));
691 # Computes the arc secant asec(z) = acos(1 / z).
695 _divbyzero "asec($z)", $z if ($z == 0);
702 # Computes the arc cosecant sec(z) = asin(1 / z).
706 _divbyzero "acsc($z)", $z if ($z == 0);
715 sub acosec { Math::Complex::acsc(@_) }
720 # Computes the arc cotangent acot(z) = -i/2 log((i+z) / (z-i))
724 $z = cplx($z, 0) unless ref $z;
725 _divbyzero "acot($z)", "$z - i" if ($z == i);
726 return i/-2 * log((i + $z) / ($z - i));
734 sub acotan { Math::Complex::acot(@_) }
739 # Computes the hyperbolic cosine cosh(z) = (exp(z) + exp(-z))/2.
748 my ($x, $y) = @{$z->cartesian};
751 return cplx(0.5 * ($ex + $ex_1), 0) if $real;
752 return (ref $z)->make(cos($y) * ($ex + $ex_1)/2,
753 sin($y) * ($ex - $ex_1)/2);
759 # Computes the hyperbolic sine sinh(z) = (exp(z) - exp(-z))/2.
768 my ($x, $y) = @{$z->cartesian};
771 return cplx(0.5 * ($ex - $ex_1), 0) if $real;
772 return (ref $z)->make(cos($y) * ($ex - $ex_1)/2,
773 sin($y) * ($ex + $ex_1)/2);
779 # Computes the hyperbolic tangent tanh(z) = sinh(z) / cosh(z).
784 _divbyzero "tanh($z)", "cosh($z)" if ($cz == 0);
785 return sinh($z) / $cz;
791 # Computes the hyperbolic secant sech(z) = 1 / cosh(z).
796 _divbyzero "sech($z)", "cosh($z)" if ($cz == 0);
803 # Computes the hyperbolic cosecant csch(z) = 1 / sinh(z).
808 _divbyzero "csch($z)", "sinh($z)" if ($sz == 0);
817 sub cosech { Math::Complex::csch(@_) }
822 # Computes the hyperbolic cotangent coth(z) = cosh(z) / sinh(z).
827 _divbyzero "coth($z)", "sinh($z)" if ($sz == 0);
828 return cosh($z) / $sz;
836 sub cotanh { Math::Complex::coth(@_) }
841 # Computes the arc hyperbolic cosine acosh(z) = log(z + sqrt(z*z-1)).
845 $z = cplx($z, 0) unless ref $z;
846 return log($z + sqrt($z*$z - 1));
852 # Computes the arc hyperbolic sine asinh(z) = log(z + sqrt(z*z-1))
856 $z = cplx($z, 0) unless ref $z;
857 return log($z + sqrt($z*$z + 1));
863 # Computes the arc hyperbolic tangent atanh(z) = 1/2 log((1+z) / (1-z)).
867 _divbyzero 'atanh(1)', "1 - $z" if ($z == 1);
868 $z = cplx($z, 0) unless ref $z;
869 my $cz = (1 + $z) / (1 - $z);
876 # Computes the hyperbolic arc secant asech(z) = acosh(1 / z).
880 _divbyzero 'asech(0)', $z if ($z == 0);
881 return acosh(1 / $z);
887 # Computes the hyperbolic arc cosecant acsch(z) = asinh(1 / z).
891 _divbyzero 'acsch(0)', $z if ($z == 0);
892 return asinh(1 / $z);
900 sub acosech { Math::Complex::acsch(@_) }
905 # Computes the arc hyperbolic cotangent acoth(z) = 1/2 log((1+z) / (z-1)).
909 _divbyzero 'acoth(1)', "$z - 1" if ($z == 1);
910 $z = cplx($z, 0) unless ref $z;
911 my $cz = (1 + $z) / ($z - 1);
920 sub acotanh { Math::Complex::acoth(@_) }
925 # Compute atan(z1/z2).
928 my ($z1, $z2, $inverted) = @_;
929 my ($re1, $im1) = ref $z1 ? @{$z1->cartesian} : ($z1, 0);
930 my ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
932 if (defined $inverted && $inverted) { # atan(z2/z1)
933 return pi * ($re2 > 0 ? 1 : -1) if $re1 == 0 && $im1 == 0;
936 return pi * ($re1 > 0 ? 1 : -1) if $re2 == 0 && $im2 == 0;
946 # Set (fetch if no argument) display format for all complex numbers that
947 # don't happen to have overrriden it via ->display_format
949 # When called as a method, this actually sets the display format for
950 # the current object.
952 # Valid object formats are 'c' and 'p' for cartesian and polar. The first
953 # letter is used actually, so the type can be fully spelled out for clarity.
959 if (ref $self) { # Called as a method
961 } else { # Regular procedure call
967 return defined $self->{display} ? $self->{display} : $display
968 unless defined $format;
969 return $self->{display} = $format;
972 return $display unless defined $format;
973 return $display = $format;
979 # Show nicely formatted complex number under its cartesian or polar form,
980 # depending on the current display format:
982 # . If a specific display format has been recorded for this object, use it.
983 # . Otherwise, use the generic current default for all complex numbers,
984 # which is a package global variable.
991 $format = $z->{display} if defined $z->{display};
993 return $z->stringify_polar if $format =~ /^p/i;
994 return $z->stringify_cartesian;
998 # ->stringify_cartesian
1000 # Stringify as a cartesian representation 'a+bi'.
1002 sub stringify_cartesian {
1004 my ($x, $y) = @{$z->cartesian};
1007 $x = int($x + ($x < 0 ? -1 : 1) * 1e-14)
1008 if int(abs($x)) != int(abs($x) + 1e-14);
1009 $y = int($y + ($y < 0 ? -1 : 1) * 1e-14)
1010 if int(abs($y)) != int(abs($y) + 1e-14);
1012 $re = "$x" if abs($x) >= 1e-14;
1013 if ($y == 1) { $im = 'i' }
1014 elsif ($y == -1) { $im = '-i' }
1015 elsif (abs($y) >= 1e-14) { $im = $y . "i" }
1018 $str = $re if defined $re;
1019 $str .= "+$im" if defined $im;
1022 $str = '0' unless $str;
1030 # Stringify as a polar representation '[r,t]'.
1032 sub stringify_polar {
1034 my ($r, $t) = @{$z->polar};
1038 return '[0,0]' if $r <= $eps;
1042 $nt = ($nt - int($nt)) * $tpi;
1043 $nt += $tpi if $nt < 0; # Range [0, 2pi]
1045 if (abs($nt) <= $eps) { $theta = 0 }
1046 elsif (abs(pi-$nt) <= $eps) { $theta = 'pi' }
1048 if (defined $theta) {
1049 $r = int($r + ($r < 0 ? -1 : 1) * $eps)
1050 if int(abs($r)) != int(abs($r) + $eps);
1051 $theta = int($theta + ($theta < 0 ? -1 : 1) * $eps)
1052 if ($theta ne 'pi' and
1053 int(abs($theta)) != int(abs($theta) + $eps));
1054 return "\[$r,$theta\]";
1058 # Okay, number is not a real. Try to identify pi/n and friends...
1061 $nt -= $tpi if $nt > pi;
1064 for ($k = 1, $kpi = pi; $k < 10; $k++, $kpi += pi) {
1065 $n = int($kpi / $nt + ($nt > 0 ? 1 : -1) * 0.5);
1066 if (abs($kpi/$n - $nt) <= $eps) {
1067 $theta = ($nt < 0 ? '-':'').
1068 ($k == 1 ? 'pi':"${k}pi").'/'.abs($n);
1073 $theta = $nt unless defined $theta;
1075 $r = int($r + ($r < 0 ? -1 : 1) * $eps)
1076 if int(abs($r)) != int(abs($r) + $eps);
1077 $theta = int($theta + ($theta < 0 ? -1 : 1) * $eps)
1078 if ($theta !~ m(^-?\d*pi/\d+$) and
1079 int(abs($theta)) != int(abs($theta) + $eps));
1081 return "\[$r,$theta\]";
1089 Math::Complex - complex numbers and associated mathematical functions
1095 $z = Math::Complex->make(5, 6);
1097 $j = cplxe(1, 2*pi/3);
1101 This package lets you create and manipulate complex numbers. By default,
1102 I<Perl> limits itself to real numbers, but an extra C<use> statement brings
1103 full complex support, along with a full set of mathematical functions
1104 typically associated with and/or extended to complex numbers.
1106 If you wonder what complex numbers are, they were invented to be able to solve
1107 the following equation:
1111 and by definition, the solution is noted I<i> (engineers use I<j> instead since
1112 I<i> usually denotes an intensity, but the name does not matter). The number
1113 I<i> is a pure I<imaginary> number.
1115 The arithmetics with pure imaginary numbers works just like you would expect
1116 it with real numbers... you just have to remember that
1122 5i + 7i = i * (5 + 7) = 12i
1123 4i - 3i = i * (4 - 3) = i
1128 Complex numbers are numbers that have both a real part and an imaginary
1129 part, and are usually noted:
1133 where C<a> is the I<real> part and C<b> is the I<imaginary> part. The
1134 arithmetic with complex numbers is straightforward. You have to
1135 keep track of the real and the imaginary parts, but otherwise the
1136 rules used for real numbers just apply:
1138 (4 + 3i) + (5 - 2i) = (4 + 5) + i(3 - 2) = 9 + i
1139 (2 + i) * (4 - i) = 2*4 + 4i -2i -i*i = 8 + 2i + 1 = 9 + 2i
1141 A graphical representation of complex numbers is possible in a plane
1142 (also called the I<complex plane>, but it's really a 2D plane).
1147 is the point whose coordinates are (a, b). Actually, it would
1148 be the vector originating from (0, 0) to (a, b). It follows that the addition
1149 of two complex numbers is a vectorial addition.
1151 Since there is a bijection between a point in the 2D plane and a complex
1152 number (i.e. the mapping is unique and reciprocal), a complex number
1153 can also be uniquely identified with polar coordinates:
1157 where C<rho> is the distance to the origin, and C<theta> the angle between
1158 the vector and the I<x> axis. There is a notation for this using the
1159 exponential form, which is:
1161 rho * exp(i * theta)
1163 where I<i> is the famous imaginary number introduced above. Conversion
1164 between this form and the cartesian form C<a + bi> is immediate:
1166 a = rho * cos(theta)
1167 b = rho * sin(theta)
1169 which is also expressed by this formula:
1171 z = rho * exp(i * theta) = rho * (cos theta + i * sin theta)
1173 In other words, it's the projection of the vector onto the I<x> and I<y>
1174 axes. Mathematicians call I<rho> the I<norm> or I<modulus> and I<theta>
1175 the I<argument> of the complex number. The I<norm> of C<z> will be
1178 The polar notation (also known as the trigonometric
1179 representation) is much more handy for performing multiplications and
1180 divisions of complex numbers, whilst the cartesian notation is better
1181 suited for additions and substractions. Real numbers are on the I<x>
1182 axis, and therefore I<theta> is zero.
1184 All the common operations that can be performed on a real number have
1185 been defined to work on complex numbers as well, and are merely
1186 I<extensions> of the operations defined on real numbers. This means
1187 they keep their natural meaning when there is no imaginary part, provided
1188 the number is within their definition set.
1190 For instance, the C<sqrt> routine which computes the square root of
1191 its argument is only defined for positive real numbers and yields a
1192 positive real number (it is an application from B<R+> to B<R+>).
1193 If we allow it to return a complex number, then it can be extended to
1194 negative real numbers to become an application from B<R> to B<C> (the
1195 set of complex numbers):
1197 sqrt(x) = x >= 0 ? sqrt(x) : sqrt(-x)*i
1199 It can also be extended to be an application from B<C> to B<C>,
1200 whilst its restriction to B<R> behaves as defined above by using
1201 the following definition:
1203 sqrt(z = [r,t]) = sqrt(r) * exp(i * t/2)
1205 Indeed, a negative real number can be noted C<[x,pi]>
1206 (the modulus I<x> is always positive, so C<[x,pi]> is really C<-x>, a
1208 and the above definition states that
1210 sqrt([x,pi]) = sqrt(x) * exp(i*pi/2) = [sqrt(x),pi/2] = sqrt(x)*i
1212 which is exactly what we had defined for negative real numbers above.
1214 All the common mathematical functions defined on real numbers that
1215 are extended to complex numbers share that same property of working
1216 I<as usual> when the imaginary part is zero (otherwise, it would not
1217 be called an extension, would it?).
1219 A I<new> operation possible on a complex number that is
1220 the identity for real numbers is called the I<conjugate>, and is noted
1221 with an horizontal bar above the number, or C<~z> here.
1228 z * ~z = (a + bi) * (a - bi) = a*a + b*b
1230 We saw that the norm of C<z> was noted C<abs(z)> and was defined as the
1231 distance to the origin, also known as:
1233 rho = abs(z) = sqrt(a*a + b*b)
1237 z * ~z = abs(z) ** 2
1239 If z is a pure real number (i.e. C<b == 0>), then the above yields:
1243 which is true (C<abs> has the regular meaning for real number, i.e. stands
1244 for the absolute value). This example explains why the norm of C<z> is
1245 noted C<abs(z)>: it extends the C<abs> function to complex numbers, yet
1246 is the regular C<abs> we know when the complex number actually has no
1247 imaginary part... This justifies I<a posteriori> our use of the C<abs>
1248 notation for the norm.
1252 Given the following notations:
1254 z1 = a + bi = r1 * exp(i * t1)
1255 z2 = c + di = r2 * exp(i * t2)
1256 z = <any complex or real number>
1258 the following (overloaded) operations are supported on complex numbers:
1260 z1 + z2 = (a + c) + i(b + d)
1261 z1 - z2 = (a - c) + i(b - d)
1262 z1 * z2 = (r1 * r2) * exp(i * (t1 + t2))
1263 z1 / z2 = (r1 / r2) * exp(i * (t1 - t2))
1264 z1 ** z2 = exp(z2 * log z1)
1266 abs(z1) = r1 = sqrt(a*a + b*b)
1267 sqrt(z1) = sqrt(r1) * exp(i * t1/2)
1268 exp(z1) = exp(a) * exp(i * b)
1269 log(z1) = log(r1) + i*t1
1270 sin(z1) = 1/2i (exp(i * z1) - exp(-i * z1))
1271 cos(z1) = 1/2 (exp(i * z1) + exp(-i * z1))
1273 atan2(z1, z2) = atan(z1/z2)
1275 The following extra operations are supported on both real and complex
1282 cbrt(z) = z ** (1/3)
1283 log10(z) = log(z) / log(10)
1284 logn(z, n) = log(z) / log(n)
1286 tan(z) = sin(z) / cos(z)
1292 asin(z) = -i * log(i*z + sqrt(1-z*z))
1293 acos(z) = -i * log(z + sqrt(z*z-1))
1294 atan(z) = i/2 * log((i+z) / (i-z))
1296 acsc(z) = asin(1 / z)
1297 asec(z) = acos(1 / z)
1298 acot(z) = -i/2 * log((i+z) / (z-i))
1300 sinh(z) = 1/2 (exp(z) - exp(-z))
1301 cosh(z) = 1/2 (exp(z) + exp(-z))
1302 tanh(z) = sinh(z) / cosh(z) = (exp(z) - exp(-z)) / (exp(z) + exp(-z))
1304 csch(z) = 1 / sinh(z)
1305 sech(z) = 1 / cosh(z)
1306 coth(z) = 1 / tanh(z)
1308 asinh(z) = log(z + sqrt(z*z+1))
1309 acosh(z) = log(z + sqrt(z*z-1))
1310 atanh(z) = 1/2 * log((1+z) / (1-z))
1312 acsch(z) = asinh(1 / z)
1313 asech(z) = acosh(1 / z)
1314 acoth(z) = atanh(1 / z) = 1/2 * log((1+z) / (z-1))
1316 I<log>, I<csc>, I<cot>, I<acsc>, I<acot>, I<csch>, I<coth>,
1317 I<acosech>, I<acotanh>, have aliases I<ln>, I<cosec>, I<cotan>,
1318 I<acosec>, I<acotan>, I<cosech>, I<cotanh>, I<acosech>, I<acotanh>,
1321 The I<root> function is available to compute all the I<n>
1322 roots of some complex, where I<n> is a strictly positive integer.
1323 There are exactly I<n> such roots, returned as a list. Getting the
1324 number mathematicians call C<j> such that:
1328 is a simple matter of writing:
1330 $j = ((root(1, 3))[1];
1332 The I<k>th root for C<z = [r,t]> is given by:
1334 (root(z, n))[k] = r**(1/n) * exp(i * (t + 2*k*pi)/n)
1336 The I<spaceship> comparison operator, E<lt>=E<gt>, is also defined. In
1337 order to ensure its restriction to real numbers is conform to what you
1338 would expect, the comparison is run on the real part of the complex
1339 number first, and imaginary parts are compared only when the real
1344 To create a complex number, use either:
1346 $z = Math::Complex->make(3, 4);
1349 if you know the cartesian form of the number, or
1353 if you like. To create a number using the trigonometric form, use either:
1355 $z = Math::Complex->emake(5, pi/3);
1356 $x = cplxe(5, pi/3);
1358 instead. The first argument is the modulus, the second is the angle
1359 (in radians, the full circle is 2*pi). (Mnmemonic: C<e> is used as a
1360 notation for complex numbers in the trigonometric form).
1362 It is possible to write:
1364 $x = cplxe(-3, pi/4);
1366 but that will be silently converted into C<[3,-3pi/4]>, since the modulus
1367 must be positive (it represents the distance to the origin in the complex
1370 =head1 STRINGIFICATION
1372 When printed, a complex number is usually shown under its cartesian
1373 form I<a+bi>, but there are legitimate cases where the polar format
1374 I<[r,t]> is more appropriate.
1376 By calling the routine C<Math::Complex::display_format> and supplying either
1377 C<"polar"> or C<"cartesian">, you override the default display format,
1378 which is C<"cartesian">. Not supplying any argument returns the current
1381 This default can be overridden on a per-number basis by calling the
1382 C<display_format> method instead. As before, not supplying any argument
1383 returns the current display format for this number. Otherwise whatever you
1384 specify will be the new display format for I<this> particular number.
1390 Math::Complex::display_format('polar');
1391 $j = ((root(1, 3))[1];
1392 print "j = $j\n"; # Prints "j = [1,2pi/3]
1393 $j->display_format('cartesian');
1394 print "j = $j\n"; # Prints "j = -0.5+0.866025403784439i"
1396 The polar format attempts to emphasize arguments like I<k*pi/n>
1397 (where I<n> is a positive integer and I<k> an integer within [-9,+9]).
1401 Thanks to overloading, the handling of arithmetics with complex numbers
1402 is simple and almost transparent.
1404 Here are some examples:
1408 $j = cplxe(1, 2*pi/3); # $j ** 3 == 1
1409 print "j = $j, j**3 = ", $j ** 3, "\n";
1410 print "1 + j + j**2 = ", 1 + $j + $j**2, "\n";
1412 $z = -16 + 0*i; # Force it to be a complex
1413 print "sqrt($z) = ", sqrt($z), "\n";
1415 $k = exp(i * 2*pi/3);
1416 print "$j - $k = ", $j - $k, "\n";
1418 =head1 ERRORS DUE TO DIVISION BY ZERO
1420 The division (/) and the following functions
1439 cannot be computed for all arguments because that would mean dividing
1440 by zero. These situations cause fatal runtime errors looking like this
1442 cot(0): Division by zero.
1443 (Because in the definition of cot(0), the divisor sin(0) is 0)
1446 For the C<csc>, C<cot>, C<asec>, C<acsc>, C<csch>, C<coth>, C<asech>,
1447 C<acsch>, the argument cannot be C<0> (zero). For the C<atanh>,
1448 C<acoth>, the argument cannot be C<1> (one). For the C<atan>, C<acot>,
1449 the argument cannot be C<i> (the imaginary unit). For the C<tan>,
1450 C<sec>, C<tanh>, C<sech>, the argument cannot be I<pi/2 + k * pi>, where
1451 I<k> is any integer.
1455 Saying C<use Math::Complex;> exports many mathematical routines in the
1456 caller environment and even overrides some (C<sin>, C<cos>, C<sqrt>,
1457 C<log>, C<exp>). This is construed as a feature by the Authors,
1460 The code is not optimized for speed, although we try to use the cartesian
1461 form for addition-like operators and the trigonometric form for all
1462 multiplication-like operators.
1464 The arg() routine does not ensure the angle is within the range [-pi,+pi]
1465 (a side effect caused by multiplication and division using the trigonometric
1468 All routines expect to be given real or complex numbers. Don't attempt to
1469 use BigFloat, since Perl has currently no rule to disambiguate a '+'
1470 operation (for instance) between two overloaded entities.
1474 Raphael Manfredi <F<Raphael_Manfredi@grenoble.hp.com>> and
1475 Jarkko Hietaniemi <F<jhi@iki.fi>>.