3 # Complex numbers and associated mathematical functions
4 # -- Raphael Manfredi, September 1996
5 # -- Jarkko Hietaniemi, March 1997
8 package Math::Complex; @ISA = qw(Exporter);
12 use vars qw(@EXPORT $package $display
13 $pi $i $ilog10 $logn %logn);
20 cosec csc sec cotan cot
22 acosec acsc asec acotan acot
24 cosech csch sech cotanh coth
26 acosech acsch asech acotanh acoth
53 $package = 'Math::Complex'; # Package name
54 $display = 'cartesian'; # Default display format
57 # Object attributes (internal):
58 # cartesian [real, imaginary] -- cartesian form
59 # polar [rho, theta] -- polar form
60 # c_dirty cartesian form not up-to-date
61 # p_dirty polar form not up-to-date
62 # display display format (package's global when not set)
68 # Create a new complex number (cartesian form)
71 my $self = bless {}, shift;
73 $self->{'cartesian'} = [$re, $im];
82 # Create a new complex number (exponential form)
85 my $self = bless {}, shift;
86 my ($rho, $theta) = @_;
87 $theta += pi() if $rho < 0;
88 $self->{'polar'} = [abs($rho), $theta];
94 sub new { &make } # For backward compatibility only.
99 # Creates a complex number from a (re, im) tuple.
100 # This avoids the burden of writing Math::Complex->make(re, im).
104 return $package->make($re, defined $im ? $im : 0);
110 # Creates a complex number from a (rho, theta) tuple.
111 # This avoids the burden of writing Math::Complex->emake(rho, theta).
114 my ($rho, $theta) = @_;
115 return $package->emake($rho, defined $theta ? $theta : 0);
121 # The number defined as 2 * pi = 360 degrees
124 $pi = 4 * atan2(1, 1) unless $pi;
131 # The number defined as i*i = -1;
134 $i = bless {} unless $i; # There can be only one i
135 $i->{'cartesian'} = [0, 1];
136 $i->{'polar'} = [1, pi/2];
143 # Attribute access/set routines
146 sub cartesian {$_[0]->{c_dirty} ?
147 $_[0]->update_cartesian : $_[0]->{'cartesian'}}
148 sub polar {$_[0]->{p_dirty} ?
149 $_[0]->update_polar : $_[0]->{'polar'}}
151 sub set_cartesian { $_[0]->{p_dirty}++; $_[0]->{'cartesian'} = $_[1] }
152 sub set_polar { $_[0]->{c_dirty}++; $_[0]->{'polar'} = $_[1] }
157 # Recompute and return the cartesian form, given accurate polar form.
159 sub update_cartesian {
161 my ($r, $t) = @{$self->{'polar'}};
162 $self->{c_dirty} = 0;
163 return $self->{'cartesian'} = [$r * cos $t, $r * sin $t];
170 # Recompute and return the polar form, given accurate cartesian form.
174 my ($x, $y) = @{$self->{'cartesian'}};
175 $self->{p_dirty} = 0;
176 return $self->{'polar'} = [0, 0] if $x == 0 && $y == 0;
177 return $self->{'polar'} = [sqrt($x*$x + $y*$y), atan2($y, $x)];
186 my ($z1, $z2, $regular) = @_;
187 $z2 = cplx($z2, 0) unless ref $z2;
188 my ($re1, $im1) = @{$z1->cartesian};
189 my ($re2, $im2) = @{$z2->cartesian};
190 unless (defined $regular) {
191 $z1->set_cartesian([$re1 + $re2, $im1 + $im2]);
194 return (ref $z1)->make($re1 + $re2, $im1 + $im2);
203 my ($z1, $z2, $inverted) = @_;
204 $z2 = cplx($z2, 0) unless ref $z2;
205 my ($re1, $im1) = @{$z1->cartesian};
206 my ($re2, $im2) = @{$z2->cartesian};
207 unless (defined $inverted) {
208 $z1->set_cartesian([$re1 - $re2, $im1 - $im2]);
212 (ref $z1)->make($re2 - $re1, $im2 - $im1) :
213 (ref $z1)->make($re1 - $re2, $im1 - $im2);
222 my ($z1, $z2, $regular) = @_;
223 my ($r1, $t1) = @{$z1->polar};
224 my ($r2, $t2) = ref $z2 ?
225 @{$z2->polar} : (abs($z2), $z2 >= 0 ? 0 : pi);
226 unless (defined $regular) {
227 $z1->set_polar([$r1 * $r2, $t1 + $t2]);
230 return (ref $z1)->emake($r1 * $r2, $t1 + $t2);
236 # Die on division by zero.
239 warn $package . '::' . "$_[0]: Division by zero.\n";
240 warn "(Because in the definition of $_[0], $_[1] is 0)\n"
243 my $dmess = "Died at $up[1] line $up[2].\n";
253 my ($z1, $z2, $inverted) = @_;
254 my ($r1, $t1) = @{$z1->polar};
255 my ($r2, $t2) = ref $z2 ?
256 @{$z2->polar} : (abs($z2), $z2 >= 0 ? 0 : pi);
257 unless (defined $inverted) {
258 divbyzero "$z1/0" if ($r2 == 0);
259 $z1->set_polar([$r1 / $r2, $t1 - $t2]);
263 divbyzero "$z2/0" if ($r1 == 0);
264 return (ref $z1)->emake($r2 / $r1, $t2 - $t1);
266 divbyzero "$z1/0" if ($r2 == 0);
267 return (ref $z1)->emake($r1 / $r2, $t1 - $t2);
274 # Computes z1**z2 = exp(z2 * log z1)).
277 my ($z1, $z2, $inverted) = @_;
278 return exp($z1 * log $z2) if defined $inverted && $inverted;
279 return exp($z2 * log $z1);
285 # Computes z1 <=> z2.
286 # Sorts on the real part first, then on the imaginary part. Thus 2-4i > 3+8i.
289 my ($z1, $z2, $inverted) = @_;
290 $z2 = cplx($z2, 0) unless ref $z2;
291 my ($re1, $im1) = @{$z1->cartesian};
292 my ($re2, $im2) = @{$z2->cartesian};
293 my $sgn = $inverted ? -1 : 1;
294 return $sgn * ($re1 <=> $re2) if $re1 != $re2;
295 return $sgn * ($im1 <=> $im2);
306 my ($r, $t) = @{$z->polar};
307 return (ref $z)->emake($r, pi + $t);
309 my ($re, $im) = @{$z->cartesian};
310 return (ref $z)->make(-$re, -$im);
316 # Compute complex's conjugate.
321 my ($r, $t) = @{$z->polar};
322 return (ref $z)->emake($r, -$t);
324 my ($re, $im) = @{$z->cartesian};
325 return (ref $z)->make($re, -$im);
331 # Compute complex's norm (rho).
335 return abs($z) unless ref $z;
336 my ($r, $t) = @{$z->polar};
343 # Compute complex's argument (theta).
347 return ($z < 0 ? pi : 0) unless ref $z;
348 my ($r, $t) = @{$z->polar};
359 $z = cplx($z, 0) unless ref $z;
360 my ($r, $t) = @{$z->polar};
361 return (ref $z)->emake(sqrt($r), $t/2);
367 # Compute cbrt(z) (cubic root).
371 return cplx($z, 0) ** (1/3) unless ref $z;
372 my ($r, $t) = @{$z->polar};
373 return (ref $z)->emake($r**(1/3), $t/3);
379 # Computes all nth root for z, returning an array whose size is n.
380 # `n' must be a positive integer.
382 # The roots are given by (for k = 0..n-1):
384 # z^(1/n) = r^(1/n) (cos ((t+2 k pi)/n) + i sin ((t+2 k pi)/n))
389 return undef unless $n > 0;
390 my ($r, $t) = ref $z ? @{$z->polar} : (abs($z), $z >= 0 ? 0 : pi);
393 my $theta_inc = 2 * pi / $n;
394 my $rho = $r ** (1/$n);
396 my $complex = ref($z) || $package;
397 for ($k = 0, $theta = $t / $n; $k < $n; $k++, $theta += $theta_inc) {
398 push(@root, $complex->emake($rho, $theta));
410 return $z unless ref $z;
411 my ($re, $im) = @{$z->cartesian};
422 return 0 unless ref $z;
423 my ($re, $im) = @{$z->cartesian};
434 $z = cplx($z, 0) unless ref $z;
435 my ($x, $y) = @{$z->cartesian};
436 return (ref $z)->emake(exp($x), $y);
446 $z = cplx($z, 0) unless ref $z;
447 my ($r, $t) = @{$z->polar};
448 my ($x, $y) = @{$z->cartesian};
449 $t -= 2 * pi if ($t > pi() and $x < 0);
450 $t += 2 * pi if ($t < -pi() and $x < 0);
451 return (ref $z)->make(log($r), $t);
459 sub ln { Math::Complex::log(@_) }
468 my $ilog10 = 1 / log(10) unless defined $ilog10;
469 return log(cplx($z, 0)) * $ilog10 unless ref $z;
470 my ($r, $t) = @{$z->polar};
471 return (ref $z)->make(log($r) * $ilog10, $t * $ilog10);
477 # Compute logn(z,n) = log(z) / log(n)
481 $z = cplx($z, 0) unless ref $z;
482 my $logn = $logn{$n};
483 $logn = $logn{$n} = log($n) unless defined $logn; # Cache log(n)
484 return log($z) / $logn;
490 # Compute cos(z) = (exp(iz) + exp(-iz))/2.
494 my ($x, $y) = @{$z->cartesian};
497 return (ref $z)->make(cos($x) * ($ey + $ey_1)/2,
498 sin($x) * ($ey_1 - $ey)/2);
504 # Compute sin(z) = (exp(iz) - exp(-iz))/2.
508 my ($x, $y) = @{$z->cartesian};
511 return (ref $z)->make(sin($x) * ($ey + $ey_1)/2,
512 cos($x) * ($ey - $ey_1)/2);
518 # Compute tan(z) = sin(z) / cos(z).
523 divbyzero "tan($z)", "cos($z)" if ($cz == 0);
524 return sin($z) / $cz;
530 # Computes the secant sec(z) = 1 / cos(z).
535 divbyzero "sec($z)", "cos($z)" if ($cz == 0);
542 # Computes the cosecant csc(z) = 1 / sin(z).
547 divbyzero "csc($z)", "sin($z)" if ($sz == 0);
556 sub cosec { Math::Complex::csc(@_) }
561 # Computes cot(z) = 1 / tan(z).
566 divbyzero "cot($z)", "sin($z)" if ($sz == 0);
567 return cos($z) / $sz;
575 sub cotan { Math::Complex::cot(@_) }
580 # Computes the arc cosine acos(z) = -i log(z + sqrt(z*z-1)).
584 $z = cplx($z, 0) unless ref $z;
585 return ~i * log($z + (Re($z) * Im($z) > 0 ? 1 : -1) * sqrt($z*$z - 1));
591 # Computes the arc sine asin(z) = -i log(iz + sqrt(1-z*z)).
595 $z = cplx($z, 0) unless ref $z;
596 return ~i * log(i * $z + sqrt(1 - $z*$z));
602 # Computes the arc tangent atan(z) = i/2 log((i+z) / (i-z)).
606 divbyzero "atan($z)", "i - $z" if ($z == i);
607 return i/2*log((i + $z) / (i - $z));
613 # Computes the arc secant asec(z) = acos(1 / z).
623 # Computes the arc cosecant sec(z) = asin(1 / z).
633 # Alias for acosec().
635 sub acsc { Math::Complex::acosec(@_) }
640 # Computes the arc cotangent acot(z) = -i/2 log((i+z) / (z-i))
644 divbyzero "acot($z)", "$z - i" if ($z == i);
645 return i/-2 * log((i + $z) / ($z - i));
653 sub acotan { Math::Complex::acot(@_) }
658 # Computes the hyperbolic cosine cosh(z) = (exp(z) + exp(-z))/2.
662 $z = cplx($z, 0) unless ref $z;
663 my ($x, $y) = @{$z->cartesian};
666 return ($ex + $ex_1)/2 unless ref $z;
667 return (ref $z)->make(cos($y) * ($ex + $ex_1)/2,
668 sin($y) * ($ex - $ex_1)/2);
674 # Computes the hyperbolic sine sinh(z) = (exp(z) - exp(-z))/2.
678 $z = cplx($z, 0) unless ref $z;
679 my ($x, $y) = @{$z->cartesian};
682 return ($ex - $ex_1)/2 unless ref $z;
683 return (ref $z)->make(cos($y) * ($ex - $ex_1)/2,
684 sin($y) * ($ex + $ex_1)/2);
690 # Computes the hyperbolic tangent tanh(z) = sinh(z) / cosh(z).
695 divbyzero "tanh($z)", "cosh($z)" if ($cz == 0);
696 return sinh($z) / $cz;
702 # Computes the hyperbolic secant sech(z) = 1 / cosh(z).
707 divbyzero "sech($z)", "cosh($z)" if ($cz == 0);
714 # Computes the hyperbolic cosecant csch(z) = 1 / sinh(z).
719 divbyzero "csch($z)", "sinh($z)" if ($sz == 0);
728 sub cosech { Math::Complex::csch(@_) }
733 # Computes the hyperbolic cotangent coth(z) = cosh(z) / sinh(z).
738 divbyzero "coth($z)", "sinh($z)" if ($sz == 0);
739 return cosh($z) / $sz;
747 sub cotanh { Math::Complex::coth(@_) }
752 # Computes the arc hyperbolic cosine acosh(z) = log(z + sqrt(z*z-1)).
756 $z = cplx($z, 0) unless ref $z; # asinh(-2)
757 return log($z + sqrt($z*$z - 1));
763 # Computes the arc hyperbolic sine asinh(z) = log(z + sqrt(z*z-1))
767 $z = cplx($z, 0) unless ref $z; # asinh(-2)
768 return log($z + sqrt($z*$z + 1));
774 # Computes the arc hyperbolic tangent atanh(z) = 1/2 log((1+z) / (1-z)).
778 $z = cplx($z, 0) unless ref $z; # atanh(-2)
779 divbyzero 'atanh(1)', "1 - $z" if ($z == 1);
780 my $cz = (1 + $z) / (1 - $z);
787 # Computes the hyperbolic arc secant asech(z) = acosh(1 / z).
791 divbyzero 'asech(0)', $z if ($z == 0);
792 return acosh(1 / $z);
798 # Computes the hyperbolic arc cosecant acsch(z) = asinh(1 / z).
802 divbyzero 'acsch(0)', $z if ($z == 0);
803 return asinh(1 / $z);
811 sub acosech { Math::Complex::acsch(@_) }
816 # Computes the arc hyperbolic cotangent acoth(z) = 1/2 log((1+z) / (z-1)).
820 $z = cplx($z, 0) unless ref $z; # acoth(-2)
821 divbyzero 'acoth(1)', "$z - 1" if ($z == 1);
822 my $cz = (1 + $z) / ($z - 1);
831 sub acotanh { Math::Complex::acoth(@_) }
836 # Compute atan(z1/z2).
839 my ($z1, $z2, $inverted) = @_;
840 my ($re1, $im1) = @{$z1->cartesian};
841 my ($re2, $im2) = @{$z2->cartesian};
843 if (defined $inverted && $inverted) { # atan(z2/z1)
844 return pi * ($re2 > 0 ? 1 : -1) if $re1 == 0 && $im1 == 0;
847 return pi * ($re1 > 0 ? 1 : -1) if $re2 == 0 && $im2 == 0;
857 # Set (fetch if no argument) display format for all complex numbers that
858 # don't happen to have overrriden it via ->display_format
860 # When called as a method, this actually sets the display format for
861 # the current object.
863 # Valid object formats are 'c' and 'p' for cartesian and polar. The first
864 # letter is used actually, so the type can be fully spelled out for clarity.
870 if (ref $self) { # Called as a method
872 } else { # Regular procedure call
878 return defined $self->{display} ? $self->{display} : $display
879 unless defined $format;
880 return $self->{display} = $format;
883 return $display unless defined $format;
884 return $display = $format;
890 # Show nicely formatted complex number under its cartesian or polar form,
891 # depending on the current display format:
893 # . If a specific display format has been recorded for this object, use it.
894 # . Otherwise, use the generic current default for all complex numbers,
895 # which is a package global variable.
902 $format = $z->{display} if defined $z->{display};
904 return $z->stringify_polar if $format =~ /^p/i;
905 return $z->stringify_cartesian;
909 # ->stringify_cartesian
911 # Stringify as a cartesian representation 'a+bi'.
913 sub stringify_cartesian {
915 my ($x, $y) = @{$z->cartesian};
918 $x = int($x + ($x < 0 ? -1 : 1) * 1e-14)
919 if int(abs($x)) != int(abs($x) + 1e-14);
920 $y = int($y + ($y < 0 ? -1 : 1) * 1e-14)
921 if int(abs($y)) != int(abs($y) + 1e-14);
923 $re = "$x" if abs($x) >= 1e-14;
924 if ($y == 1) { $im = 'i' }
925 elsif ($y == -1) { $im = '-i' }
926 elsif (abs($y) >= 1e-14) { $im = $y . "i" }
929 $str = $re if defined $re;
930 $str .= "+$im" if defined $im;
933 $str = '0' unless $str;
941 # Stringify as a polar representation '[r,t]'.
943 sub stringify_polar {
945 my ($r, $t) = @{$z->polar};
949 return '[0,0]' if $r <= $eps;
953 $nt = ($nt - int($nt)) * $tpi;
954 $nt += $tpi if $nt < 0; # Range [0, 2pi]
956 if (abs($nt) <= $eps) { $theta = 0 }
957 elsif (abs(pi-$nt) <= $eps) { $theta = 'pi' }
959 if (defined $theta) {
960 $r = int($r + ($r < 0 ? -1 : 1) * $eps)
961 if int(abs($r)) != int(abs($r) + $eps);
962 $theta = int($theta + ($theta < 0 ? -1 : 1) * $eps)
963 if ($theta ne 'pi' and
964 int(abs($theta)) != int(abs($theta) + $eps));
965 return "\[$r,$theta\]";
969 # Okay, number is not a real. Try to identify pi/n and friends...
972 $nt -= $tpi if $nt > pi;
975 for ($k = 1, $kpi = pi; $k < 10; $k++, $kpi += pi) {
976 $n = int($kpi / $nt + ($nt > 0 ? 1 : -1) * 0.5);
977 if (abs($kpi/$n - $nt) <= $eps) {
978 $theta = ($nt < 0 ? '-':'').
979 ($k == 1 ? 'pi':"${k}pi").'/'.abs($n);
984 $theta = $nt unless defined $theta;
986 $r = int($r + ($r < 0 ? -1 : 1) * $eps)
987 if int(abs($r)) != int(abs($r) + $eps);
988 $theta = int($theta + ($theta < 0 ? -1 : 1) * $eps)
989 if ($theta !~ m(^-?\d*pi/\d+$) and
990 int(abs($theta)) != int(abs($theta) + $eps));
992 return "\[$r,$theta\]";
1000 Math::Complex - complex numbers and associated mathematical functions
1005 $z = Math::Complex->make(5, 6);
1007 $j = cplxe(1, 2*pi/3);
1011 This package lets you create and manipulate complex numbers. By default,
1012 I<Perl> limits itself to real numbers, but an extra C<use> statement brings
1013 full complex support, along with a full set of mathematical functions
1014 typically associated with and/or extended to complex numbers.
1016 If you wonder what complex numbers are, they were invented to be able to solve
1017 the following equation:
1021 and by definition, the solution is noted I<i> (engineers use I<j> instead since
1022 I<i> usually denotes an intensity, but the name does not matter). The number
1023 I<i> is a pure I<imaginary> number.
1025 The arithmetics with pure imaginary numbers works just like you would expect
1026 it with real numbers... you just have to remember that
1032 5i + 7i = i * (5 + 7) = 12i
1033 4i - 3i = i * (4 - 3) = i
1038 Complex numbers are numbers that have both a real part and an imaginary
1039 part, and are usually noted:
1043 where C<a> is the I<real> part and C<b> is the I<imaginary> part. The
1044 arithmetic with complex numbers is straightforward. You have to
1045 keep track of the real and the imaginary parts, but otherwise the
1046 rules used for real numbers just apply:
1048 (4 + 3i) + (5 - 2i) = (4 + 5) + i(3 - 2) = 9 + i
1049 (2 + i) * (4 - i) = 2*4 + 4i -2i -i*i = 8 + 2i + 1 = 9 + 2i
1051 A graphical representation of complex numbers is possible in a plane
1052 (also called the I<complex plane>, but it's really a 2D plane).
1057 is the point whose coordinates are (a, b). Actually, it would
1058 be the vector originating from (0, 0) to (a, b). It follows that the addition
1059 of two complex numbers is a vectorial addition.
1061 Since there is a bijection between a point in the 2D plane and a complex
1062 number (i.e. the mapping is unique and reciprocal), a complex number
1063 can also be uniquely identified with polar coordinates:
1067 where C<rho> is the distance to the origin, and C<theta> the angle between
1068 the vector and the I<x> axis. There is a notation for this using the
1069 exponential form, which is:
1071 rho * exp(i * theta)
1073 where I<i> is the famous imaginary number introduced above. Conversion
1074 between this form and the cartesian form C<a + bi> is immediate:
1076 a = rho * cos(theta)
1077 b = rho * sin(theta)
1079 which is also expressed by this formula:
1081 z = rho * exp(i * theta) = rho * (cos theta + i * sin theta)
1083 In other words, it's the projection of the vector onto the I<x> and I<y>
1084 axes. Mathematicians call I<rho> the I<norm> or I<modulus> and I<theta>
1085 the I<argument> of the complex number. The I<norm> of C<z> will be
1088 The polar notation (also known as the trigonometric
1089 representation) is much more handy for performing multiplications and
1090 divisions of complex numbers, whilst the cartesian notation is better
1091 suited for additions and substractions. Real numbers are on the I<x>
1092 axis, and therefore I<theta> is zero.
1094 All the common operations that can be performed on a real number have
1095 been defined to work on complex numbers as well, and are merely
1096 I<extensions> of the operations defined on real numbers. This means
1097 they keep their natural meaning when there is no imaginary part, provided
1098 the number is within their definition set.
1100 For instance, the C<sqrt> routine which computes the square root of
1101 its argument is only defined for positive real numbers and yields a
1102 positive real number (it is an application from B<R+> to B<R+>).
1103 If we allow it to return a complex number, then it can be extended to
1104 negative real numbers to become an application from B<R> to B<C> (the
1105 set of complex numbers):
1107 sqrt(x) = x >= 0 ? sqrt(x) : sqrt(-x)*i
1109 It can also be extended to be an application from B<C> to B<C>,
1110 whilst its restriction to B<R> behaves as defined above by using
1111 the following definition:
1113 sqrt(z = [r,t]) = sqrt(r) * exp(i * t/2)
1115 Indeed, a negative real number can be noted C<[x,pi]>
1116 (the modulus I<x> is always positive, so C<[x,pi]> is really C<-x>, a
1118 and the above definition states that
1120 sqrt([x,pi]) = sqrt(x) * exp(i*pi/2) = [sqrt(x),pi/2] = sqrt(x)*i
1122 which is exactly what we had defined for negative real numbers above.
1124 All the common mathematical functions defined on real numbers that
1125 are extended to complex numbers share that same property of working
1126 I<as usual> when the imaginary part is zero (otherwise, it would not
1127 be called an extension, would it?).
1129 A I<new> operation possible on a complex number that is
1130 the identity for real numbers is called the I<conjugate>, and is noted
1131 with an horizontal bar above the number, or C<~z> here.
1138 z * ~z = (a + bi) * (a - bi) = a*a + b*b
1140 We saw that the norm of C<z> was noted C<abs(z)> and was defined as the
1141 distance to the origin, also known as:
1143 rho = abs(z) = sqrt(a*a + b*b)
1147 z * ~z = abs(z) ** 2
1149 If z is a pure real number (i.e. C<b == 0>), then the above yields:
1153 which is true (C<abs> has the regular meaning for real number, i.e. stands
1154 for the absolute value). This example explains why the norm of C<z> is
1155 noted C<abs(z)>: it extends the C<abs> function to complex numbers, yet
1156 is the regular C<abs> we know when the complex number actually has no
1157 imaginary part... This justifies I<a posteriori> our use of the C<abs>
1158 notation for the norm.
1162 Given the following notations:
1164 z1 = a + bi = r1 * exp(i * t1)
1165 z2 = c + di = r2 * exp(i * t2)
1166 z = <any complex or real number>
1168 the following (overloaded) operations are supported on complex numbers:
1170 z1 + z2 = (a + c) + i(b + d)
1171 z1 - z2 = (a - c) + i(b - d)
1172 z1 * z2 = (r1 * r2) * exp(i * (t1 + t2))
1173 z1 / z2 = (r1 / r2) * exp(i * (t1 - t2))
1174 z1 ** z2 = exp(z2 * log z1)
1176 abs(z1) = r1 = sqrt(a*a + b*b)
1177 sqrt(z1) = sqrt(r1) * exp(i * t1/2)
1178 exp(z1) = exp(a) * exp(i * b)
1179 log(z1) = log(r1) + i*t1
1180 sin(z1) = 1/2i (exp(i * z1) - exp(-i * z1))
1181 cos(z1) = 1/2 (exp(i * z1) + exp(-i * z1))
1183 atan2(z1, z2) = atan(z1/z2)
1185 The following extra operations are supported on both real and complex
1192 cbrt(z) = z ** (1/3)
1193 log10(z) = log(z) / log(10)
1194 logn(z, n) = log(z) / log(n)
1196 tan(z) = sin(z) / cos(z)
1202 asin(z) = -i * log(i*z + sqrt(1-z*z))
1203 acos(z) = -i * log(z + sqrt(z*z-1))
1204 atan(z) = i/2 * log((i+z) / (i-z))
1206 acsc(z) = asin(1 / z)
1207 asec(z) = acos(1 / z)
1208 acot(z) = -i/2 * log((i+z) / (z-i))
1210 sinh(z) = 1/2 (exp(z) - exp(-z))
1211 cosh(z) = 1/2 (exp(z) + exp(-z))
1212 tanh(z) = sinh(z) / cosh(z) = (exp(z) - exp(-z)) / (exp(z) + exp(-z))
1214 csch(z) = 1 / sinh(z)
1215 sech(z) = 1 / cosh(z)
1216 coth(z) = 1 / tanh(z)
1218 asinh(z) = log(z + sqrt(z*z+1))
1219 acosh(z) = log(z + sqrt(z*z-1))
1220 atanh(z) = 1/2 * log((1+z) / (1-z))
1222 acsch(z) = asinh(1 / z)
1223 asech(z) = acosh(1 / z)
1224 acoth(z) = atanh(1 / z) = 1/2 * log((1+z) / (z-1))
1226 I<log>, I<csc>, I<cot>, I<acsc>, I<acot>, I<csch>, I<coth>,
1227 I<acosech>, I<acotanh>, have aliases I<ln>, I<cosec>, I<cotan>,
1228 I<acosec>, I<acotan>, I<cosech>, I<cotanh>, I<acosech>, I<acotanh>,
1231 The I<root> function is available to compute all the I<n>
1232 roots of some complex, where I<n> is a strictly positive integer.
1233 There are exactly I<n> such roots, returned as a list. Getting the
1234 number mathematicians call C<j> such that:
1238 is a simple matter of writing:
1240 $j = ((root(1, 3))[1];
1242 The I<k>th root for C<z = [r,t]> is given by:
1244 (root(z, n))[k] = r**(1/n) * exp(i * (t + 2*k*pi)/n)
1246 The I<spaceship> comparison operator is also defined. In order to
1247 ensure its restriction to real numbers is conform to what you would
1248 expect, the comparison is run on the real part of the complex number
1249 first, and imaginary parts are compared only when the real parts
1254 To create a complex number, use either:
1256 $z = Math::Complex->make(3, 4);
1259 if you know the cartesian form of the number, or
1263 if you like. To create a number using the trigonometric form, use either:
1265 $z = Math::Complex->emake(5, pi/3);
1266 $x = cplxe(5, pi/3);
1268 instead. The first argument is the modulus, the second is the angle
1269 (in radians, the full circle is 2*pi). (Mnmemonic: C<e> is used as a
1270 notation for complex numbers in the trigonometric form).
1272 It is possible to write:
1274 $x = cplxe(-3, pi/4);
1276 but that will be silently converted into C<[3,-3pi/4]>, since the modulus
1277 must be positive (it represents the distance to the origin in the complex
1280 =head1 STRINGIFICATION
1282 When printed, a complex number is usually shown under its cartesian
1283 form I<a+bi>, but there are legitimate cases where the polar format
1284 I<[r,t]> is more appropriate.
1286 By calling the routine C<Math::Complex::display_format> and supplying either
1287 C<"polar"> or C<"cartesian">, you override the default display format,
1288 which is C<"cartesian">. Not supplying any argument returns the current
1291 This default can be overridden on a per-number basis by calling the
1292 C<display_format> method instead. As before, not supplying any argument
1293 returns the current display format for this number. Otherwise whatever you
1294 specify will be the new display format for I<this> particular number.
1300 Math::Complex::display_format('polar');
1301 $j = ((root(1, 3))[1];
1302 print "j = $j\n"; # Prints "j = [1,2pi/3]
1303 $j->display_format('cartesian');
1304 print "j = $j\n"; # Prints "j = -0.5+0.866025403784439i"
1306 The polar format attempts to emphasize arguments like I<k*pi/n>
1307 (where I<n> is a positive integer and I<k> an integer within [-9,+9]).
1311 Thanks to overloading, the handling of arithmetics with complex numbers
1312 is simple and almost transparent.
1314 Here are some examples:
1318 $j = cplxe(1, 2*pi/3); # $j ** 3 == 1
1319 print "j = $j, j**3 = ", $j ** 3, "\n";
1320 print "1 + j + j**2 = ", 1 + $j + $j**2, "\n";
1322 $z = -16 + 0*i; # Force it to be a complex
1323 print "sqrt($z) = ", sqrt($z), "\n";
1325 $k = exp(i * 2*pi/3);
1326 print "$j - $k = ", $j - $k, "\n";
1330 Saying C<use Math::Complex;> exports many mathematical routines in the caller
1331 environment. This is construed as a feature by the Author, actually... ;-)
1333 The code is not optimized for speed, although we try to use the cartesian
1334 form for addition-like operators and the trigonometric form for all
1335 multiplication-like operators.
1337 The arg() routine does not ensure the angle is within the range [-pi,+pi]
1338 (a side effect caused by multiplication and division using the trigonometric
1341 All routines expect to be given real or complex numbers. Don't attempt to
1342 use BigFloat, since Perl has currently no rule to disambiguate a '+'
1343 operation (for instance) between two overloaded entities.
1347 Raphael Manfredi <F<Raphael_Manfredi@grenoble.hp.com>>
1348 Jarkko Hietaniemi <F<jhi@iki.fi>>