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);
516 # Die on division by zero.
519 my $mess = "$_[0]: Logarithm of zero.\n";
522 $mess .= "(Because in the definition of $_[0], the argument ";
523 $mess .= "$_[1] " unless ($_[1] eq '0');
529 $mess .= "Died at $up[1] line $up[2].\n";
541 $z = cplx($z, 0) unless ref $z;
542 my ($x, $y) = @{$z->cartesian};
543 my ($r, $t) = @{$z->polar};
544 $t -= 2 * pi if ($t > pi() and $x < 0);
545 $t += 2 * pi if ($t < -pi() and $x < 0);
546 return (ref $z)->make(log($r), $t);
554 sub ln { Math::Complex::log(@_) }
565 return log(cplx($z, 0)) * log10inv unless ref $z;
566 my ($r, $t) = @{$z->polar};
567 return (ref $z)->make(log($r) * log10inv, $t * log10inv);
573 # Compute logn(z,n) = log(z) / log(n)
577 $z = cplx($z, 0) unless ref $z;
578 my $logn = $logn{$n};
579 $logn = $logn{$n} = log($n) unless defined $logn; # Cache log(n)
580 return log($z) / $logn;
586 # Compute cos(z) = (exp(iz) + exp(-iz))/2.
590 $z = cplx($z, 0) unless ref $z;
591 my ($x, $y) = @{$z->cartesian};
594 return (ref $z)->make(cos($x) * ($ey + $ey_1)/2,
595 sin($x) * ($ey_1 - $ey)/2);
601 # Compute sin(z) = (exp(iz) - exp(-iz))/2.
605 $z = cplx($z, 0) unless ref $z;
606 my ($x, $y) = @{$z->cartesian};
609 return (ref $z)->make(sin($x) * ($ey + $ey_1)/2,
610 cos($x) * ($ey - $ey_1)/2);
616 # Compute tan(z) = sin(z) / cos(z).
621 _divbyzero "tan($z)", "cos($z)" if ($cz == 0);
622 return sin($z) / $cz;
628 # Computes the secant sec(z) = 1 / cos(z).
633 _divbyzero "sec($z)", "cos($z)" if ($cz == 0);
640 # Computes the cosecant csc(z) = 1 / sin(z).
645 _divbyzero "csc($z)", "sin($z)" if ($sz == 0);
654 sub cosec { Math::Complex::csc(@_) }
659 # Computes cot(z) = 1 / tan(z).
664 _divbyzero "cot($z)", "sin($z)" if ($sz == 0);
665 return cos($z) / $sz;
673 sub cotan { Math::Complex::cot(@_) }
678 # Computes the arc cosine acos(z) = -i log(z + sqrt(z*z-1)).
682 $z = cplx($z, 0) unless ref $z;
683 my ($re, $im) = @{$z->cartesian};
684 return atan2(sqrt(1 - $re * $re), $re)
685 if ($im == 0 and abs($re) <= 1.0);
686 my $acos = ~i * log($z + sqrt($z*$z - 1));
688 (abs($re) < 1 && abs($im) < 1) ||
689 (abs($re) > 1 && abs($im) > 1
690 && !($re > 1 && $im > 1)
691 && !($re < -1 && $im < -1))) {
692 # this rule really, REALLY, must be simpler
701 # Computes the arc sine asin(z) = -i log(iz + sqrt(1-z*z)).
705 $z = cplx($z, 0) unless ref $z;
706 my ($re, $im) = @{$z->cartesian};
707 return atan2($re, sqrt(1 - $re * $re))
708 if ($im == 0 and abs($re) <= 1.0);
709 return ~i * log(i * $z + sqrt(1 - $z*$z));
715 # Computes the arc tangent atan(z) = i/2 log((i+z) / (i-z)).
719 $z = cplx($z, 0) unless ref $z;
720 _divbyzero "atan(i)" if ( $z == i);
721 _divbyzero "atan(-i)" if (-$z == i);
722 return i/2*log((i + $z) / (i - $z));
728 # Computes the arc secant asec(z) = acos(1 / z).
732 _divbyzero "asec($z)", $z if ($z == 0);
733 $z = cplx($z, 0) unless ref $z;
734 my ($re, $im) = @{$z->cartesian};
735 if ($im == 0 && abs($re) >= 1.0) {
737 return atan2(sqrt(1 - $ire * $ire), $ire);
739 my $asec = acos(1 / $z);
740 return ~$asec if $re < 0 && $re > -1 && $im == 0;
741 return -$asec if $im && !($re > 0 && $im > 0) && !($re < 0 && $im < 0);
748 # Computes the arc cosecant acsc(z) = asin(1 / z).
752 _divbyzero "acsc($z)", $z if ($z == 0);
753 $z = cplx($z, 0) unless ref $z;
754 my ($re, $im) = @{$z->cartesian};
755 if ($im == 0 && abs($re) >= 1.0) {
757 return atan2($ire, sqrt(1 - $ire * $ire));
759 my $acsc = asin(1 / $z);
760 return ~$acsc if $re < 0 && $re > -1 && $im == 0;
769 sub acosec { Math::Complex::acsc(@_) }
774 # Computes the arc cotangent acot(z) = atan(1 / z)
778 _divbyzero "acot($z)" if ($z == 0);
779 $z = cplx($z, 0) unless ref $z;
780 _divbyzero "acot(i)", if ( $z == i);
781 _divbyzero "acot(-i)" if (-$z == i);
790 sub acotan { Math::Complex::acot(@_) }
795 # Computes the hyperbolic cosine cosh(z) = (exp(z) + exp(-z))/2.
804 my ($x, $y) = @{$z->cartesian};
807 return cplx(0.5 * ($ex + $ex_1), 0) if $real;
808 return (ref $z)->make(cos($y) * ($ex + $ex_1)/2,
809 sin($y) * ($ex - $ex_1)/2);
815 # Computes the hyperbolic sine sinh(z) = (exp(z) - exp(-z))/2.
824 my ($x, $y) = @{$z->cartesian};
827 return cplx(0.5 * ($ex - $ex_1), 0) if $real;
828 return (ref $z)->make(cos($y) * ($ex - $ex_1)/2,
829 sin($y) * ($ex + $ex_1)/2);
835 # Computes the hyperbolic tangent tanh(z) = sinh(z) / cosh(z).
840 _divbyzero "tanh($z)", "cosh($z)" if ($cz == 0);
841 return sinh($z) / $cz;
847 # Computes the hyperbolic secant sech(z) = 1 / cosh(z).
852 _divbyzero "sech($z)", "cosh($z)" if ($cz == 0);
859 # Computes the hyperbolic cosecant csch(z) = 1 / sinh(z).
864 _divbyzero "csch($z)", "sinh($z)" if ($sz == 0);
873 sub cosech { Math::Complex::csch(@_) }
878 # Computes the hyperbolic cotangent coth(z) = cosh(z) / sinh(z).
883 _divbyzero "coth($z)", "sinh($z)" if ($sz == 0);
884 return cosh($z) / $sz;
892 sub cotanh { Math::Complex::coth(@_) }
897 # Computes the arc hyperbolic cosine acosh(z) = log(z +- sqrt(z*z-1)).
901 $z = cplx($z, 0) unless ref $z;
902 my ($re, $im) = @{$z->cartesian};
903 return log($re + sqrt(cplx($re*$re - 1, 0)))
904 if ($im == 0 && $re < 0);
905 return log($z + sqrt($z*$z - 1));
911 # Computes the arc hyperbolic sine asinh(z) = log(z + sqrt(z*z-1))
915 $z = cplx($z, 0) unless ref $z;
916 return log($z + sqrt($z*$z + 1));
922 # Computes the arc hyperbolic tangent atanh(z) = 1/2 log((1+z) / (1-z)).
926 _divbyzero 'atanh(1)', "1 - $z" if ($z == 1);
927 _logofzero 'atanh(-1)' if ($z == -1);
928 $z = cplx($z, 0) unless ref $z;
929 my ($re, $im) = @{$z->cartesian};
930 if ($im == 0 && $re > 1) {
931 return cplx(atanh(1 / $re), pi/2);
933 return log((1 + $z) / (1 - $z)) / 2;
939 # Computes the hyperbolic arc secant asech(z) = acosh(1 / z).
943 _divbyzero 'asech(0)', $z if ($z == 0);
944 $z = cplx($z, 0) unless ref $z;
945 my ($re, $im) = @{$z->cartesian};
946 if ($im == 0 && $re < 0) {
948 return log($ire + sqrt(cplx($ire*$ire - 1, 0)));
950 return acosh(1 / $z);
956 # Computes the hyperbolic arc cosecant acsch(z) = asinh(1 / z).
960 _divbyzero 'acsch(0)', $z if ($z == 0);
961 return asinh(1 / $z);
969 sub acosech { Math::Complex::acsch(@_) }
974 # Computes the arc hyperbolic cotangent acoth(z) = 1/2 log((1+z) / (z-1)).
978 _divbyzero 'acoth(1)', "$z - 1" if ($z == 1);
979 _logofzero 'acoth(-1)' if ($z == -1);
980 $z = cplx($z, 0) unless ref $z;
981 my ($re, $im) = @{$z->cartesian};
982 if ($im == 0 and abs($re) < 1) {
983 return cplx(acoth(1/$re) , pi/2);
985 return log((1 + $z) / ($z - 1)) / 2;
993 sub acotanh { Math::Complex::acoth(@_) }
998 # Compute atan(z1/z2).
1001 my ($z1, $z2, $inverted) = @_;
1002 my ($re1, $im1) = ref $z1 ? @{$z1->cartesian} : ($z1, 0);
1003 my ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
1005 if (defined $inverted && $inverted) { # atan(z2/z1)
1006 return pi * ($re2 > 0 ? 1 : -1) if $re1 == 0 && $im1 == 0;
1009 return pi * ($re1 > 0 ? 1 : -1) if $re2 == 0 && $im2 == 0;
1019 # Set (fetch if no argument) display format for all complex numbers that
1020 # don't happen to have overrriden it via ->display_format
1022 # When called as a method, this actually sets the display format for
1023 # the current object.
1025 # Valid object formats are 'c' and 'p' for cartesian and polar. The first
1026 # letter is used actually, so the type can be fully spelled out for clarity.
1028 sub display_format {
1032 if (ref $self) { # Called as a method
1034 } else { # Regular procedure call
1039 if (defined $self) {
1040 return defined $self->{display} ? $self->{display} : $display
1041 unless defined $format;
1042 return $self->{display} = $format;
1045 return $display unless defined $format;
1046 return $display = $format;
1052 # Show nicely formatted complex number under its cartesian or polar form,
1053 # depending on the current display format:
1055 # . If a specific display format has been recorded for this object, use it.
1056 # . Otherwise, use the generic current default for all complex numbers,
1057 # which is a package global variable.
1064 $format = $z->{display} if defined $z->{display};
1066 return $z->stringify_polar if $format =~ /^p/i;
1067 return $z->stringify_cartesian;
1071 # ->stringify_cartesian
1073 # Stringify as a cartesian representation 'a+bi'.
1075 sub stringify_cartesian {
1077 my ($x, $y) = @{$z->cartesian};
1080 $x = int($x + ($x < 0 ? -1 : 1) * 1e-14)
1081 if int(abs($x)) != int(abs($x) + 1e-14);
1082 $y = int($y + ($y < 0 ? -1 : 1) * 1e-14)
1083 if int(abs($y)) != int(abs($y) + 1e-14);
1085 $re = "$x" if abs($x) >= 1e-14;
1086 if ($y == 1) { $im = 'i' }
1087 elsif ($y == -1) { $im = '-i' }
1088 elsif (abs($y) >= 1e-14) { $im = $y . "i" }
1091 $str = $re if defined $re;
1092 $str .= "+$im" if defined $im;
1095 $str = '0' unless $str;
1103 # Stringify as a polar representation '[r,t]'.
1105 sub stringify_polar {
1107 my ($r, $t) = @{$z->polar};
1111 return '[0,0]' if $r <= $eps;
1115 $nt = ($nt - int($nt)) * $tpi;
1116 $nt += $tpi if $nt < 0; # Range [0, 2pi]
1118 if (abs($nt) <= $eps) { $theta = 0 }
1119 elsif (abs(pi-$nt) <= $eps) { $theta = 'pi' }
1121 if (defined $theta) {
1122 $r = int($r + ($r < 0 ? -1 : 1) * $eps)
1123 if int(abs($r)) != int(abs($r) + $eps);
1124 $theta = int($theta + ($theta < 0 ? -1 : 1) * $eps)
1125 if ($theta ne 'pi' and
1126 int(abs($theta)) != int(abs($theta) + $eps));
1127 return "\[$r,$theta\]";
1131 # Okay, number is not a real. Try to identify pi/n and friends...
1134 $nt -= $tpi if $nt > pi;
1137 for ($k = 1, $kpi = pi; $k < 10; $k++, $kpi += pi) {
1138 $n = int($kpi / $nt + ($nt > 0 ? 1 : -1) * 0.5);
1139 if (abs($kpi/$n - $nt) <= $eps) {
1140 $theta = ($nt < 0 ? '-':'').
1141 ($k == 1 ? 'pi':"${k}pi").'/'.abs($n);
1146 $theta = $nt unless defined $theta;
1148 $r = int($r + ($r < 0 ? -1 : 1) * $eps)
1149 if int(abs($r)) != int(abs($r) + $eps);
1150 $theta = int($theta + ($theta < 0 ? -1 : 1) * $eps)
1151 if ($theta !~ m(^-?\d*pi/\d+$) and
1152 int(abs($theta)) != int(abs($theta) + $eps));
1154 return "\[$r,$theta\]";
1162 Math::Complex - complex numbers and associated mathematical functions
1168 $z = Math::Complex->make(5, 6);
1170 $j = cplxe(1, 2*pi/3);
1174 This package lets you create and manipulate complex numbers. By default,
1175 I<Perl> limits itself to real numbers, but an extra C<use> statement brings
1176 full complex support, along with a full set of mathematical functions
1177 typically associated with and/or extended to complex numbers.
1179 If you wonder what complex numbers are, they were invented to be able to solve
1180 the following equation:
1184 and by definition, the solution is noted I<i> (engineers use I<j> instead since
1185 I<i> usually denotes an intensity, but the name does not matter). The number
1186 I<i> is a pure I<imaginary> number.
1188 The arithmetics with pure imaginary numbers works just like you would expect
1189 it with real numbers... you just have to remember that
1195 5i + 7i = i * (5 + 7) = 12i
1196 4i - 3i = i * (4 - 3) = i
1201 Complex numbers are numbers that have both a real part and an imaginary
1202 part, and are usually noted:
1206 where C<a> is the I<real> part and C<b> is the I<imaginary> part. The
1207 arithmetic with complex numbers is straightforward. You have to
1208 keep track of the real and the imaginary parts, but otherwise the
1209 rules used for real numbers just apply:
1211 (4 + 3i) + (5 - 2i) = (4 + 5) + i(3 - 2) = 9 + i
1212 (2 + i) * (4 - i) = 2*4 + 4i -2i -i*i = 8 + 2i + 1 = 9 + 2i
1214 A graphical representation of complex numbers is possible in a plane
1215 (also called the I<complex plane>, but it's really a 2D plane).
1220 is the point whose coordinates are (a, b). Actually, it would
1221 be the vector originating from (0, 0) to (a, b). It follows that the addition
1222 of two complex numbers is a vectorial addition.
1224 Since there is a bijection between a point in the 2D plane and a complex
1225 number (i.e. the mapping is unique and reciprocal), a complex number
1226 can also be uniquely identified with polar coordinates:
1230 where C<rho> is the distance to the origin, and C<theta> the angle between
1231 the vector and the I<x> axis. There is a notation for this using the
1232 exponential form, which is:
1234 rho * exp(i * theta)
1236 where I<i> is the famous imaginary number introduced above. Conversion
1237 between this form and the cartesian form C<a + bi> is immediate:
1239 a = rho * cos(theta)
1240 b = rho * sin(theta)
1242 which is also expressed by this formula:
1244 z = rho * exp(i * theta) = rho * (cos theta + i * sin theta)
1246 In other words, it's the projection of the vector onto the I<x> and I<y>
1247 axes. Mathematicians call I<rho> the I<norm> or I<modulus> and I<theta>
1248 the I<argument> of the complex number. The I<norm> of C<z> will be
1251 The polar notation (also known as the trigonometric
1252 representation) is much more handy for performing multiplications and
1253 divisions of complex numbers, whilst the cartesian notation is better
1254 suited for additions and substractions. Real numbers are on the I<x>
1255 axis, and therefore I<theta> is zero.
1257 All the common operations that can be performed on a real number have
1258 been defined to work on complex numbers as well, and are merely
1259 I<extensions> of the operations defined on real numbers. This means
1260 they keep their natural meaning when there is no imaginary part, provided
1261 the number is within their definition set.
1263 For instance, the C<sqrt> routine which computes the square root of
1264 its argument is only defined for positive real numbers and yields a
1265 positive real number (it is an application from B<R+> to B<R+>).
1266 If we allow it to return a complex number, then it can be extended to
1267 negative real numbers to become an application from B<R> to B<C> (the
1268 set of complex numbers):
1270 sqrt(x) = x >= 0 ? sqrt(x) : sqrt(-x)*i
1272 It can also be extended to be an application from B<C> to B<C>,
1273 whilst its restriction to B<R> behaves as defined above by using
1274 the following definition:
1276 sqrt(z = [r,t]) = sqrt(r) * exp(i * t/2)
1278 Indeed, a negative real number can be noted C<[x,pi]>
1279 (the modulus I<x> is always positive, so C<[x,pi]> is really C<-x>, a
1281 and the above definition states that
1283 sqrt([x,pi]) = sqrt(x) * exp(i*pi/2) = [sqrt(x),pi/2] = sqrt(x)*i
1285 which is exactly what we had defined for negative real numbers above.
1287 All the common mathematical functions defined on real numbers that
1288 are extended to complex numbers share that same property of working
1289 I<as usual> when the imaginary part is zero (otherwise, it would not
1290 be called an extension, would it?).
1292 A I<new> operation possible on a complex number that is
1293 the identity for real numbers is called the I<conjugate>, and is noted
1294 with an horizontal bar above the number, or C<~z> here.
1301 z * ~z = (a + bi) * (a - bi) = a*a + b*b
1303 We saw that the norm of C<z> was noted C<abs(z)> and was defined as the
1304 distance to the origin, also known as:
1306 rho = abs(z) = sqrt(a*a + b*b)
1310 z * ~z = abs(z) ** 2
1312 If z is a pure real number (i.e. C<b == 0>), then the above yields:
1316 which is true (C<abs> has the regular meaning for real number, i.e. stands
1317 for the absolute value). This example explains why the norm of C<z> is
1318 noted C<abs(z)>: it extends the C<abs> function to complex numbers, yet
1319 is the regular C<abs> we know when the complex number actually has no
1320 imaginary part... This justifies I<a posteriori> our use of the C<abs>
1321 notation for the norm.
1325 Given the following notations:
1327 z1 = a + bi = r1 * exp(i * t1)
1328 z2 = c + di = r2 * exp(i * t2)
1329 z = <any complex or real number>
1331 the following (overloaded) operations are supported on complex numbers:
1333 z1 + z2 = (a + c) + i(b + d)
1334 z1 - z2 = (a - c) + i(b - d)
1335 z1 * z2 = (r1 * r2) * exp(i * (t1 + t2))
1336 z1 / z2 = (r1 / r2) * exp(i * (t1 - t2))
1337 z1 ** z2 = exp(z2 * log z1)
1339 abs(z1) = r1 = sqrt(a*a + b*b)
1340 sqrt(z1) = sqrt(r1) * exp(i * t1/2)
1341 exp(z1) = exp(a) * exp(i * b)
1342 log(z1) = log(r1) + i*t1
1343 sin(z1) = 1/2i (exp(i * z1) - exp(-i * z1))
1344 cos(z1) = 1/2 (exp(i * z1) + exp(-i * z1))
1346 atan2(z1, z2) = atan(z1/z2)
1348 The following extra operations are supported on both real and complex
1355 cbrt(z) = z ** (1/3)
1356 log10(z) = log(z) / log(10)
1357 logn(z, n) = log(z) / log(n)
1359 tan(z) = sin(z) / cos(z)
1365 asin(z) = -i * log(i*z + sqrt(1-z*z))
1366 acos(z) = -i * log(z + sqrt(z*z-1))
1367 atan(z) = i/2 * log((i+z) / (i-z))
1369 acsc(z) = asin(1 / z)
1370 asec(z) = acos(1 / z)
1371 acot(z) = atan(1 / z) = -i/2 * log((i+z) / (z-i))
1373 sinh(z) = 1/2 (exp(z) - exp(-z))
1374 cosh(z) = 1/2 (exp(z) + exp(-z))
1375 tanh(z) = sinh(z) / cosh(z) = (exp(z) - exp(-z)) / (exp(z) + exp(-z))
1377 csch(z) = 1 / sinh(z)
1378 sech(z) = 1 / cosh(z)
1379 coth(z) = 1 / tanh(z)
1381 asinh(z) = log(z + sqrt(z*z+1))
1382 acosh(z) = log(z + sqrt(z*z-1))
1383 atanh(z) = 1/2 * log((1+z) / (1-z))
1385 acsch(z) = asinh(1 / z)
1386 asech(z) = acosh(1 / z)
1387 acoth(z) = atanh(1 / z) = 1/2 * log((1+z) / (z-1))
1389 I<log>, I<csc>, I<cot>, I<acsc>, I<acot>, I<csch>, I<coth>,
1390 I<acosech>, I<acotanh>, have aliases I<ln>, I<cosec>, I<cotan>,
1391 I<acosec>, I<acotan>, I<cosech>, I<cotanh>, I<acosech>, I<acotanh>,
1394 The I<root> function is available to compute all the I<n>
1395 roots of some complex, where I<n> is a strictly positive integer.
1396 There are exactly I<n> such roots, returned as a list. Getting the
1397 number mathematicians call C<j> such that:
1401 is a simple matter of writing:
1403 $j = ((root(1, 3))[1];
1405 The I<k>th root for C<z = [r,t]> is given by:
1407 (root(z, n))[k] = r**(1/n) * exp(i * (t + 2*k*pi)/n)
1409 The I<spaceship> comparison operator, E<lt>=E<gt>, is also defined. In
1410 order to ensure its restriction to real numbers is conform to what you
1411 would expect, the comparison is run on the real part of the complex
1412 number first, and imaginary parts are compared only when the real
1417 To create a complex number, use either:
1419 $z = Math::Complex->make(3, 4);
1422 if you know the cartesian form of the number, or
1426 if you like. To create a number using the trigonometric form, use either:
1428 $z = Math::Complex->emake(5, pi/3);
1429 $x = cplxe(5, pi/3);
1431 instead. The first argument is the modulus, the second is the angle
1432 (in radians, the full circle is 2*pi). (Mnmemonic: C<e> is used as a
1433 notation for complex numbers in the trigonometric form).
1435 It is possible to write:
1437 $x = cplxe(-3, pi/4);
1439 but that will be silently converted into C<[3,-3pi/4]>, since the modulus
1440 must be positive (it represents the distance to the origin in the complex
1443 =head1 STRINGIFICATION
1445 When printed, a complex number is usually shown under its cartesian
1446 form I<a+bi>, but there are legitimate cases where the polar format
1447 I<[r,t]> is more appropriate.
1449 By calling the routine C<Math::Complex::display_format> and supplying either
1450 C<"polar"> or C<"cartesian">, you override the default display format,
1451 which is C<"cartesian">. Not supplying any argument returns the current
1454 This default can be overridden on a per-number basis by calling the
1455 C<display_format> method instead. As before, not supplying any argument
1456 returns the current display format for this number. Otherwise whatever you
1457 specify will be the new display format for I<this> particular number.
1463 Math::Complex::display_format('polar');
1464 $j = ((root(1, 3))[1];
1465 print "j = $j\n"; # Prints "j = [1,2pi/3]
1466 $j->display_format('cartesian');
1467 print "j = $j\n"; # Prints "j = -0.5+0.866025403784439i"
1469 The polar format attempts to emphasize arguments like I<k*pi/n>
1470 (where I<n> is a positive integer and I<k> an integer within [-9,+9]).
1474 Thanks to overloading, the handling of arithmetics with complex numbers
1475 is simple and almost transparent.
1477 Here are some examples:
1481 $j = cplxe(1, 2*pi/3); # $j ** 3 == 1
1482 print "j = $j, j**3 = ", $j ** 3, "\n";
1483 print "1 + j + j**2 = ", 1 + $j + $j**2, "\n";
1485 $z = -16 + 0*i; # Force it to be a complex
1486 print "sqrt($z) = ", sqrt($z), "\n";
1488 $k = exp(i * 2*pi/3);
1489 print "$j - $k = ", $j - $k, "\n";
1491 =head1 ERRORS DUE TO DIVISION BY ZERO
1493 The division (/) and the following functions
1512 cannot be computed for all arguments because that would mean dividing
1513 by zero or taking logarithm of zero. These situations cause fatal
1514 runtime errors looking like this
1516 cot(0): Division by zero.
1517 (Because in the definition of cot(0), the divisor sin(0) is 0)
1522 atanh(-1): Logarithm of zero.
1525 For the C<csc>, C<cot>, C<asec>, C<acsc>, C<acot>, C<csch>, C<coth>,
1526 C<asech>, C<acsch>, the argument cannot be C<0> (zero). For the
1527 C<atanh>, C<acoth>, the argument cannot be C<1> (one). For the
1528 C<atanh>, C<acoth>, the argument cannot be C<-1> (minus one). For the
1529 C<atan>, C<acot>, the argument cannot be C<i> (the imaginary unit).
1530 For the C<atan>, C<acoth>, the argument cannot be C<-i> (the negative
1531 imaginary unit). For the C<tan>, C<sec>, C<tanh>, C<sech>, the
1532 argument cannot be I<pi/2 + k * pi>, where I<k> is any integer.
1536 Saying C<use Math::Complex;> exports many mathematical routines in the
1537 caller environment and even overrides some (C<sin>, C<cos>, C<sqrt>,
1538 C<log>, C<exp>). This is construed as a feature by the Authors,
1541 The code is not optimized for speed, although we try to use the cartesian
1542 form for addition-like operators and the trigonometric form for all
1543 multiplication-like operators.
1545 The arg() routine does not ensure the angle is within the range [-pi,+pi]
1546 (a side effect caused by multiplication and division using the trigonometric
1549 All routines expect to be given real or complex numbers. Don't attempt to
1550 use BigFloat, since Perl has currently no rule to disambiguate a '+'
1551 operation (for instance) between two overloaded entities.
1555 Raphael Manfredi <F<Raphael_Manfredi@grenoble.hp.com>> and
1556 Jarkko Hietaniemi <F<jhi@iki.fi>>.