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) = @_;
323 _zerotozero if ($z1 == 0 and $z2 == 0);
324 return 1 if ($z1 == 1);
325 return 0 if ($z2 == 0);
326 $z2 = cplx($z2) unless ref $z2;
327 unless (defined $inverted) {
328 my $z3 = exp($z2 * log $z1);
329 $z1->set_cartesian([@{$z3->cartesian}]);
332 return exp($z2 * log $z1) unless $inverted;
333 return exp($z1 * log $z2);
339 # Computes z1 <=> z2.
340 # Sorts on the real part first, then on the imaginary part. Thus 2-4i > 3+8i.
343 my ($z1, $z2, $inverted) = @_;
344 my ($re1, $im1) = ref $z1 ? @{$z1->cartesian} : ($z1, 0);
345 my ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
346 my $sgn = $inverted ? -1 : 1;
347 return $sgn * ($re1 <=> $re2) if $re1 != $re2;
348 return $sgn * ($im1 <=> $im2);
359 my ($r, $t) = @{$z->polar};
360 return (ref $z)->emake($r, pi + $t);
362 my ($re, $im) = @{$z->cartesian};
363 return (ref $z)->make(-$re, -$im);
369 # Compute complex's conjugate.
374 my ($r, $t) = @{$z->polar};
375 return (ref $z)->emake($r, -$t);
377 my ($re, $im) = @{$z->cartesian};
378 return (ref $z)->make($re, -$im);
384 # Compute complex's norm (rho).
388 return abs($z) unless ref $z;
389 my ($r, $t) = @{$z->polar};
396 # Compute complex's argument (theta).
400 return ($z < 0 ? pi : 0) unless ref $z;
401 my ($r, $t) = @{$z->polar};
412 $z = cplx($z, 0) unless ref $z;
413 my ($r, $t) = @{$z->polar};
414 return (ref $z)->emake(sqrt($r), $t/2);
420 # Compute cbrt(z) (cubic root).
424 return cplx($z, 0) ** (1/3) unless ref $z;
425 my ($r, $t) = @{$z->polar};
426 return (ref $z)->emake($r**(1/3), $t/3);
435 my $mess = "Root $_[0] not defined, root must be positive integer.\n";
439 $mess .= "Died at $up[1] line $up[2].\n";
447 # Computes all nth root for z, returning an array whose size is n.
448 # `n' must be a positive integer.
450 # The roots are given by (for k = 0..n-1):
452 # z^(1/n) = r^(1/n) (cos ((t+2 k pi)/n) + i sin ((t+2 k pi)/n))
456 _rootbad($n) if ($n < 1 or int($n) != $n);
457 my ($r, $t) = ref $z ? @{$z->polar} : (abs($z), $z >= 0 ? 0 : pi);
460 my $theta_inc = 2 * pi / $n;
461 my $rho = $r ** (1/$n);
463 my $complex = ref($z) || $package;
464 for ($k = 0, $theta = $t / $n; $k < $n; $k++, $theta += $theta_inc) {
465 push(@root, $complex->emake($rho, $theta));
477 return $z unless ref $z;
478 my ($re, $im) = @{$z->cartesian};
489 return 0 unless ref $z;
490 my ($re, $im) = @{$z->cartesian};
501 $z = cplx($z, 0) unless ref $z;
502 my ($x, $y) = @{$z->cartesian};
503 return (ref $z)->emake(exp($x), $y);
513 $z = cplx($z, 0) unless ref $z;
514 my ($x, $y) = @{$z->cartesian};
515 my ($r, $t) = @{$z->polar};
516 $t -= 2 * pi if ($t > pi() and $x < 0);
517 $t += 2 * pi if ($t < -pi() and $x < 0);
518 return (ref $z)->make(log($r), $t);
526 sub ln { Math::Complex::log(@_) }
537 return log(cplx($z, 0)) * log10inv unless ref $z;
538 my ($r, $t) = @{$z->polar};
539 return (ref $z)->make(log($r) * log10inv, $t * log10inv);
545 # Compute logn(z,n) = log(z) / log(n)
549 $z = cplx($z, 0) unless ref $z;
550 my $logn = $logn{$n};
551 $logn = $logn{$n} = log($n) unless defined $logn; # Cache log(n)
552 return log($z) / $logn;
558 # Compute cos(z) = (exp(iz) + exp(-iz))/2.
562 $z = cplx($z, 0) unless ref $z;
563 my ($x, $y) = @{$z->cartesian};
566 return (ref $z)->make(cos($x) * ($ey + $ey_1)/2,
567 sin($x) * ($ey_1 - $ey)/2);
573 # Compute sin(z) = (exp(iz) - exp(-iz))/2.
577 $z = cplx($z, 0) unless ref $z;
578 my ($x, $y) = @{$z->cartesian};
581 return (ref $z)->make(sin($x) * ($ey + $ey_1)/2,
582 cos($x) * ($ey - $ey_1)/2);
588 # Compute tan(z) = sin(z) / cos(z).
593 _divbyzero "tan($z)", "cos($z)" if ($cz == 0);
594 return sin($z) / $cz;
600 # Computes the secant sec(z) = 1 / cos(z).
605 _divbyzero "sec($z)", "cos($z)" if ($cz == 0);
612 # Computes the cosecant csc(z) = 1 / sin(z).
617 _divbyzero "csc($z)", "sin($z)" if ($sz == 0);
626 sub cosec { Math::Complex::csc(@_) }
631 # Computes cot(z) = 1 / tan(z).
636 _divbyzero "cot($z)", "sin($z)" if ($sz == 0);
637 return cos($z) / $sz;
645 sub cotan { Math::Complex::cot(@_) }
650 # Computes the arc cosine acos(z) = -i log(z + sqrt(z*z-1)).
654 $z = cplx($z, 0) unless ref $z;
655 return ~i * log($z + (Re($z) * Im($z) > 0 ? 1 : -1) * sqrt($z*$z - 1));
661 # Computes the arc sine asin(z) = -i log(iz + sqrt(1-z*z)).
665 $z = cplx($z, 0) unless ref $z;
666 return ~i * log(i * $z + sqrt(1 - $z*$z));
672 # Computes the arc tangent atan(z) = i/2 log((i+z) / (i-z)).
676 $z = cplx($z, 0) unless ref $z;
677 _divbyzero "atan($z)", "i - $z" if ($z == i);
678 return i/2*log((i + $z) / (i - $z));
684 # Computes the arc secant asec(z) = acos(1 / z).
688 _divbyzero "asec($z)", $z if ($z == 0);
695 # Computes the arc cosecant sec(z) = asin(1 / z).
699 _divbyzero "acsc($z)", $z if ($z == 0);
708 sub acosec { Math::Complex::acsc(@_) }
713 # Computes the arc cotangent acot(z) = -i/2 log((i+z) / (z-i))
717 $z = cplx($z, 0) unless ref $z;
718 _divbyzero "acot($z)", "$z - i" if ($z == i);
719 return i/-2 * log((i + $z) / ($z - i));
727 sub acotan { Math::Complex::acot(@_) }
732 # Computes the hyperbolic cosine cosh(z) = (exp(z) + exp(-z))/2.
741 my ($x, $y) = @{$z->cartesian};
744 return ($ex + $ex_1)/2 if $real;
745 return (ref $z)->make(cos($y) * ($ex + $ex_1)/2,
746 sin($y) * ($ex - $ex_1)/2);
752 # Computes the hyperbolic sine sinh(z) = (exp(z) - exp(-z))/2.
761 my ($x, $y) = @{$z->cartesian};
764 return ($ex - $ex_1)/2 if $real;
765 return (ref $z)->make(cos($y) * ($ex - $ex_1)/2,
766 sin($y) * ($ex + $ex_1)/2);
772 # Computes the hyperbolic tangent tanh(z) = sinh(z) / cosh(z).
777 _divbyzero "tanh($z)", "cosh($z)" if ($cz == 0);
778 return sinh($z) / $cz;
784 # Computes the hyperbolic secant sech(z) = 1 / cosh(z).
789 _divbyzero "sech($z)", "cosh($z)" if ($cz == 0);
796 # Computes the hyperbolic cosecant csch(z) = 1 / sinh(z).
801 _divbyzero "csch($z)", "sinh($z)" if ($sz == 0);
810 sub cosech { Math::Complex::csch(@_) }
815 # Computes the hyperbolic cotangent coth(z) = cosh(z) / sinh(z).
820 _divbyzero "coth($z)", "sinh($z)" if ($sz == 0);
821 return cosh($z) / $sz;
829 sub cotanh { Math::Complex::coth(@_) }
834 # Computes the arc hyperbolic cosine acosh(z) = log(z + sqrt(z*z-1)).
838 $z = cplx($z, 0) unless ref $z;
839 return log($z + sqrt($z*$z - 1));
845 # Computes the arc hyperbolic sine asinh(z) = log(z + sqrt(z*z-1))
849 $z = cplx($z, 0) unless ref $z;
850 return log($z + sqrt($z*$z + 1));
856 # Computes the arc hyperbolic tangent atanh(z) = 1/2 log((1+z) / (1-z)).
860 _divbyzero 'atanh(1)', "1 - $z" if ($z == 1);
861 $z = cplx($z, 0) unless ref $z;
862 my $cz = (1 + $z) / (1 - $z);
869 # Computes the hyperbolic arc secant asech(z) = acosh(1 / z).
873 _divbyzero 'asech(0)', $z if ($z == 0);
874 return acosh(1 / $z);
880 # Computes the hyperbolic arc cosecant acsch(z) = asinh(1 / z).
884 _divbyzero 'acsch(0)', $z if ($z == 0);
885 return asinh(1 / $z);
893 sub acosech { Math::Complex::acsch(@_) }
898 # Computes the arc hyperbolic cotangent acoth(z) = 1/2 log((1+z) / (z-1)).
902 _divbyzero 'acoth(1)', "$z - 1" if ($z == 1);
903 $z = cplx($z, 0) unless ref $z;
904 my $cz = (1 + $z) / ($z - 1);
913 sub acotanh { Math::Complex::acoth(@_) }
918 # Compute atan(z1/z2).
921 my ($z1, $z2, $inverted) = @_;
922 my ($re1, $im1) = ref $z1 ? @{$z1->cartesian} : ($z1, 0);
923 my ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
925 if (defined $inverted && $inverted) { # atan(z2/z1)
926 return pi * ($re2 > 0 ? 1 : -1) if $re1 == 0 && $im1 == 0;
929 return pi * ($re1 > 0 ? 1 : -1) if $re2 == 0 && $im2 == 0;
939 # Set (fetch if no argument) display format for all complex numbers that
940 # don't happen to have overrriden it via ->display_format
942 # When called as a method, this actually sets the display format for
943 # the current object.
945 # Valid object formats are 'c' and 'p' for cartesian and polar. The first
946 # letter is used actually, so the type can be fully spelled out for clarity.
952 if (ref $self) { # Called as a method
954 } else { # Regular procedure call
960 return defined $self->{display} ? $self->{display} : $display
961 unless defined $format;
962 return $self->{display} = $format;
965 return $display unless defined $format;
966 return $display = $format;
972 # Show nicely formatted complex number under its cartesian or polar form,
973 # depending on the current display format:
975 # . If a specific display format has been recorded for this object, use it.
976 # . Otherwise, use the generic current default for all complex numbers,
977 # which is a package global variable.
984 $format = $z->{display} if defined $z->{display};
986 return $z->stringify_polar if $format =~ /^p/i;
987 return $z->stringify_cartesian;
991 # ->stringify_cartesian
993 # Stringify as a cartesian representation 'a+bi'.
995 sub stringify_cartesian {
997 my ($x, $y) = @{$z->cartesian};
1000 $x = int($x + ($x < 0 ? -1 : 1) * 1e-14)
1001 if int(abs($x)) != int(abs($x) + 1e-14);
1002 $y = int($y + ($y < 0 ? -1 : 1) * 1e-14)
1003 if int(abs($y)) != int(abs($y) + 1e-14);
1005 $re = "$x" if abs($x) >= 1e-14;
1006 if ($y == 1) { $im = 'i' }
1007 elsif ($y == -1) { $im = '-i' }
1008 elsif (abs($y) >= 1e-14) { $im = $y . "i" }
1011 $str = $re if defined $re;
1012 $str .= "+$im" if defined $im;
1015 $str = '0' unless $str;
1023 # Stringify as a polar representation '[r,t]'.
1025 sub stringify_polar {
1027 my ($r, $t) = @{$z->polar};
1031 return '[0,0]' if $r <= $eps;
1035 $nt = ($nt - int($nt)) * $tpi;
1036 $nt += $tpi if $nt < 0; # Range [0, 2pi]
1038 if (abs($nt) <= $eps) { $theta = 0 }
1039 elsif (abs(pi-$nt) <= $eps) { $theta = 'pi' }
1041 if (defined $theta) {
1042 $r = int($r + ($r < 0 ? -1 : 1) * $eps)
1043 if int(abs($r)) != int(abs($r) + $eps);
1044 $theta = int($theta + ($theta < 0 ? -1 : 1) * $eps)
1045 if ($theta ne 'pi' and
1046 int(abs($theta)) != int(abs($theta) + $eps));
1047 return "\[$r,$theta\]";
1051 # Okay, number is not a real. Try to identify pi/n and friends...
1054 $nt -= $tpi if $nt > pi;
1057 for ($k = 1, $kpi = pi; $k < 10; $k++, $kpi += pi) {
1058 $n = int($kpi / $nt + ($nt > 0 ? 1 : -1) * 0.5);
1059 if (abs($kpi/$n - $nt) <= $eps) {
1060 $theta = ($nt < 0 ? '-':'').
1061 ($k == 1 ? 'pi':"${k}pi").'/'.abs($n);
1066 $theta = $nt unless defined $theta;
1068 $r = int($r + ($r < 0 ? -1 : 1) * $eps)
1069 if int(abs($r)) != int(abs($r) + $eps);
1070 $theta = int($theta + ($theta < 0 ? -1 : 1) * $eps)
1071 if ($theta !~ m(^-?\d*pi/\d+$) and
1072 int(abs($theta)) != int(abs($theta) + $eps));
1074 return "\[$r,$theta\]";
1082 Math::Complex - complex numbers and associated mathematical functions
1088 $z = Math::Complex->make(5, 6);
1090 $j = cplxe(1, 2*pi/3);
1094 This package lets you create and manipulate complex numbers. By default,
1095 I<Perl> limits itself to real numbers, but an extra C<use> statement brings
1096 full complex support, along with a full set of mathematical functions
1097 typically associated with and/or extended to complex numbers.
1099 If you wonder what complex numbers are, they were invented to be able to solve
1100 the following equation:
1104 and by definition, the solution is noted I<i> (engineers use I<j> instead since
1105 I<i> usually denotes an intensity, but the name does not matter). The number
1106 I<i> is a pure I<imaginary> number.
1108 The arithmetics with pure imaginary numbers works just like you would expect
1109 it with real numbers... you just have to remember that
1115 5i + 7i = i * (5 + 7) = 12i
1116 4i - 3i = i * (4 - 3) = i
1121 Complex numbers are numbers that have both a real part and an imaginary
1122 part, and are usually noted:
1126 where C<a> is the I<real> part and C<b> is the I<imaginary> part. The
1127 arithmetic with complex numbers is straightforward. You have to
1128 keep track of the real and the imaginary parts, but otherwise the
1129 rules used for real numbers just apply:
1131 (4 + 3i) + (5 - 2i) = (4 + 5) + i(3 - 2) = 9 + i
1132 (2 + i) * (4 - i) = 2*4 + 4i -2i -i*i = 8 + 2i + 1 = 9 + 2i
1134 A graphical representation of complex numbers is possible in a plane
1135 (also called the I<complex plane>, but it's really a 2D plane).
1140 is the point whose coordinates are (a, b). Actually, it would
1141 be the vector originating from (0, 0) to (a, b). It follows that the addition
1142 of two complex numbers is a vectorial addition.
1144 Since there is a bijection between a point in the 2D plane and a complex
1145 number (i.e. the mapping is unique and reciprocal), a complex number
1146 can also be uniquely identified with polar coordinates:
1150 where C<rho> is the distance to the origin, and C<theta> the angle between
1151 the vector and the I<x> axis. There is a notation for this using the
1152 exponential form, which is:
1154 rho * exp(i * theta)
1156 where I<i> is the famous imaginary number introduced above. Conversion
1157 between this form and the cartesian form C<a + bi> is immediate:
1159 a = rho * cos(theta)
1160 b = rho * sin(theta)
1162 which is also expressed by this formula:
1164 z = rho * exp(i * theta) = rho * (cos theta + i * sin theta)
1166 In other words, it's the projection of the vector onto the I<x> and I<y>
1167 axes. Mathematicians call I<rho> the I<norm> or I<modulus> and I<theta>
1168 the I<argument> of the complex number. The I<norm> of C<z> will be
1171 The polar notation (also known as the trigonometric
1172 representation) is much more handy for performing multiplications and
1173 divisions of complex numbers, whilst the cartesian notation is better
1174 suited for additions and substractions. Real numbers are on the I<x>
1175 axis, and therefore I<theta> is zero.
1177 All the common operations that can be performed on a real number have
1178 been defined to work on complex numbers as well, and are merely
1179 I<extensions> of the operations defined on real numbers. This means
1180 they keep their natural meaning when there is no imaginary part, provided
1181 the number is within their definition set.
1183 For instance, the C<sqrt> routine which computes the square root of
1184 its argument is only defined for positive real numbers and yields a
1185 positive real number (it is an application from B<R+> to B<R+>).
1186 If we allow it to return a complex number, then it can be extended to
1187 negative real numbers to become an application from B<R> to B<C> (the
1188 set of complex numbers):
1190 sqrt(x) = x >= 0 ? sqrt(x) : sqrt(-x)*i
1192 It can also be extended to be an application from B<C> to B<C>,
1193 whilst its restriction to B<R> behaves as defined above by using
1194 the following definition:
1196 sqrt(z = [r,t]) = sqrt(r) * exp(i * t/2)
1198 Indeed, a negative real number can be noted C<[x,pi]>
1199 (the modulus I<x> is always positive, so C<[x,pi]> is really C<-x>, a
1201 and the above definition states that
1203 sqrt([x,pi]) = sqrt(x) * exp(i*pi/2) = [sqrt(x),pi/2] = sqrt(x)*i
1205 which is exactly what we had defined for negative real numbers above.
1207 All the common mathematical functions defined on real numbers that
1208 are extended to complex numbers share that same property of working
1209 I<as usual> when the imaginary part is zero (otherwise, it would not
1210 be called an extension, would it?).
1212 A I<new> operation possible on a complex number that is
1213 the identity for real numbers is called the I<conjugate>, and is noted
1214 with an horizontal bar above the number, or C<~z> here.
1221 z * ~z = (a + bi) * (a - bi) = a*a + b*b
1223 We saw that the norm of C<z> was noted C<abs(z)> and was defined as the
1224 distance to the origin, also known as:
1226 rho = abs(z) = sqrt(a*a + b*b)
1230 z * ~z = abs(z) ** 2
1232 If z is a pure real number (i.e. C<b == 0>), then the above yields:
1236 which is true (C<abs> has the regular meaning for real number, i.e. stands
1237 for the absolute value). This example explains why the norm of C<z> is
1238 noted C<abs(z)>: it extends the C<abs> function to complex numbers, yet
1239 is the regular C<abs> we know when the complex number actually has no
1240 imaginary part... This justifies I<a posteriori> our use of the C<abs>
1241 notation for the norm.
1245 Given the following notations:
1247 z1 = a + bi = r1 * exp(i * t1)
1248 z2 = c + di = r2 * exp(i * t2)
1249 z = <any complex or real number>
1251 the following (overloaded) operations are supported on complex numbers:
1253 z1 + z2 = (a + c) + i(b + d)
1254 z1 - z2 = (a - c) + i(b - d)
1255 z1 * z2 = (r1 * r2) * exp(i * (t1 + t2))
1256 z1 / z2 = (r1 / r2) * exp(i * (t1 - t2))
1257 z1 ** z2 = exp(z2 * log z1)
1259 abs(z1) = r1 = sqrt(a*a + b*b)
1260 sqrt(z1) = sqrt(r1) * exp(i * t1/2)
1261 exp(z1) = exp(a) * exp(i * b)
1262 log(z1) = log(r1) + i*t1
1263 sin(z1) = 1/2i (exp(i * z1) - exp(-i * z1))
1264 cos(z1) = 1/2 (exp(i * z1) + exp(-i * z1))
1266 atan2(z1, z2) = atan(z1/z2)
1268 The following extra operations are supported on both real and complex
1275 cbrt(z) = z ** (1/3)
1276 log10(z) = log(z) / log(10)
1277 logn(z, n) = log(z) / log(n)
1279 tan(z) = sin(z) / cos(z)
1285 asin(z) = -i * log(i*z + sqrt(1-z*z))
1286 acos(z) = -i * log(z + sqrt(z*z-1))
1287 atan(z) = i/2 * log((i+z) / (i-z))
1289 acsc(z) = asin(1 / z)
1290 asec(z) = acos(1 / z)
1291 acot(z) = -i/2 * log((i+z) / (z-i))
1293 sinh(z) = 1/2 (exp(z) - exp(-z))
1294 cosh(z) = 1/2 (exp(z) + exp(-z))
1295 tanh(z) = sinh(z) / cosh(z) = (exp(z) - exp(-z)) / (exp(z) + exp(-z))
1297 csch(z) = 1 / sinh(z)
1298 sech(z) = 1 / cosh(z)
1299 coth(z) = 1 / tanh(z)
1301 asinh(z) = log(z + sqrt(z*z+1))
1302 acosh(z) = log(z + sqrt(z*z-1))
1303 atanh(z) = 1/2 * log((1+z) / (1-z))
1305 acsch(z) = asinh(1 / z)
1306 asech(z) = acosh(1 / z)
1307 acoth(z) = atanh(1 / z) = 1/2 * log((1+z) / (z-1))
1309 I<log>, I<csc>, I<cot>, I<acsc>, I<acot>, I<csch>, I<coth>,
1310 I<acosech>, I<acotanh>, have aliases I<ln>, I<cosec>, I<cotan>,
1311 I<acosec>, I<acotan>, I<cosech>, I<cotanh>, I<acosech>, I<acotanh>,
1314 The I<root> function is available to compute all the I<n>
1315 roots of some complex, where I<n> is a strictly positive integer.
1316 There are exactly I<n> such roots, returned as a list. Getting the
1317 number mathematicians call C<j> such that:
1321 is a simple matter of writing:
1323 $j = ((root(1, 3))[1];
1325 The I<k>th root for C<z = [r,t]> is given by:
1327 (root(z, n))[k] = r**(1/n) * exp(i * (t + 2*k*pi)/n)
1329 The I<spaceship> comparison operator is also defined. In order to
1330 ensure its restriction to real numbers is conform to what you would
1331 expect, the comparison is run on the real part of the complex number
1332 first, and imaginary parts are compared only when the real parts
1337 To create a complex number, use either:
1339 $z = Math::Complex->make(3, 4);
1342 if you know the cartesian form of the number, or
1346 if you like. To create a number using the trigonometric form, use either:
1348 $z = Math::Complex->emake(5, pi/3);
1349 $x = cplxe(5, pi/3);
1351 instead. The first argument is the modulus, the second is the angle
1352 (in radians, the full circle is 2*pi). (Mnmemonic: C<e> is used as a
1353 notation for complex numbers in the trigonometric form).
1355 It is possible to write:
1357 $x = cplxe(-3, pi/4);
1359 but that will be silently converted into C<[3,-3pi/4]>, since the modulus
1360 must be positive (it represents the distance to the origin in the complex
1363 =head1 STRINGIFICATION
1365 When printed, a complex number is usually shown under its cartesian
1366 form I<a+bi>, but there are legitimate cases where the polar format
1367 I<[r,t]> is more appropriate.
1369 By calling the routine C<Math::Complex::display_format> and supplying either
1370 C<"polar"> or C<"cartesian">, you override the default display format,
1371 which is C<"cartesian">. Not supplying any argument returns the current
1374 This default can be overridden on a per-number basis by calling the
1375 C<display_format> method instead. As before, not supplying any argument
1376 returns the current display format for this number. Otherwise whatever you
1377 specify will be the new display format for I<this> particular number.
1383 Math::Complex::display_format('polar');
1384 $j = ((root(1, 3))[1];
1385 print "j = $j\n"; # Prints "j = [1,2pi/3]
1386 $j->display_format('cartesian');
1387 print "j = $j\n"; # Prints "j = -0.5+0.866025403784439i"
1389 The polar format attempts to emphasize arguments like I<k*pi/n>
1390 (where I<n> is a positive integer and I<k> an integer within [-9,+9]).
1394 Thanks to overloading, the handling of arithmetics with complex numbers
1395 is simple and almost transparent.
1397 Here are some examples:
1401 $j = cplxe(1, 2*pi/3); # $j ** 3 == 1
1402 print "j = $j, j**3 = ", $j ** 3, "\n";
1403 print "1 + j + j**2 = ", 1 + $j + $j**2, "\n";
1405 $z = -16 + 0*i; # Force it to be a complex
1406 print "sqrt($z) = ", sqrt($z), "\n";
1408 $k = exp(i * 2*pi/3);
1409 print "$j - $k = ", $j - $k, "\n";
1411 =head1 ERRORS DUE TO DIVISION BY ZERO
1413 The division (/) and the following functions
1432 cannot be computed for all arguments because that would mean dividing
1433 by zero. These situations cause fatal runtime errors looking like this
1435 cot(0): Division by zero.
1436 (Because in the definition of cot(0), the divisor sin(0) is 0)
1439 For the C<csc>, C<cot>, C<asec>, C<acsc>, C<csch>, C<coth>, C<asech>,
1440 C<acsch>, the argument cannot be C<0> (zero). For the C<atanh>,
1441 C<acoth>, the argument cannot be C<1> (one). For the C<atan>, C<acot>,
1442 the argument cannot be C<i> (the imaginary unit). For the C<tan>,
1443 C<sec>, C<tanh>, C<sech>, the argument cannot be I<pi/2 + k * pi>, where
1444 I<k> is any integer.
1448 Saying C<use Math::Complex;> exports many mathematical routines in the
1449 caller environment and even overrides some (C<sin>, C<cos>, C<sqrt>,
1450 C<log>, C<exp>). This is construed as a feature by the Authors,
1453 The code is not optimized for speed, although we try to use the cartesian
1454 form for addition-like operators and the trigonometric form for all
1455 multiplication-like operators.
1457 The arg() routine does not ensure the angle is within the range [-pi,+pi]
1458 (a side effect caused by multiplication and division using the trigonometric
1461 All routines expect to be given real or complex numbers. Don't attempt to
1462 use BigFloat, since Perl has currently no rule to disambiguate a '+'
1463 operation (for instance) between two overloaded entities.
1467 Raphael Manfredi <F<Raphael_Manfredi@grenoble.hp.com>>
1468 Jarkko Hietaniemi <F<jhi@iki.fi>>