3 # Complex numbers and associated mathematical functions
4 # -- Raphael Manfredi, September 1996
5 # -- Jarkko Hietaniemi, March 1997
12 use vars qw($VERSION @ISA
15 $pi $i $ilog10 $logn %logn);
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 $pi = 4 * atan2(1, 1) unless $pi;
146 # The number defined as i*i = -1;
149 $i = bless {} unless $i; # There can be only one i
150 $i->{'cartesian'} = [0, 1];
151 $i->{'polar'} = [1, pi/2];
158 # Attribute access/set routines
161 sub cartesian {$_[0]->{c_dirty} ?
162 $_[0]->update_cartesian : $_[0]->{'cartesian'}}
163 sub polar {$_[0]->{p_dirty} ?
164 $_[0]->update_polar : $_[0]->{'polar'}}
166 sub set_cartesian { $_[0]->{p_dirty}++; $_[0]->{'cartesian'} = $_[1] }
167 sub set_polar { $_[0]->{c_dirty}++; $_[0]->{'polar'} = $_[1] }
172 # Recompute and return the cartesian form, given accurate polar form.
174 sub update_cartesian {
176 my ($r, $t) = @{$self->{'polar'}};
177 $self->{c_dirty} = 0;
178 return $self->{'cartesian'} = [$r * cos $t, $r * sin $t];
185 # Recompute and return the polar form, given accurate cartesian form.
189 my ($x, $y) = @{$self->{'cartesian'}};
190 $self->{p_dirty} = 0;
191 return $self->{'polar'} = [0, 0] if $x == 0 && $y == 0;
192 return $self->{'polar'} = [sqrt($x*$x + $y*$y), atan2($y, $x)];
201 my ($z1, $z2, $regular) = @_;
202 $z2 = cplx($z2, 0) unless ref $z2;
203 my ($re1, $im1) = @{$z1->cartesian};
204 my ($re2, $im2) = @{$z2->cartesian};
205 unless (defined $regular) {
206 $z1->set_cartesian([$re1 + $re2, $im1 + $im2]);
209 return (ref $z1)->make($re1 + $re2, $im1 + $im2);
218 my ($z1, $z2, $inverted) = @_;
219 $z2 = cplx($z2, 0) unless ref $z2;
220 my ($re1, $im1) = @{$z1->cartesian};
221 my ($re2, $im2) = @{$z2->cartesian};
222 unless (defined $inverted) {
223 $z1->set_cartesian([$re1 - $re2, $im1 - $im2]);
227 (ref $z1)->make($re2 - $re1, $im2 - $im1) :
228 (ref $z1)->make($re1 - $re2, $im1 - $im2);
237 my ($z1, $z2, $regular) = @_;
238 my ($r1, $t1) = @{$z1->polar};
239 my ($r2, $t2) = ref $z2 ?
240 @{$z2->polar} : (abs($z2), $z2 >= 0 ? 0 : pi);
241 unless (defined $regular) {
242 $z1->set_polar([$r1 * $r2, $t1 + $t2]);
245 return (ref $z1)->emake($r1 * $r2, $t1 + $t2);
251 # Die on division by zero.
254 warn "$_[0]: Division by zero.\n";
255 warn "(Because in the definition of $_[0], $_[1] is 0)\n"
258 my $dmess = "Died at $up[1] line $up[2].\n";
268 my ($z1, $z2, $inverted) = @_;
269 my ($r1, $t1) = @{$z1->polar};
270 my ($r2, $t2) = ref $z2 ?
271 @{$z2->polar} : (abs($z2), $z2 >= 0 ? 0 : pi);
272 unless (defined $inverted) {
273 divbyzero "$z1/0" if ($r2 == 0);
274 $z1->set_polar([$r1 / $r2, $t1 - $t2]);
278 divbyzero "$z2/0" if ($r1 == 0);
279 return (ref $z1)->emake($r2 / $r1, $t2 - $t1);
281 divbyzero "$z1/0" if ($r2 == 0);
282 return (ref $z1)->emake($r1 / $r2, $t1 - $t2);
289 # Computes z1**z2 = exp(z2 * log z1)).
292 my ($z1, $z2, $inverted) = @_;
293 return exp($z1 * log $z2) if defined $inverted && $inverted;
294 return exp($z2 * log $z1);
300 # Computes z1 <=> z2.
301 # Sorts on the real part first, then on the imaginary part. Thus 2-4i > 3+8i.
304 my ($z1, $z2, $inverted) = @_;
305 $z2 = cplx($z2, 0) unless ref $z2;
306 my ($re1, $im1) = @{$z1->cartesian};
307 my ($re2, $im2) = @{$z2->cartesian};
308 my $sgn = $inverted ? -1 : 1;
309 return $sgn * ($re1 <=> $re2) if $re1 != $re2;
310 return $sgn * ($im1 <=> $im2);
321 my ($r, $t) = @{$z->polar};
322 return (ref $z)->emake($r, pi + $t);
324 my ($re, $im) = @{$z->cartesian};
325 return (ref $z)->make(-$re, -$im);
331 # Compute complex's conjugate.
336 my ($r, $t) = @{$z->polar};
337 return (ref $z)->emake($r, -$t);
339 my ($re, $im) = @{$z->cartesian};
340 return (ref $z)->make($re, -$im);
346 # Compute complex's norm (rho).
350 return abs($z) unless ref $z;
351 my ($r, $t) = @{$z->polar};
358 # Compute complex's argument (theta).
362 return ($z < 0 ? pi : 0) unless ref $z;
363 my ($r, $t) = @{$z->polar};
374 $z = cplx($z, 0) unless ref $z;
375 my ($r, $t) = @{$z->polar};
376 return (ref $z)->emake(sqrt($r), $t/2);
382 # Compute cbrt(z) (cubic root).
386 return cplx($z, 0) ** (1/3) unless ref $z;
387 my ($r, $t) = @{$z->polar};
388 return (ref $z)->emake($r**(1/3), $t/3);
394 # Computes all nth root for z, returning an array whose size is n.
395 # `n' must be a positive integer.
397 # The roots are given by (for k = 0..n-1):
399 # z^(1/n) = r^(1/n) (cos ((t+2 k pi)/n) + i sin ((t+2 k pi)/n))
404 return undef unless $n > 0;
405 my ($r, $t) = ref $z ? @{$z->polar} : (abs($z), $z >= 0 ? 0 : pi);
408 my $theta_inc = 2 * pi / $n;
409 my $rho = $r ** (1/$n);
411 my $complex = ref($z) || $package;
412 for ($k = 0, $theta = $t / $n; $k < $n; $k++, $theta += $theta_inc) {
413 push(@root, $complex->emake($rho, $theta));
425 return $z unless ref $z;
426 my ($re, $im) = @{$z->cartesian};
437 return 0 unless ref $z;
438 my ($re, $im) = @{$z->cartesian};
449 $z = cplx($z, 0) unless ref $z;
450 my ($x, $y) = @{$z->cartesian};
451 return (ref $z)->emake(exp($x), $y);
461 $z = cplx($z, 0) unless ref $z;
462 my ($r, $t) = @{$z->polar};
463 my ($x, $y) = @{$z->cartesian};
464 $t -= 2 * pi if ($t > pi() and $x < 0);
465 $t += 2 * pi if ($t < -pi() and $x < 0);
466 return (ref $z)->make(log($r), $t);
474 sub ln { Math::Complex::log(@_) }
483 my $ilog10 = 1 / log(10) unless defined $ilog10;
484 return log(cplx($z, 0)) * $ilog10 unless ref $z;
485 my ($r, $t) = @{$z->polar};
486 return (ref $z)->make(log($r) * $ilog10, $t * $ilog10);
492 # Compute logn(z,n) = log(z) / log(n)
496 $z = cplx($z, 0) unless ref $z;
497 my $logn = $logn{$n};
498 $logn = $logn{$n} = log($n) unless defined $logn; # Cache log(n)
499 return log($z) / $logn;
505 # Compute cos(z) = (exp(iz) + exp(-iz))/2.
509 my ($x, $y) = @{$z->cartesian};
512 return (ref $z)->make(cos($x) * ($ey + $ey_1)/2,
513 sin($x) * ($ey_1 - $ey)/2);
519 # Compute sin(z) = (exp(iz) - exp(-iz))/2.
523 my ($x, $y) = @{$z->cartesian};
526 return (ref $z)->make(sin($x) * ($ey + $ey_1)/2,
527 cos($x) * ($ey - $ey_1)/2);
533 # Compute tan(z) = sin(z) / cos(z).
538 divbyzero "tan($z)", "cos($z)" if ($cz == 0);
539 return sin($z) / $cz;
545 # Computes the secant sec(z) = 1 / cos(z).
550 divbyzero "sec($z)", "cos($z)" if ($cz == 0);
557 # Computes the cosecant csc(z) = 1 / sin(z).
562 divbyzero "csc($z)", "sin($z)" if ($sz == 0);
571 sub cosec { Math::Complex::csc(@_) }
576 # Computes cot(z) = 1 / tan(z).
581 divbyzero "cot($z)", "sin($z)" if ($sz == 0);
582 return cos($z) / $sz;
590 sub cotan { Math::Complex::cot(@_) }
595 # Computes the arc cosine acos(z) = -i log(z + sqrt(z*z-1)).
599 $z = cplx($z, 0) unless ref $z;
600 return ~i * log($z + (Re($z) * Im($z) > 0 ? 1 : -1) * sqrt($z*$z - 1));
606 # Computes the arc sine asin(z) = -i log(iz + sqrt(1-z*z)).
610 $z = cplx($z, 0) unless ref $z;
611 return ~i * log(i * $z + sqrt(1 - $z*$z));
617 # Computes the arc tangent atan(z) = i/2 log((i+z) / (i-z)).
621 divbyzero "atan($z)", "i - $z" if ($z == i);
622 return i/2*log((i + $z) / (i - $z));
628 # Computes the arc secant asec(z) = acos(1 / z).
638 # Computes the arc cosecant sec(z) = asin(1 / z).
648 # Alias for acosec().
650 sub acsc { Math::Complex::acosec(@_) }
655 # Computes the arc cotangent acot(z) = -i/2 log((i+z) / (z-i))
659 divbyzero "acot($z)", "$z - i" if ($z == i);
660 return i/-2 * log((i + $z) / ($z - i));
668 sub acotan { Math::Complex::acot(@_) }
673 # Computes the hyperbolic cosine cosh(z) = (exp(z) + exp(-z))/2.
677 $z = cplx($z, 0) unless ref $z;
678 my ($x, $y) = @{$z->cartesian};
681 return ($ex + $ex_1)/2 unless ref $z;
682 return (ref $z)->make(cos($y) * ($ex + $ex_1)/2,
683 sin($y) * ($ex - $ex_1)/2);
689 # Computes the hyperbolic sine sinh(z) = (exp(z) - exp(-z))/2.
693 $z = cplx($z, 0) unless ref $z;
694 my ($x, $y) = @{$z->cartesian};
697 return ($ex - $ex_1)/2 unless ref $z;
698 return (ref $z)->make(cos($y) * ($ex - $ex_1)/2,
699 sin($y) * ($ex + $ex_1)/2);
705 # Computes the hyperbolic tangent tanh(z) = sinh(z) / cosh(z).
710 divbyzero "tanh($z)", "cosh($z)" if ($cz == 0);
711 return sinh($z) / $cz;
717 # Computes the hyperbolic secant sech(z) = 1 / cosh(z).
722 divbyzero "sech($z)", "cosh($z)" if ($cz == 0);
729 # Computes the hyperbolic cosecant csch(z) = 1 / sinh(z).
734 divbyzero "csch($z)", "sinh($z)" if ($sz == 0);
743 sub cosech { Math::Complex::csch(@_) }
748 # Computes the hyperbolic cotangent coth(z) = cosh(z) / sinh(z).
753 divbyzero "coth($z)", "sinh($z)" if ($sz == 0);
754 return cosh($z) / $sz;
762 sub cotanh { Math::Complex::coth(@_) }
767 # Computes the arc hyperbolic cosine acosh(z) = log(z + sqrt(z*z-1)).
771 $z = cplx($z, 0) unless ref $z; # asinh(-2)
772 return log($z + sqrt($z*$z - 1));
778 # Computes the arc hyperbolic sine asinh(z) = log(z + sqrt(z*z-1))
782 $z = cplx($z, 0) unless ref $z; # asinh(-2)
783 return log($z + sqrt($z*$z + 1));
789 # Computes the arc hyperbolic tangent atanh(z) = 1/2 log((1+z) / (1-z)).
793 $z = cplx($z, 0) unless ref $z; # atanh(-2)
794 divbyzero 'atanh(1)', "1 - $z" if ($z == 1);
795 my $cz = (1 + $z) / (1 - $z);
802 # Computes the hyperbolic arc secant asech(z) = acosh(1 / z).
806 divbyzero 'asech(0)', $z if ($z == 0);
807 return acosh(1 / $z);
813 # Computes the hyperbolic arc cosecant acsch(z) = asinh(1 / z).
817 divbyzero 'acsch(0)', $z if ($z == 0);
818 return asinh(1 / $z);
826 sub acosech { Math::Complex::acsch(@_) }
831 # Computes the arc hyperbolic cotangent acoth(z) = 1/2 log((1+z) / (z-1)).
835 $z = cplx($z, 0) unless ref $z; # acoth(-2)
836 divbyzero 'acoth(1)', "$z - 1" if ($z == 1);
837 my $cz = (1 + $z) / ($z - 1);
846 sub acotanh { Math::Complex::acoth(@_) }
851 # Compute atan(z1/z2).
854 my ($z1, $z2, $inverted) = @_;
855 my ($re1, $im1) = @{$z1->cartesian};
856 my ($re2, $im2) = @{$z2->cartesian};
858 if (defined $inverted && $inverted) { # atan(z2/z1)
859 return pi * ($re2 > 0 ? 1 : -1) if $re1 == 0 && $im1 == 0;
862 return pi * ($re1 > 0 ? 1 : -1) if $re2 == 0 && $im2 == 0;
872 # Set (fetch if no argument) display format for all complex numbers that
873 # don't happen to have overrriden it via ->display_format
875 # When called as a method, this actually sets the display format for
876 # the current object.
878 # Valid object formats are 'c' and 'p' for cartesian and polar. The first
879 # letter is used actually, so the type can be fully spelled out for clarity.
885 if (ref $self) { # Called as a method
887 } else { # Regular procedure call
893 return defined $self->{display} ? $self->{display} : $display
894 unless defined $format;
895 return $self->{display} = $format;
898 return $display unless defined $format;
899 return $display = $format;
905 # Show nicely formatted complex number under its cartesian or polar form,
906 # depending on the current display format:
908 # . If a specific display format has been recorded for this object, use it.
909 # . Otherwise, use the generic current default for all complex numbers,
910 # which is a package global variable.
917 $format = $z->{display} if defined $z->{display};
919 return $z->stringify_polar if $format =~ /^p/i;
920 return $z->stringify_cartesian;
924 # ->stringify_cartesian
926 # Stringify as a cartesian representation 'a+bi'.
928 sub stringify_cartesian {
930 my ($x, $y) = @{$z->cartesian};
933 $x = int($x + ($x < 0 ? -1 : 1) * 1e-14)
934 if int(abs($x)) != int(abs($x) + 1e-14);
935 $y = int($y + ($y < 0 ? -1 : 1) * 1e-14)
936 if int(abs($y)) != int(abs($y) + 1e-14);
938 $re = "$x" if abs($x) >= 1e-14;
939 if ($y == 1) { $im = 'i' }
940 elsif ($y == -1) { $im = '-i' }
941 elsif (abs($y) >= 1e-14) { $im = $y . "i" }
944 $str = $re if defined $re;
945 $str .= "+$im" if defined $im;
948 $str = '0' unless $str;
956 # Stringify as a polar representation '[r,t]'.
958 sub stringify_polar {
960 my ($r, $t) = @{$z->polar};
964 return '[0,0]' if $r <= $eps;
968 $nt = ($nt - int($nt)) * $tpi;
969 $nt += $tpi if $nt < 0; # Range [0, 2pi]
971 if (abs($nt) <= $eps) { $theta = 0 }
972 elsif (abs(pi-$nt) <= $eps) { $theta = 'pi' }
974 if (defined $theta) {
975 $r = int($r + ($r < 0 ? -1 : 1) * $eps)
976 if int(abs($r)) != int(abs($r) + $eps);
977 $theta = int($theta + ($theta < 0 ? -1 : 1) * $eps)
978 if ($theta ne 'pi' and
979 int(abs($theta)) != int(abs($theta) + $eps));
980 return "\[$r,$theta\]";
984 # Okay, number is not a real. Try to identify pi/n and friends...
987 $nt -= $tpi if $nt > pi;
990 for ($k = 1, $kpi = pi; $k < 10; $k++, $kpi += pi) {
991 $n = int($kpi / $nt + ($nt > 0 ? 1 : -1) * 0.5);
992 if (abs($kpi/$n - $nt) <= $eps) {
993 $theta = ($nt < 0 ? '-':'').
994 ($k == 1 ? 'pi':"${k}pi").'/'.abs($n);
999 $theta = $nt unless defined $theta;
1001 $r = int($r + ($r < 0 ? -1 : 1) * $eps)
1002 if int(abs($r)) != int(abs($r) + $eps);
1003 $theta = int($theta + ($theta < 0 ? -1 : 1) * $eps)
1004 if ($theta !~ m(^-?\d*pi/\d+$) and
1005 int(abs($theta)) != int(abs($theta) + $eps));
1007 return "\[$r,$theta\]";
1015 Math::Complex - complex numbers and associated mathematical functions
1021 $z = Math::Complex->make(5, 6);
1023 $j = cplxe(1, 2*pi/3);
1027 This package lets you create and manipulate complex numbers. By default,
1028 I<Perl> limits itself to real numbers, but an extra C<use> statement brings
1029 full complex support, along with a full set of mathematical functions
1030 typically associated with and/or extended to complex numbers.
1032 If you wonder what complex numbers are, they were invented to be able to solve
1033 the following equation:
1037 and by definition, the solution is noted I<i> (engineers use I<j> instead since
1038 I<i> usually denotes an intensity, but the name does not matter). The number
1039 I<i> is a pure I<imaginary> number.
1041 The arithmetics with pure imaginary numbers works just like you would expect
1042 it with real numbers... you just have to remember that
1048 5i + 7i = i * (5 + 7) = 12i
1049 4i - 3i = i * (4 - 3) = i
1054 Complex numbers are numbers that have both a real part and an imaginary
1055 part, and are usually noted:
1059 where C<a> is the I<real> part and C<b> is the I<imaginary> part. The
1060 arithmetic with complex numbers is straightforward. You have to
1061 keep track of the real and the imaginary parts, but otherwise the
1062 rules used for real numbers just apply:
1064 (4 + 3i) + (5 - 2i) = (4 + 5) + i(3 - 2) = 9 + i
1065 (2 + i) * (4 - i) = 2*4 + 4i -2i -i*i = 8 + 2i + 1 = 9 + 2i
1067 A graphical representation of complex numbers is possible in a plane
1068 (also called the I<complex plane>, but it's really a 2D plane).
1073 is the point whose coordinates are (a, b). Actually, it would
1074 be the vector originating from (0, 0) to (a, b). It follows that the addition
1075 of two complex numbers is a vectorial addition.
1077 Since there is a bijection between a point in the 2D plane and a complex
1078 number (i.e. the mapping is unique and reciprocal), a complex number
1079 can also be uniquely identified with polar coordinates:
1083 where C<rho> is the distance to the origin, and C<theta> the angle between
1084 the vector and the I<x> axis. There is a notation for this using the
1085 exponential form, which is:
1087 rho * exp(i * theta)
1089 where I<i> is the famous imaginary number introduced above. Conversion
1090 between this form and the cartesian form C<a + bi> is immediate:
1092 a = rho * cos(theta)
1093 b = rho * sin(theta)
1095 which is also expressed by this formula:
1097 z = rho * exp(i * theta) = rho * (cos theta + i * sin theta)
1099 In other words, it's the projection of the vector onto the I<x> and I<y>
1100 axes. Mathematicians call I<rho> the I<norm> or I<modulus> and I<theta>
1101 the I<argument> of the complex number. The I<norm> of C<z> will be
1104 The polar notation (also known as the trigonometric
1105 representation) is much more handy for performing multiplications and
1106 divisions of complex numbers, whilst the cartesian notation is better
1107 suited for additions and substractions. Real numbers are on the I<x>
1108 axis, and therefore I<theta> is zero.
1110 All the common operations that can be performed on a real number have
1111 been defined to work on complex numbers as well, and are merely
1112 I<extensions> of the operations defined on real numbers. This means
1113 they keep their natural meaning when there is no imaginary part, provided
1114 the number is within their definition set.
1116 For instance, the C<sqrt> routine which computes the square root of
1117 its argument is only defined for positive real numbers and yields a
1118 positive real number (it is an application from B<R+> to B<R+>).
1119 If we allow it to return a complex number, then it can be extended to
1120 negative real numbers to become an application from B<R> to B<C> (the
1121 set of complex numbers):
1123 sqrt(x) = x >= 0 ? sqrt(x) : sqrt(-x)*i
1125 It can also be extended to be an application from B<C> to B<C>,
1126 whilst its restriction to B<R> behaves as defined above by using
1127 the following definition:
1129 sqrt(z = [r,t]) = sqrt(r) * exp(i * t/2)
1131 Indeed, a negative real number can be noted C<[x,pi]>
1132 (the modulus I<x> is always positive, so C<[x,pi]> is really C<-x>, a
1134 and the above definition states that
1136 sqrt([x,pi]) = sqrt(x) * exp(i*pi/2) = [sqrt(x),pi/2] = sqrt(x)*i
1138 which is exactly what we had defined for negative real numbers above.
1140 All the common mathematical functions defined on real numbers that
1141 are extended to complex numbers share that same property of working
1142 I<as usual> when the imaginary part is zero (otherwise, it would not
1143 be called an extension, would it?).
1145 A I<new> operation possible on a complex number that is
1146 the identity for real numbers is called the I<conjugate>, and is noted
1147 with an horizontal bar above the number, or C<~z> here.
1154 z * ~z = (a + bi) * (a - bi) = a*a + b*b
1156 We saw that the norm of C<z> was noted C<abs(z)> and was defined as the
1157 distance to the origin, also known as:
1159 rho = abs(z) = sqrt(a*a + b*b)
1163 z * ~z = abs(z) ** 2
1165 If z is a pure real number (i.e. C<b == 0>), then the above yields:
1169 which is true (C<abs> has the regular meaning for real number, i.e. stands
1170 for the absolute value). This example explains why the norm of C<z> is
1171 noted C<abs(z)>: it extends the C<abs> function to complex numbers, yet
1172 is the regular C<abs> we know when the complex number actually has no
1173 imaginary part... This justifies I<a posteriori> our use of the C<abs>
1174 notation for the norm.
1178 Given the following notations:
1180 z1 = a + bi = r1 * exp(i * t1)
1181 z2 = c + di = r2 * exp(i * t2)
1182 z = <any complex or real number>
1184 the following (overloaded) operations are supported on complex numbers:
1186 z1 + z2 = (a + c) + i(b + d)
1187 z1 - z2 = (a - c) + i(b - d)
1188 z1 * z2 = (r1 * r2) * exp(i * (t1 + t2))
1189 z1 / z2 = (r1 / r2) * exp(i * (t1 - t2))
1190 z1 ** z2 = exp(z2 * log z1)
1192 abs(z1) = r1 = sqrt(a*a + b*b)
1193 sqrt(z1) = sqrt(r1) * exp(i * t1/2)
1194 exp(z1) = exp(a) * exp(i * b)
1195 log(z1) = log(r1) + i*t1
1196 sin(z1) = 1/2i (exp(i * z1) - exp(-i * z1))
1197 cos(z1) = 1/2 (exp(i * z1) + exp(-i * z1))
1199 atan2(z1, z2) = atan(z1/z2)
1201 The following extra operations are supported on both real and complex
1208 cbrt(z) = z ** (1/3)
1209 log10(z) = log(z) / log(10)
1210 logn(z, n) = log(z) / log(n)
1212 tan(z) = sin(z) / cos(z)
1218 asin(z) = -i * log(i*z + sqrt(1-z*z))
1219 acos(z) = -i * log(z + sqrt(z*z-1))
1220 atan(z) = i/2 * log((i+z) / (i-z))
1222 acsc(z) = asin(1 / z)
1223 asec(z) = acos(1 / z)
1224 acot(z) = -i/2 * log((i+z) / (z-i))
1226 sinh(z) = 1/2 (exp(z) - exp(-z))
1227 cosh(z) = 1/2 (exp(z) + exp(-z))
1228 tanh(z) = sinh(z) / cosh(z) = (exp(z) - exp(-z)) / (exp(z) + exp(-z))
1230 csch(z) = 1 / sinh(z)
1231 sech(z) = 1 / cosh(z)
1232 coth(z) = 1 / tanh(z)
1234 asinh(z) = log(z + sqrt(z*z+1))
1235 acosh(z) = log(z + sqrt(z*z-1))
1236 atanh(z) = 1/2 * log((1+z) / (1-z))
1238 acsch(z) = asinh(1 / z)
1239 asech(z) = acosh(1 / z)
1240 acoth(z) = atanh(1 / z) = 1/2 * log((1+z) / (z-1))
1242 I<log>, I<csc>, I<cot>, I<acsc>, I<acot>, I<csch>, I<coth>,
1243 I<acosech>, I<acotanh>, have aliases I<ln>, I<cosec>, I<cotan>,
1244 I<acosec>, I<acotan>, I<cosech>, I<cotanh>, I<acosech>, I<acotanh>,
1247 The I<root> function is available to compute all the I<n>
1248 roots of some complex, where I<n> is a strictly positive integer.
1249 There are exactly I<n> such roots, returned as a list. Getting the
1250 number mathematicians call C<j> such that:
1254 is a simple matter of writing:
1256 $j = ((root(1, 3))[1];
1258 The I<k>th root for C<z = [r,t]> is given by:
1260 (root(z, n))[k] = r**(1/n) * exp(i * (t + 2*k*pi)/n)
1262 The I<spaceship> comparison operator is also defined. In order to
1263 ensure its restriction to real numbers is conform to what you would
1264 expect, the comparison is run on the real part of the complex number
1265 first, and imaginary parts are compared only when the real parts
1270 To create a complex number, use either:
1272 $z = Math::Complex->make(3, 4);
1275 if you know the cartesian form of the number, or
1279 if you like. To create a number using the trigonometric form, use either:
1281 $z = Math::Complex->emake(5, pi/3);
1282 $x = cplxe(5, pi/3);
1284 instead. The first argument is the modulus, the second is the angle
1285 (in radians, the full circle is 2*pi). (Mnmemonic: C<e> is used as a
1286 notation for complex numbers in the trigonometric form).
1288 It is possible to write:
1290 $x = cplxe(-3, pi/4);
1292 but that will be silently converted into C<[3,-3pi/4]>, since the modulus
1293 must be positive (it represents the distance to the origin in the complex
1296 =head1 STRINGIFICATION
1298 When printed, a complex number is usually shown under its cartesian
1299 form I<a+bi>, but there are legitimate cases where the polar format
1300 I<[r,t]> is more appropriate.
1302 By calling the routine C<Math::Complex::display_format> and supplying either
1303 C<"polar"> or C<"cartesian">, you override the default display format,
1304 which is C<"cartesian">. Not supplying any argument returns the current
1307 This default can be overridden on a per-number basis by calling the
1308 C<display_format> method instead. As before, not supplying any argument
1309 returns the current display format for this number. Otherwise whatever you
1310 specify will be the new display format for I<this> particular number.
1316 Math::Complex::display_format('polar');
1317 $j = ((root(1, 3))[1];
1318 print "j = $j\n"; # Prints "j = [1,2pi/3]
1319 $j->display_format('cartesian');
1320 print "j = $j\n"; # Prints "j = -0.5+0.866025403784439i"
1322 The polar format attempts to emphasize arguments like I<k*pi/n>
1323 (where I<n> is a positive integer and I<k> an integer within [-9,+9]).
1327 Thanks to overloading, the handling of arithmetics with complex numbers
1328 is simple and almost transparent.
1330 Here are some examples:
1334 $j = cplxe(1, 2*pi/3); # $j ** 3 == 1
1335 print "j = $j, j**3 = ", $j ** 3, "\n";
1336 print "1 + j + j**2 = ", 1 + $j + $j**2, "\n";
1338 $z = -16 + 0*i; # Force it to be a complex
1339 print "sqrt($z) = ", sqrt($z), "\n";
1341 $k = exp(i * 2*pi/3);
1342 print "$j - $k = ", $j - $k, "\n";
1346 The division (/) and the following functions
1363 cannot be computed for all arguments because that would mean dividing
1364 by zero. These situations cause fatal runtime errors looking like this
1366 cot(0): Division by zero.
1367 (Because in the definition of cot(0), sin(0) is 0)
1372 Saying C<use Math::Complex;> exports many mathematical routines in the caller
1373 environment. This is construed as a feature by the Author, actually... ;-)
1375 The code is not optimized for speed, although we try to use the cartesian
1376 form for addition-like operators and the trigonometric form for all
1377 multiplication-like operators.
1379 The arg() routine does not ensure the angle is within the range [-pi,+pi]
1380 (a side effect caused by multiplication and division using the trigonometric
1383 All routines expect to be given real or complex numbers. Don't attempt to
1384 use BigFloat, since Perl has currently no rule to disambiguate a '+'
1385 operation (for instance) between two overloaded entities.
1389 Raphael Manfredi <F<Raphael_Manfredi@grenoble.hp.com>>
1390 Jarkko Hietaniemi <F<jhi@iki.fi>>