This is a live mirror of the Perl 5 development currently hosted at https://github.com/perl/perl5
Math-Complex: customized VERSION bump
[perl5.git] / cpan / Math-Complex / lib / Math / Complex.pm
CommitLineData
66730be0
RM
1#
2# Complex numbers and associated mathematical functions
b42d0ec9
JH
3# -- Raphael Manfredi Since Sep 1996
4# -- Jarkko Hietaniemi Since Mar 1997
5# -- Daniel S. Lewart Since Sep 1997
fb73857a 6#
a0d0e21e 7
5aabfad6 8package Math::Complex;
a0d0e21e 9
9853179e 10{ use 5.006; }
e7deb4af
SH
11use strict;
12
2c9daa73 13our $VERSION = 1.59_01;
8b4fe368
SP
14
15use Config;
476757f7 16
0720d38b
JH
17our ($Inf, $ExpInf);
18our ($vax_float, $has_inf, $has_nan);
19
9fbe1b12 20BEGIN {
0720d38b
JH
21 $vax_float = (pack("d",1) =~ /^[\x80\x10]\x40/);
22 $has_inf = !$vax_float;
23 $has_nan = !$vax_float;
24
25 unless ($has_inf) {
26 # For example in vax, there is no Inf,
27 # and just mentioning the DBL_MAX (1.70141183460469229e+38)
28 # causes SIGFPE.
29
30 # These are pretty useless without a real infinity,
31 # but setting them makes for less warnings about their
32 # undefined values.
33 $Inf = "Inf";
34 $ExpInf = "Inf";
35 return;
36 }
37
38 my %DBL_MAX = # These are IEEE 754 maxima.
8b4fe368
SP
39 (
40 4 => '1.70141183460469229e+38',
41 8 => '1.7976931348623157e+308',
57dd0abb
SP
42 # AFAICT the 10, 12, and 16-byte long doubles
43 # all have the same maximum.
8b4fe368 44 10 => '1.1897314953572317650857593266280070162E+4932',
57dd0abb 45 12 => '1.1897314953572317650857593266280070162E+4932',
d3056722 46 16 => '1.1897314953572317650857593266280070162E+4932',
8b4fe368 47 );
0720d38b 48
57dd0abb
SP
49 my $nvsize = $Config{nvsize} ||
50 ($Config{uselongdouble} && $Config{longdblsize}) ||
51 $Config{doublesize};
52 die "Math::Complex: Could not figure out nvsize\n"
53 unless defined $nvsize;
54 die "Math::Complex: Cannot not figure out max nv (nvsize = $nvsize)\n"
55 unless defined $DBL_MAX{$nvsize};
8b4fe368 56 my $DBL_MAX = eval $DBL_MAX{$nvsize};
57dd0abb
SP
57 die "Math::Complex: Could not figure out max nv (nvsize = $nvsize)\n"
58 unless defined $DBL_MAX;
8b4fe368 59 my $BIGGER_THAN_THIS = 1e30; # Must find something bigger than this.
1515bec6 60 if ($^O eq 'unicosmk') {
8b4fe368 61 $Inf = $DBL_MAX;
1515bec6 62 } else {
0720d38b 63 local $SIG{FPE} = sub { };
618e05e9 64 local $!;
8b4fe368 65 # We do want an arithmetic overflow, Inf INF inf Infinity.
1515bec6 66 for my $t (
8b4fe368 67 'exp(99999)', # Enough even with 128-bit long doubles.
1515bec6
SP
68 'inf',
69 'Inf',
70 'INF',
71 'infinity',
72 'Infinity',
73 'INFINITY',
8b4fe368 74 '1e99999',
1515bec6 75 ) {
1515bec6
SP
76 local $^W = 0;
77 my $i = eval "$t+1.0";
9f4452f7 78 if (defined $i && $i > $BIGGER_THAN_THIS) {
1515bec6
SP
79 $Inf = $i;
80 last;
81 }
0720d38b 82 }
8b4fe368
SP
83 $Inf = $DBL_MAX unless defined $Inf; # Oh well, close enough.
84 die "Math::Complex: Could not get Infinity"
85 unless $Inf > $BIGGER_THAN_THIS;
0720d38b
JH
86 $ExpInf = eval 'exp(99999)';
87 }
9f4452f7 88 # print "# On this machine, Inf = '$Inf'\n";
9fbe1b12 89}
fb73857a 90
e7deb4af
SH
91use Scalar::Util qw(set_prototype);
92
93use warnings;
94no warnings 'syntax'; # To avoid the (_) warnings.
95
96BEGIN {
97 # For certain functions that we override, in 5.10 or better
98 # we can set a smarter prototype that will handle the lexical $_
99 # (also a 5.10+ feature).
100 if ($] >= 5.010000) {
101 set_prototype \&abs, '_';
102 set_prototype \&cos, '_';
103 set_prototype \&exp, '_';
104 set_prototype \&log, '_';
105 set_prototype \&sin, '_';
106 set_prototype \&sqrt, '_';
107 }
108}
fb73857a 109
9fbe1b12
JH
110my $i;
111my %LOGN;
0c721ce2 112
91cb744f 113# Regular expression for floating point numbers.
bf5f1b4c
JH
114# These days we could use Scalar::Util::lln(), I guess.
115my $gre = qr'\s*([\+\-]?(?:(?:(?:\d+(?:_\d+)*(?:\.\d*(?:_\d+)*)?|\.\d+(?:_\d+)*)(?:[eE][\+\-]?\d+(?:_\d+)*)?))|inf)'i;
91cb744f 116
9fbe1b12 117require Exporter;
0c721ce2 118
9853179e 119our @ISA = qw(Exporter);
5aabfad6 120
5aabfad6
PP
121my @trig = qw(
122 pi
fb73857a 123 tan
5aabfad6
PP
124 csc cosec sec cot cotan
125 asin acos atan
126 acsc acosec asec acot acotan
127 sinh cosh tanh
128 csch cosech sech coth cotanh
129 asinh acosh atanh
130 acsch acosech asech acoth acotanh
131 );
132
9853179e 133our @EXPORT = (qw(
b42d0ec9 134 i Re Im rho theta arg
fb73857a 135 sqrt log ln
5aabfad6
PP
136 log10 logn cbrt root
137 cplx cplxe
bf5f1b4c 138 atan2
5aabfad6
PP
139 ),
140 @trig);
141
1515bec6 142my @pi = qw(pi pi2 pi4 pip2 pip4 Inf);
affad850 143
9853179e 144our @EXPORT_OK = @pi;
bf5f1b4c 145
9853179e 146our %EXPORT_TAGS = (
5aabfad6 147 'trig' => [@trig],
affad850 148 'pi' => [@pi],
66730be0 149);
a0d0e21e 150
a5f75d66 151use overload
3194a40a
Z
152 '=' => \&_copy,
153 '+=' => \&_plus,
affad850 154 '+' => \&_plus,
3194a40a 155 '-=' => \&_minus,
affad850 156 '-' => \&_minus,
3194a40a 157 '*=' => \&_multiply,
affad850 158 '*' => \&_multiply,
3194a40a 159 '/=' => \&_divide,
affad850 160 '/' => \&_divide,
3194a40a 161 '**=' => \&_power,
affad850
SP
162 '**' => \&_power,
163 '==' => \&_numeq,
164 '<=>' => \&_spaceship,
165 'neg' => \&_negate,
166 '~' => \&_conjugate,
66730be0
RM
167 'abs' => \&abs,
168 'sqrt' => \&sqrt,
169 'exp' => \&exp,
170 'log' => \&log,
171 'sin' => \&sin,
172 'cos' => \&cos,
173 'atan2' => \&atan2,
affad850 174 '""' => \&_stringify;
66730be0
RM
175
176#
b42d0ec9 177# Package "privates"
66730be0
RM
178#
179
16357284
JH
180my %DISPLAY_FORMAT = ('style' => 'cartesian',
181 'polar_pretty_print' => 1);
182my $eps = 1e-14; # Epsilon
66730be0
RM
183
184#
185# Object attributes (internal):
186# cartesian [real, imaginary] -- cartesian form
187# polar [rho, theta] -- polar form
188# c_dirty cartesian form not up-to-date
189# p_dirty polar form not up-to-date
190# display display format (package's global when not set)
191#
192
b42d0ec9
JH
193# Die on bad *make() arguments.
194
195sub _cannot_make {
bf5f1b4c 196 die "@{[(caller(1))[3]]}: Cannot take $_[0] of '$_[1]'.\n";
b42d0ec9
JH
197}
198
bf5f1b4c 199sub _make {
91cb744f 200 my $arg = shift;
bf5f1b4c 201 my ($p, $q);
91cb744f 202
bf5f1b4c
JH
203 if ($arg =~ /^$gre$/) {
204 ($p, $q) = ($1, 0);
205 } elsif ($arg =~ /^(?:$gre)?$gre\s*i\s*$/) {
91cb744f 206 ($p, $q) = ($1 || 0, $2);
bf5f1b4c 207 } elsif ($arg =~ /^\s*\(\s*$gre\s*(?:,\s*$gre\s*)?\)\s*$/) {
91cb744f 208 ($p, $q) = ($1, $2 || 0);
91cb744f
JH
209 }
210
bf5f1b4c 211 if (defined $p) {
91cb744f 212 $p =~ s/^\+//;
0720d38b 213 $p =~ s/^(-?)inf$/"${1}9**9**9"/e if $has_inf;
91cb744f 214 $q =~ s/^\+//;
0720d38b 215 $q =~ s/^(-?)inf$/"${1}9**9**9"/e if $has_inf;
91cb744f
JH
216 }
217
bf5f1b4c
JH
218 return ($p, $q);
219}
220
221sub _emake {
222 my $arg = shift;
223 my ($p, $q);
224
225 if ($arg =~ /^\s*\[\s*$gre\s*(?:,\s*$gre\s*)?\]\s*$/) {
226 ($p, $q) = ($1, $2 || 0);
227 } elsif ($arg =~ m!^\s*\[\s*$gre\s*(?:,\s*([-+]?\d*\s*)?pi(?:/\s*(\d+))?\s*)?\]\s*$!) {
228 ($p, $q) = ($1, ($2 eq '-' ? -1 : ($2 || 1)) * pi() / ($3 || 1));
229 } elsif ($arg =~ /^\s*\[\s*$gre\s*\]\s*$/) {
230 ($p, $q) = ($1, 0);
231 } elsif ($arg =~ /^\s*$gre\s*$/) {
232 ($p, $q) = ($1, 0);
233 }
234
235 if (defined $p) {
236 $p =~ s/^\+//;
237 $q =~ s/^\+//;
0720d38b
JH
238 $p =~ s/^(-?)inf$/"${1}9**9**9"/e if $has_inf;
239 $q =~ s/^(-?)inf$/"${1}9**9**9"/e if $has_inf;
bf5f1b4c
JH
240 }
241
242 return ($p, $q);
91cb744f
JH
243}
244
3194a40a
Z
245sub _copy {
246 my $self = shift;
247 my $clone = {%$self};
248 if ($self->{'cartesian'}) {
249 $clone->{'cartesian'} = [@{$self->{'cartesian'}}];
250 }
251 if ($self->{'polar'}) {
252 $clone->{'polar'} = [@{$self->{'polar'}}];
253 }
254 bless $clone,__PACKAGE__;
255 return $clone;
256}
257
66730be0
RM
258#
259# ->make
260#
261# Create a new complex number (cartesian form)
262#
263sub make {
bf5f1b4c
JH
264 my $self = bless {}, shift;
265 my ($re, $im);
266 if (@_ == 0) {
267 ($re, $im) = (0, 0);
268 } elsif (@_ == 1) {
269 return (ref $self)->emake($_[0])
270 if ($_[0] =~ /^\s*\[/);
271 ($re, $im) = _make($_[0]);
272 } elsif (@_ == 2) {
273 ($re, $im) = @_;
274 }
275 if (defined $re) {
91cb744f 276 _cannot_make("real part", $re) unless $re =~ /^$gre$/;
bf5f1b4c
JH
277 }
278 $im ||= 0;
279 _cannot_make("imaginary part", $im) unless $im =~ /^$gre$/;
affad850 280 $self->_set_cartesian([$re, $im ]);
bf5f1b4c
JH
281 $self->display_format('cartesian');
282
283 return $self;
66730be0
RM
284}
285
286#
287# ->emake
288#
289# Create a new complex number (exponential form)
290#
291sub emake {
bf5f1b4c
JH
292 my $self = bless {}, shift;
293 my ($rho, $theta);
294 if (@_ == 0) {
295 ($rho, $theta) = (0, 0);
296 } elsif (@_ == 1) {
297 return (ref $self)->make($_[0])
298 if ($_[0] =~ /^\s*\(/ || $_[0] =~ /i\s*$/);
299 ($rho, $theta) = _emake($_[0]);
300 } elsif (@_ == 2) {
301 ($rho, $theta) = @_;
302 }
303 if (defined $rho && defined $theta) {
fb73857a
PP
304 if ($rho < 0) {
305 $rho = -$rho;
306 $theta = ($theta <= 0) ? $theta + pi() : $theta - pi();
307 }
bf5f1b4c
JH
308 }
309 if (defined $rho) {
91cb744f 310 _cannot_make("rho", $rho) unless $rho =~ /^$gre$/;
bf5f1b4c
JH
311 }
312 $theta ||= 0;
313 _cannot_make("theta", $theta) unless $theta =~ /^$gre$/;
affad850 314 $self->_set_polar([$rho, $theta]);
bf5f1b4c
JH
315 $self->display_format('polar');
316
317 return $self;
66730be0
RM
318}
319
320sub new { &make } # For backward compatibility only.
321
322#
323# cplx
324#
325# Creates a complex number from a (re, im) tuple.
326# This avoids the burden of writing Math::Complex->make(re, im).
327#
328sub cplx {
91cb744f 329 return __PACKAGE__->make(@_);
66730be0
RM
330}
331
332#
333# cplxe
334#
335# Creates a complex number from a (rho, theta) tuple.
336# This avoids the burden of writing Math::Complex->emake(rho, theta).
337#
338sub cplxe {
91cb744f 339 return __PACKAGE__->emake(@_);
66730be0
RM
340}
341
342#
343# pi
344#
fb73857a 345# The number defined as pi = 180 degrees
66730be0 346#
6570f784 347sub pi () { 4 * CORE::atan2(1, 1) }
5cd24f17
PP
348
349#
affad850 350# pi2
5cd24f17 351#
fb73857a
PP
352# The full circle
353#
affad850
SP
354sub pi2 () { 2 * pi }
355
356#
357# pi4
358#
359# The full circle twice.
360#
361sub pi4 () { 4 * pi }
fb73857a 362
5cd24f17 363#
fb73857a
PP
364# pip2
365#
366# The quarter circle
367#
6570f784 368sub pip2 () { pi / 2 }
5cd24f17 369
fb73857a 370#
affad850 371# pip4
d09ae4e6 372#
affad850 373# The eighth circle.
d09ae4e6 374#
affad850 375sub pip4 () { pi / 4 }
d09ae4e6
JH
376
377#
affad850 378# _uplog10
fb73857a
PP
379#
380# Used in log10().
381#
affad850 382sub _uplog10 () { 1 / CORE::log(10) }
66730be0
RM
383
384#
385# i
386#
387# The number defined as i*i = -1;
388#
389sub i () {
5cd24f17
PP
390 return $i if ($i);
391 $i = bless {};
40da2db3 392 $i->{'cartesian'} = [0, 1];
fb73857a 393 $i->{'polar'} = [1, pip2];
66730be0
RM
394 $i->{c_dirty} = 0;
395 $i->{p_dirty} = 0;
396 return $i;
397}
398
399#
affad850 400# _ip2
1fa12f56
JH
401#
402# Half of i.
403#
affad850 404sub _ip2 () { i / 2 }
1fa12f56
JH
405
406#
66730be0
RM
407# Attribute access/set routines
408#
409
affad850
SP
410sub _cartesian {$_[0]->{c_dirty} ?
411 $_[0]->_update_cartesian : $_[0]->{'cartesian'}}
412sub _polar {$_[0]->{p_dirty} ?
413 $_[0]->_update_polar : $_[0]->{'polar'}}
66730be0 414
affad850
SP
415sub _set_cartesian { $_[0]->{p_dirty}++; $_[0]->{c_dirty} = 0;
416 $_[0]->{'cartesian'} = $_[1] }
417sub _set_polar { $_[0]->{c_dirty}++; $_[0]->{p_dirty} = 0;
418 $_[0]->{'polar'} = $_[1] }
66730be0
RM
419
420#
affad850 421# ->_update_cartesian
66730be0
RM
422#
423# Recompute and return the cartesian form, given accurate polar form.
424#
affad850 425sub _update_cartesian {
66730be0 426 my $self = shift;
40da2db3 427 my ($r, $t) = @{$self->{'polar'}};
66730be0 428 $self->{c_dirty} = 0;
a8693bd3 429 return $self->{'cartesian'} = [$r * CORE::cos($t), $r * CORE::sin($t)];
66730be0
RM
430}
431
432#
433#
affad850 434# ->_update_polar
66730be0
RM
435#
436# Recompute and return the polar form, given accurate cartesian form.
437#
affad850 438sub _update_polar {
66730be0 439 my $self = shift;
40da2db3 440 my ($x, $y) = @{$self->{'cartesian'}};
66730be0 441 $self->{p_dirty} = 0;
40da2db3 442 return $self->{'polar'} = [0, 0] if $x == 0 && $y == 0;
1fa12f56
JH
443 return $self->{'polar'} = [CORE::sqrt($x*$x + $y*$y),
444 CORE::atan2($y, $x)];
66730be0
RM
445}
446
447#
affad850 448# (_plus)
66730be0
RM
449#
450# Computes z1+z2.
451#
affad850 452sub _plus {
66730be0 453 my ($z1, $z2, $regular) = @_;
affad850 454 my ($re1, $im1) = @{$z1->_cartesian};
0e505df1 455 $z2 = cplx($z2) unless ref $z2;
affad850 456 my ($re2, $im2) = ref $z2 ? @{$z2->_cartesian} : ($z2, 0);
66730be0 457 unless (defined $regular) {
affad850 458 $z1->_set_cartesian([$re1 + $re2, $im1 + $im2]);
66730be0
RM
459 return $z1;
460 }
461 return (ref $z1)->make($re1 + $re2, $im1 + $im2);
462}
463
464#
affad850 465# (_minus)
66730be0
RM
466#
467# Computes z1-z2.
468#
affad850 469sub _minus {
66730be0 470 my ($z1, $z2, $inverted) = @_;
affad850 471 my ($re1, $im1) = @{$z1->_cartesian};
0e505df1 472 $z2 = cplx($z2) unless ref $z2;
affad850 473 my ($re2, $im2) = @{$z2->_cartesian};
66730be0 474 unless (defined $inverted) {
affad850 475 $z1->_set_cartesian([$re1 - $re2, $im1 - $im2]);
66730be0
RM
476 return $z1;
477 }
478 return $inverted ?
479 (ref $z1)->make($re2 - $re1, $im2 - $im1) :
480 (ref $z1)->make($re1 - $re2, $im1 - $im2);
0e505df1 481
66730be0
RM
482}
483
484#
affad850 485# (_multiply)
66730be0
RM
486#
487# Computes z1*z2.
488#
affad850 489sub _multiply {
fb73857a
PP
490 my ($z1, $z2, $regular) = @_;
491 if ($z1->{p_dirty} == 0 and ref $z2 and $z2->{p_dirty} == 0) {
492 # if both polar better use polar to avoid rounding errors
affad850
SP
493 my ($r1, $t1) = @{$z1->_polar};
494 my ($r2, $t2) = @{$z2->_polar};
fb73857a 495 my $t = $t1 + $t2;
affad850
SP
496 if ($t > pi()) { $t -= pi2 }
497 elsif ($t <= -pi()) { $t += pi2 }
fb73857a 498 unless (defined $regular) {
affad850 499 $z1->_set_polar([$r1 * $r2, $t]);
66730be0 500 return $z1;
fb73857a
PP
501 }
502 return (ref $z1)->emake($r1 * $r2, $t);
503 } else {
affad850 504 my ($x1, $y1) = @{$z1->_cartesian};
fb73857a 505 if (ref $z2) {
affad850 506 my ($x2, $y2) = @{$z2->_cartesian};
fb73857a
PP
507 return (ref $z1)->make($x1*$x2-$y1*$y2, $x1*$y2+$y1*$x2);
508 } else {
509 return (ref $z1)->make($x1*$z2, $y1*$z2);
510 }
66730be0 511 }
66730be0
RM
512}
513
514#
0e505df1 515# _divbyzero
0c721ce2
JH
516#
517# Die on division by zero.
518#
0e505df1 519sub _divbyzero {
5cd24f17
PP
520 my $mess = "$_[0]: Division by zero.\n";
521
522 if (defined $_[1]) {
523 $mess .= "(Because in the definition of $_[0], the divisor ";
1fa12f56 524 $mess .= "$_[1] " unless ("$_[1]" eq '0');
5cd24f17
PP
525 $mess .= "is 0)\n";
526 }
527
0c721ce2 528 my @up = caller(1);
fb73857a 529
5cd24f17
PP
530 $mess .= "Died at $up[1] line $up[2].\n";
531
532 die $mess;
0c721ce2
JH
533}
534
535#
affad850 536# (_divide)
66730be0
RM
537#
538# Computes z1/z2.
539#
affad850 540sub _divide {
66730be0 541 my ($z1, $z2, $inverted) = @_;
fb73857a
PP
542 if ($z1->{p_dirty} == 0 and ref $z2 and $z2->{p_dirty} == 0) {
543 # if both polar better use polar to avoid rounding errors
affad850
SP
544 my ($r1, $t1) = @{$z1->_polar};
545 my ($r2, $t2) = @{$z2->_polar};
fb73857a
PP
546 my $t;
547 if ($inverted) {
0e505df1 548 _divbyzero "$z2/0" if ($r1 == 0);
fb73857a 549 $t = $t2 - $t1;
affad850
SP
550 if ($t > pi()) { $t -= pi2 }
551 elsif ($t <= -pi()) { $t += pi2 }
fb73857a
PP
552 return (ref $z1)->emake($r2 / $r1, $t);
553 } else {
0e505df1 554 _divbyzero "$z1/0" if ($r2 == 0);
fb73857a 555 $t = $t1 - $t2;
affad850
SP
556 if ($t > pi()) { $t -= pi2 }
557 elsif ($t <= -pi()) { $t += pi2 }
fb73857a
PP
558 return (ref $z1)->emake($r1 / $r2, $t);
559 }
560 } else {
561 my ($d, $x2, $y2);
562 if ($inverted) {
affad850 563 ($x2, $y2) = @{$z1->_cartesian};
fb73857a
PP
564 $d = $x2*$x2 + $y2*$y2;
565 _divbyzero "$z2/0" if $d == 0;
566 return (ref $z1)->make(($x2*$z2)/$d, -($y2*$z2)/$d);
567 } else {
affad850 568 my ($x1, $y1) = @{$z1->_cartesian};
fb73857a 569 if (ref $z2) {
affad850 570 ($x2, $y2) = @{$z2->_cartesian};
fb73857a
PP
571 $d = $x2*$x2 + $y2*$y2;
572 _divbyzero "$z1/0" if $d == 0;
573 my $u = ($x1*$x2 + $y1*$y2)/$d;
574 my $v = ($y1*$x2 - $x1*$y2)/$d;
575 return (ref $z1)->make($u, $v);
576 } else {
577 _divbyzero "$z1/0" if $z2 == 0;
578 return (ref $z1)->make($x1/$z2, $y1/$z2);
579 }
580 }
0c721ce2 581 }
66730be0
RM
582}
583
584#
affad850 585# (_power)
66730be0
RM
586#
587# Computes z1**z2 = exp(z2 * log z1)).
588#
affad850 589sub _power {
66730be0 590 my ($z1, $z2, $inverted) = @_;
ace5de91 591 if ($inverted) {
2820d885
DL
592 return 1 if $z1 == 0 || $z2 == 1;
593 return 0 if $z2 == 0 && Re($z1) > 0;
ace5de91 594 } else {
2820d885
DL
595 return 1 if $z2 == 0 || $z1 == 1;
596 return 0 if $z1 == 0 && Re($z2) > 0;
ace5de91 597 }
1fa12f56
JH
598 my $w = $inverted ? &exp($z1 * &log($z2))
599 : &exp($z2 * &log($z1));
d09ae4e6
JH
600 # If both arguments cartesian, return cartesian, else polar.
601 return $z1->{c_dirty} == 0 &&
602 (not ref $z2 or $z2->{c_dirty} == 0) ?
affad850 603 cplx(@{$w->_cartesian}) : $w;
66730be0
RM
604}
605
606#
affad850 607# (_spaceship)
66730be0
RM
608#
609# Computes z1 <=> z2.
2820d885 610# Sorts on the real part first, then on the imaginary part. Thus 2-4i < 3+8i.
66730be0 611#
affad850 612sub _spaceship {
66730be0 613 my ($z1, $z2, $inverted) = @_;
affad850
SP
614 my ($re1, $im1) = ref $z1 ? @{$z1->_cartesian} : ($z1, 0);
615 my ($re2, $im2) = ref $z2 ? @{$z2->_cartesian} : ($z2, 0);
66730be0
RM
616 my $sgn = $inverted ? -1 : 1;
617 return $sgn * ($re1 <=> $re2) if $re1 != $re2;
618 return $sgn * ($im1 <=> $im2);
619}
620
621#
affad850 622# (_numeq)
1fa12f56
JH
623#
624# Computes z1 == z2.
625#
affad850
SP
626# (Required in addition to _spaceship() because of NaNs.)
627sub _numeq {
1fa12f56 628 my ($z1, $z2, $inverted) = @_;
affad850
SP
629 my ($re1, $im1) = ref $z1 ? @{$z1->_cartesian} : ($z1, 0);
630 my ($re2, $im2) = ref $z2 ? @{$z2->_cartesian} : ($z2, 0);
1fa12f56
JH
631 return $re1 == $re2 && $im1 == $im2 ? 1 : 0;
632}
633
634#
affad850 635# (_negate)
66730be0
RM
636#
637# Computes -z.
638#
affad850 639sub _negate {
66730be0
RM
640 my ($z) = @_;
641 if ($z->{c_dirty}) {
affad850 642 my ($r, $t) = @{$z->_polar};
fb73857a
PP
643 $t = ($t <= 0) ? $t + pi : $t - pi;
644 return (ref $z)->emake($r, $t);
66730be0 645 }
affad850 646 my ($re, $im) = @{$z->_cartesian};
66730be0
RM
647 return (ref $z)->make(-$re, -$im);
648}
649
650#
affad850 651# (_conjugate)
66730be0 652#
affad850 653# Compute complex's _conjugate.
66730be0 654#
affad850 655sub _conjugate {
66730be0
RM
656 my ($z) = @_;
657 if ($z->{c_dirty}) {
affad850 658 my ($r, $t) = @{$z->_polar};
66730be0
RM
659 return (ref $z)->emake($r, -$t);
660 }
affad850 661 my ($re, $im) = @{$z->_cartesian};
66730be0
RM
662 return (ref $z)->make($re, -$im);
663}
664
665#
666# (abs)
667#
b42d0ec9 668# Compute or set complex's norm (rho).
66730be0
RM
669#
670sub abs {
e7deb4af 671 my ($z, $rho) = @_ ? @_ : $_;
1fa12f56
JH
672 unless (ref $z) {
673 if (@_ == 2) {
674 $_[0] = $_[1];
675 } else {
676 return CORE::abs($z);
677 }
678 }
b42d0ec9 679 if (defined $rho) {
affad850 680 $z->{'polar'} = [ $rho, ${$z->_polar}[1] ];
b42d0ec9
JH
681 $z->{p_dirty} = 0;
682 $z->{c_dirty} = 1;
683 return $rho;
684 } else {
affad850 685 return ${$z->_polar}[0];
b42d0ec9
JH
686 }
687}
688
689sub _theta {
690 my $theta = $_[0];
691
affad850
SP
692 if ($$theta > pi()) { $$theta -= pi2 }
693 elsif ($$theta <= -pi()) { $$theta += pi2 }
66730be0
RM
694}
695
696#
697# arg
698#
b42d0ec9 699# Compute or set complex's argument (theta).
66730be0
RM
700#
701sub arg {
b42d0ec9
JH
702 my ($z, $theta) = @_;
703 return $z unless ref $z;
704 if (defined $theta) {
705 _theta(\$theta);
affad850 706 $z->{'polar'} = [ ${$z->_polar}[0], $theta ];
b42d0ec9
JH
707 $z->{p_dirty} = 0;
708 $z->{c_dirty} = 1;
709 } else {
affad850 710 $theta = ${$z->_polar}[1];
b42d0ec9
JH
711 _theta(\$theta);
712 }
713 return $theta;
66730be0
RM
714}
715
716#
717# (sqrt)
718#
0c721ce2 719# Compute sqrt(z).
66730be0 720#
b42d0ec9
JH
721# It is quite tempting to use wantarray here so that in list context
722# sqrt() would return the two solutions. This, however, would
723# break things like
724#
725# print "sqrt(z) = ", sqrt($z), "\n";
726#
727# The two values would be printed side by side without no intervening
728# whitespace, quite confusing.
729# Therefore if you want the two solutions use the root().
730#
66730be0 731sub sqrt {
e7deb4af 732 my ($z) = @_ ? $_[0] : $_;
affad850 733 my ($re, $im) = ref $z ? @{$z->_cartesian} : ($z, 0);
1fa12f56
JH
734 return $re < 0 ? cplx(0, CORE::sqrt(-$re)) : CORE::sqrt($re)
735 if $im == 0;
affad850 736 my ($r, $t) = @{$z->_polar};
a8693bd3 737 return (ref $z)->emake(CORE::sqrt($r), $t/2);
66730be0
RM
738}
739
740#
741# cbrt
742#
0c721ce2 743# Compute cbrt(z) (cubic root).
66730be0 744#
b42d0ec9
JH
745# Why are we not returning three values? The same answer as for sqrt().
746#
66730be0
RM
747sub cbrt {
748 my ($z) = @_;
1fa12f56
JH
749 return $z < 0 ?
750 -CORE::exp(CORE::log(-$z)/3) :
751 ($z > 0 ? CORE::exp(CORE::log($z)/3): 0)
fb73857a 752 unless ref $z;
affad850 753 my ($r, $t) = @{$z->_polar};
1fa12f56 754 return 0 if $r == 0;
a8693bd3 755 return (ref $z)->emake(CORE::exp(CORE::log($r)/3), $t/3);
66730be0
RM
756}
757
758#
0e505df1
JH
759# _rootbad
760#
761# Die on bad root.
762#
763sub _rootbad {
bf5f1b4c 764 my $mess = "Root '$_[0]' illegal, root rank must be positive integer.\n";
0e505df1
JH
765
766 my @up = caller(1);
fb73857a 767
0e505df1
JH
768 $mess .= "Died at $up[1] line $up[2].\n";
769
770 die $mess;
771}
772
773#
66730be0
RM
774# root
775#
776# Computes all nth root for z, returning an array whose size is n.
777# `n' must be a positive integer.
778#
779# The roots are given by (for k = 0..n-1):
780#
781# z^(1/n) = r^(1/n) (cos ((t+2 k pi)/n) + i sin ((t+2 k pi)/n))
782#
783sub root {
bf5f1b4c 784 my ($z, $n, $k) = @_;
0e505df1 785 _rootbad($n) if ($n < 1 or int($n) != $n);
1fa12f56 786 my ($r, $t) = ref $z ?
affad850
SP
787 @{$z->_polar} : (CORE::abs($z), $z >= 0 ? 0 : pi);
788 my $theta_inc = pi2 / $n;
66730be0 789 my $rho = $r ** (1/$n);
d09ae4e6 790 my $cartesian = ref $z && $z->{c_dirty} == 0;
bf5f1b4c
JH
791 if (@_ == 2) {
792 my @root;
793 for (my $i = 0, my $theta = $t / $n;
794 $i < $n;
795 $i++, $theta += $theta_inc) {
796 my $w = cplxe($rho, $theta);
797 # Yes, $cartesian is loop invariant.
affad850 798 push @root, $cartesian ? cplx(@{$w->_cartesian}) : $w;
bf5f1b4c
JH
799 }
800 return @root;
801 } elsif (@_ == 3) {
802 my $w = cplxe($rho, $t / $n + $k * $theta_inc);
affad850 803 return $cartesian ? cplx(@{$w->_cartesian}) : $w;
a0d0e21e 804 }
a0d0e21e
LW
805}
806
66730be0
RM
807#
808# Re
809#
b42d0ec9 810# Return or set Re(z).
66730be0 811#
a0d0e21e 812sub Re {
b42d0ec9 813 my ($z, $Re) = @_;
66730be0 814 return $z unless ref $z;
b42d0ec9 815 if (defined $Re) {
affad850 816 $z->{'cartesian'} = [ $Re, ${$z->_cartesian}[1] ];
b42d0ec9
JH
817 $z->{c_dirty} = 0;
818 $z->{p_dirty} = 1;
819 } else {
affad850 820 return ${$z->_cartesian}[0];
b42d0ec9 821 }
a0d0e21e
LW
822}
823
66730be0
RM
824#
825# Im
826#
b42d0ec9 827# Return or set Im(z).
66730be0 828#
a0d0e21e 829sub Im {
b42d0ec9 830 my ($z, $Im) = @_;
178326e7 831 return 0 unless ref $z;
b42d0ec9 832 if (defined $Im) {
affad850 833 $z->{'cartesian'} = [ ${$z->_cartesian}[0], $Im ];
b42d0ec9
JH
834 $z->{c_dirty} = 0;
835 $z->{p_dirty} = 1;
836 } else {
affad850 837 return ${$z->_cartesian}[1];
b42d0ec9
JH
838 }
839}
840
841#
842# rho
843#
844# Return or set rho(w).
845#
846sub rho {
847 Math::Complex::abs(@_);
848}
849
850#
851# theta
852#
853# Return or set theta(w).
854#
855sub theta {
856 Math::Complex::arg(@_);
a0d0e21e
LW
857}
858
66730be0
RM
859#
860# (exp)
861#
862# Computes exp(z).
863#
864sub exp {
e7deb4af
SH
865 my ($z) = @_ ? @_ : $_;
866 return CORE::exp($z) unless ref $z;
867 my ($x, $y) = @{$z->_cartesian};
868 return (ref $z)->emake(CORE::exp($x), $y);
66730be0
RM
869}
870
871#
8c03c583
JH
872# _logofzero
873#
fb73857a 874# Die on logarithm of zero.
8c03c583
JH
875#
876sub _logofzero {
877 my $mess = "$_[0]: Logarithm of zero.\n";
878
879 if (defined $_[1]) {
880 $mess .= "(Because in the definition of $_[0], the argument ";
881 $mess .= "$_[1] " unless ($_[1] eq '0');
882 $mess .= "is 0)\n";
883 }
884
885 my @up = caller(1);
fb73857a 886
8c03c583
JH
887 $mess .= "Died at $up[1] line $up[2].\n";
888
889 die $mess;
890}
891
892#
66730be0
RM
893# (log)
894#
895# Compute log(z).
896#
897sub log {
e7deb4af 898 my ($z) = @_ ? @_ : $_;
fb73857a
PP
899 unless (ref $z) {
900 _logofzero("log") if $z == 0;
a8693bd3 901 return $z > 0 ? CORE::log($z) : cplx(CORE::log(-$z), pi);
fb73857a 902 }
affad850 903 my ($r, $t) = @{$z->_polar};
fb73857a 904 _logofzero("log") if $r == 0;
affad850
SP
905 if ($t > pi()) { $t -= pi2 }
906 elsif ($t <= -pi()) { $t += pi2 }
a8693bd3 907 return (ref $z)->make(CORE::log($r), $t);
66730be0
RM
908}
909
910#
0c721ce2
JH
911# ln
912#
913# Alias for log().
914#
915sub ln { Math::Complex::log(@_) }
916
917#
66730be0
RM
918# log10
919#
920# Compute log10(z).
921#
5cd24f17 922
66730be0 923sub log10 {
affad850 924 return Math::Complex::log($_[0]) * _uplog10;
66730be0
RM
925}
926
927#
928# logn
929#
930# Compute logn(z,n) = log(z) / log(n)
931#
932sub logn {
933 my ($z, $n) = @_;
0c721ce2 934 $z = cplx($z, 0) unless ref $z;
9fbe1b12
JH
935 my $logn = $LOGN{$n};
936 $logn = $LOGN{$n} = CORE::log($n) unless defined $logn; # Cache log(n)
1fa12f56 937 return &log($z) / $logn;
66730be0
RM
938}
939
940#
941# (cos)
942#
943# Compute cos(z) = (exp(iz) + exp(-iz))/2.
944#
945sub cos {
e7deb4af 946 my ($z) = @_ ? @_ : $_;
1fa12f56 947 return CORE::cos($z) unless ref $z;
affad850 948 my ($x, $y) = @{$z->_cartesian};
a8693bd3 949 my $ey = CORE::exp($y);
1fa12f56
JH
950 my $sx = CORE::sin($x);
951 my $cx = CORE::cos($x);
1515bec6 952 my $ey_1 = $ey ? 1 / $ey : Inf();
1fa12f56
JH
953 return (ref $z)->make($cx * ($ey + $ey_1)/2,
954 $sx * ($ey_1 - $ey)/2);
66730be0
RM
955}
956
957#
958# (sin)
959#
960# Compute sin(z) = (exp(iz) - exp(-iz))/2.
961#
962sub sin {
e7deb4af 963 my ($z) = @_ ? @_ : $_;
1fa12f56 964 return CORE::sin($z) unless ref $z;
affad850 965 my ($x, $y) = @{$z->_cartesian};
a8693bd3 966 my $ey = CORE::exp($y);
1fa12f56
JH
967 my $sx = CORE::sin($x);
968 my $cx = CORE::cos($x);
1515bec6 969 my $ey_1 = $ey ? 1 / $ey : Inf();
1fa12f56
JH
970 return (ref $z)->make($sx * ($ey + $ey_1)/2,
971 $cx * ($ey - $ey_1)/2);
66730be0
RM
972}
973
974#
975# tan
976#
977# Compute tan(z) = sin(z) / cos(z).
978#
979sub tan {
980 my ($z) = @_;
1fa12f56
JH
981 my $cz = &cos($z);
982 _divbyzero "tan($z)", "cos($z)" if $cz == 0;
983 return &sin($z) / $cz;
66730be0
RM
984}
985
986#
0c721ce2
JH
987# sec
988#
989# Computes the secant sec(z) = 1 / cos(z).
990#
991sub sec {
992 my ($z) = @_;
1fa12f56 993 my $cz = &cos($z);
0e505df1 994 _divbyzero "sec($z)", "cos($z)" if ($cz == 0);
0c721ce2
JH
995 return 1 / $cz;
996}
997
998#
999# csc
1000#
1001# Computes the cosecant csc(z) = 1 / sin(z).
1002#
1003sub csc {
1004 my ($z) = @_;
1fa12f56 1005 my $sz = &sin($z);
0e505df1 1006 _divbyzero "csc($z)", "sin($z)" if ($sz == 0);
0c721ce2
JH
1007 return 1 / $sz;
1008}
1009
66730be0 1010#
0c721ce2 1011# cosec
66730be0 1012#
0c721ce2
JH
1013# Alias for csc().
1014#
1015sub cosec { Math::Complex::csc(@_) }
1016
1017#
1018# cot
1019#
fb73857a 1020# Computes cot(z) = cos(z) / sin(z).
0c721ce2
JH
1021#
1022sub cot {
66730be0 1023 my ($z) = @_;
1fa12f56 1024 my $sz = &sin($z);
0e505df1 1025 _divbyzero "cot($z)", "sin($z)" if ($sz == 0);
1fa12f56 1026 return &cos($z) / $sz;
66730be0
RM
1027}
1028
1029#
0c721ce2
JH
1030# cotan
1031#
1032# Alias for cot().
1033#
1034sub cotan { Math::Complex::cot(@_) }
1035
1036#
66730be0
RM
1037# acos
1038#
1039# Computes the arc cosine acos(z) = -i log(z + sqrt(z*z-1)).
1040#
1041sub acos {
fb73857a 1042 my $z = $_[0];
1fa12f56
JH
1043 return CORE::atan2(CORE::sqrt(1-$z*$z), $z)
1044 if (! ref $z) && CORE::abs($z) <= 1;
40b904b7 1045 $z = cplx($z, 0) unless ref $z;
affad850 1046 my ($x, $y) = @{$z->_cartesian};
1fa12f56 1047 return 0 if $x == 1 && $y == 0;
a8693bd3
NIS
1048 my $t1 = CORE::sqrt(($x+1)*($x+1) + $y*$y);
1049 my $t2 = CORE::sqrt(($x-1)*($x-1) + $y*$y);
fb73857a
PP
1050 my $alpha = ($t1 + $t2)/2;
1051 my $beta = ($t1 - $t2)/2;
1052 $alpha = 1 if $alpha < 1;
1053 if ($beta > 1) { $beta = 1 }
1054 elsif ($beta < -1) { $beta = -1 }
a8693bd3
NIS
1055 my $u = CORE::atan2(CORE::sqrt(1-$beta*$beta), $beta);
1056 my $v = CORE::log($alpha + CORE::sqrt($alpha*$alpha-1));
fb73857a 1057 $v = -$v if $y > 0 || ($y == 0 && $x < -1);
40b904b7 1058 return (ref $z)->make($u, $v);
66730be0
RM
1059}
1060
1061#
1062# asin
1063#
1064# Computes the arc sine asin(z) = -i log(iz + sqrt(1-z*z)).
1065#
1066sub asin {
fb73857a 1067 my $z = $_[0];
1fa12f56
JH
1068 return CORE::atan2($z, CORE::sqrt(1-$z*$z))
1069 if (! ref $z) && CORE::abs($z) <= 1;
40b904b7 1070 $z = cplx($z, 0) unless ref $z;
affad850 1071 my ($x, $y) = @{$z->_cartesian};
1fa12f56 1072 return 0 if $x == 0 && $y == 0;
a8693bd3
NIS
1073 my $t1 = CORE::sqrt(($x+1)*($x+1) + $y*$y);
1074 my $t2 = CORE::sqrt(($x-1)*($x-1) + $y*$y);
fb73857a
PP
1075 my $alpha = ($t1 + $t2)/2;
1076 my $beta = ($t1 - $t2)/2;
1077 $alpha = 1 if $alpha < 1;
1078 if ($beta > 1) { $beta = 1 }
1079 elsif ($beta < -1) { $beta = -1 }
a8693bd3
NIS
1080 my $u = CORE::atan2($beta, CORE::sqrt(1-$beta*$beta));
1081 my $v = -CORE::log($alpha + CORE::sqrt($alpha*$alpha-1));
fb73857a 1082 $v = -$v if $y > 0 || ($y == 0 && $x < -1);
40b904b7 1083 return (ref $z)->make($u, $v);
66730be0
RM
1084}
1085
1086#
1087# atan
1088#
0c721ce2 1089# Computes the arc tangent atan(z) = i/2 log((i+z) / (i-z)).
66730be0
RM
1090#
1091sub atan {
1092 my ($z) = @_;
a8693bd3 1093 return CORE::atan2($z, 1) unless ref $z;
affad850 1094 my ($x, $y) = ref $z ? @{$z->_cartesian} : ($z, 0);
1fa12f56 1095 return 0 if $x == 0 && $y == 0;
8c03c583 1096 _divbyzero "atan(i)" if ( $z == i);
1fa12f56
JH
1097 _logofzero "atan(-i)" if (-$z == i); # -i is a bad file test...
1098 my $log = &log((i + $z) / (i - $z));
affad850 1099 return _ip2 * $log;
a0d0e21e
LW
1100}
1101
66730be0 1102#
0c721ce2
JH
1103# asec
1104#
1105# Computes the arc secant asec(z) = acos(1 / z).
1106#
1107sub asec {
1108 my ($z) = @_;
0e505df1 1109 _divbyzero "asec($z)", $z if ($z == 0);
fb73857a 1110 return acos(1 / $z);
0c721ce2
JH
1111}
1112
1113#
5cd24f17 1114# acsc
0c721ce2 1115#
8c03c583 1116# Computes the arc cosecant acsc(z) = asin(1 / z).
0c721ce2 1117#
5cd24f17 1118sub acsc {
0c721ce2 1119 my ($z) = @_;
0e505df1 1120 _divbyzero "acsc($z)", $z if ($z == 0);
fb73857a 1121 return asin(1 / $z);
0c721ce2
JH
1122}
1123
1124#
5cd24f17 1125# acosec
66730be0 1126#
5cd24f17 1127# Alias for acsc().
0c721ce2 1128#
5cd24f17 1129sub acosec { Math::Complex::acsc(@_) }
0c721ce2 1130
66730be0 1131#
0c721ce2
JH
1132# acot
1133#
8c03c583 1134# Computes the arc cotangent acot(z) = atan(1 / z)
0c721ce2
JH
1135#
1136sub acot {
66730be0 1137 my ($z) = @_;
1fa12f56
JH
1138 _divbyzero "acot(0)" if $z == 0;
1139 return ($z >= 0) ? CORE::atan2(1, $z) : CORE::atan2(-1, -$z)
1140 unless ref $z;
1141 _divbyzero "acot(i)" if ($z - i == 0);
1142 _logofzero "acot(-i)" if ($z + i == 0);
8c03c583 1143 return atan(1 / $z);
66730be0
RM
1144}
1145
1146#
0c721ce2
JH
1147# acotan
1148#
1149# Alias for acot().
1150#
1151sub acotan { Math::Complex::acot(@_) }
1152
1153#
66730be0
RM
1154# cosh
1155#
1156# Computes the hyperbolic cosine cosh(z) = (exp(z) + exp(-z))/2.
1157#
1158sub cosh {
1159 my ($z) = @_;
fb73857a 1160 my $ex;
0e505df1 1161 unless (ref $z) {
a8693bd3 1162 $ex = CORE::exp($z);
f1e71051 1163 return $ex ? ($ex == $ExpInf ? Inf() : ($ex + 1/$ex)/2) : Inf();
0e505df1 1164 }
affad850 1165 my ($x, $y) = @{$z->_cartesian};
a8693bd3 1166 $ex = CORE::exp($x);
1515bec6 1167 my $ex_1 = $ex ? 1 / $ex : Inf();
a8693bd3
NIS
1168 return (ref $z)->make(CORE::cos($y) * ($ex + $ex_1)/2,
1169 CORE::sin($y) * ($ex - $ex_1)/2);
66730be0
RM
1170}
1171
1172#
1173# sinh
1174#
1175# Computes the hyperbolic sine sinh(z) = (exp(z) - exp(-z))/2.
1176#
1177sub sinh {
1178 my ($z) = @_;
fb73857a 1179 my $ex;
0e505df1 1180 unless (ref $z) {
1fa12f56 1181 return 0 if $z == 0;
a8693bd3 1182 $ex = CORE::exp($z);
f1e71051 1183 return $ex ? ($ex == $ExpInf ? Inf() : ($ex - 1/$ex)/2) : -Inf();
0e505df1 1184 }
affad850 1185 my ($x, $y) = @{$z->_cartesian};
1fa12f56
JH
1186 my $cy = CORE::cos($y);
1187 my $sy = CORE::sin($y);
a8693bd3 1188 $ex = CORE::exp($x);
1515bec6 1189 my $ex_1 = $ex ? 1 / $ex : Inf();
5240e574
JH
1190 return (ref $z)->make(CORE::cos($y) * ($ex - $ex_1)/2,
1191 CORE::sin($y) * ($ex + $ex_1)/2);
66730be0
RM
1192}
1193
1194#
1195# tanh
1196#
1197# Computes the hyperbolic tangent tanh(z) = sinh(z) / cosh(z).
1198#
1199sub tanh {
1200 my ($z) = @_;
0c721ce2 1201 my $cz = cosh($z);
0e505df1 1202 _divbyzero "tanh($z)", "cosh($z)" if ($cz == 0);
1515bec6
SP
1203 my $sz = sinh($z);
1204 return 1 if $cz == $sz;
1205 return -1 if $cz == -$sz;
1206 return $sz / $cz;
66730be0
RM
1207}
1208
1209#
0c721ce2
JH
1210# sech
1211#
1212# Computes the hyperbolic secant sech(z) = 1 / cosh(z).
1213#
1214sub sech {
1215 my ($z) = @_;
1216 my $cz = cosh($z);
0e505df1 1217 _divbyzero "sech($z)", "cosh($z)" if ($cz == 0);
0c721ce2
JH
1218 return 1 / $cz;
1219}
1220
1221#
1222# csch
1223#
1224# Computes the hyperbolic cosecant csch(z) = 1 / sinh(z).
66730be0 1225#
0c721ce2
JH
1226sub csch {
1227 my ($z) = @_;
1228 my $sz = sinh($z);
0e505df1 1229 _divbyzero "csch($z)", "sinh($z)" if ($sz == 0);
0c721ce2
JH
1230 return 1 / $sz;
1231}
1232
1233#
1234# cosech
1235#
1236# Alias for csch().
1237#
1238sub cosech { Math::Complex::csch(@_) }
1239
66730be0 1240#
0c721ce2
JH
1241# coth
1242#
1243# Computes the hyperbolic cotangent coth(z) = cosh(z) / sinh(z).
1244#
1245sub coth {
66730be0 1246 my ($z) = @_;
0c721ce2 1247 my $sz = sinh($z);
1fa12f56 1248 _divbyzero "coth($z)", "sinh($z)" if $sz == 0;
1515bec6
SP
1249 my $cz = cosh($z);
1250 return 1 if $cz == $sz;
1251 return -1 if $cz == -$sz;
1252 return $cz / $sz;
66730be0
RM
1253}
1254
1255#
0c721ce2
JH
1256# cotanh
1257#
1258# Alias for coth().
1259#
1260sub cotanh { Math::Complex::coth(@_) }
1261
1262#
66730be0
RM
1263# acosh
1264#
f1e71051 1265# Computes the area/inverse hyperbolic cosine acosh(z) = log(z + sqrt(z*z-1)).
66730be0
RM
1266#
1267sub acosh {
1268 my ($z) = @_;
fb73857a 1269 unless (ref $z) {
fb73857a
PP
1270 $z = cplx($z, 0);
1271 }
affad850 1272 my ($re, $im) = @{$z->_cartesian};
fb73857a 1273 if ($im == 0) {
1fa12f56
JH
1274 return CORE::log($re + CORE::sqrt($re*$re - 1))
1275 if $re >= 1;
1276 return cplx(0, CORE::atan2(CORE::sqrt(1 - $re*$re), $re))
1277 if CORE::abs($re) < 1;
fb73857a 1278 }
9bc5fa8d 1279 my $t = &sqrt($z * $z - 1) + $z;
40b904b7
JH
1280 # Try Taylor if looking bad (this usually means that
1281 # $z was large negative, therefore the sqrt is really
1282 # close to abs(z), summing that with z...)
9bc5fa8d
JH
1283 $t = 1/(2 * $z) - 1/(8 * $z**3) + 1/(16 * $z**5) - 5/(128 * $z**7)
1284 if $t == 0;
1285 my $u = &log($t);
40b904b7 1286 $u->Im(-$u->Im) if $re < 0 && $im == 0;
9bc5fa8d 1287 return $re < 0 ? -$u : $u;
66730be0
RM
1288}
1289
1290#
1291# asinh
1292#
f1e71051 1293# Computes the area/inverse hyperbolic sine asinh(z) = log(z + sqrt(z*z+1))
66730be0
RM
1294#
1295sub asinh {
1296 my ($z) = @_;
1fa12f56
JH
1297 unless (ref $z) {
1298 my $t = $z + CORE::sqrt($z*$z + 1);
1299 return CORE::log($t) if $t;
1300 }
9bc5fa8d 1301 my $t = &sqrt($z * $z + 1) + $z;
40b904b7
JH
1302 # Try Taylor if looking bad (this usually means that
1303 # $z was large negative, therefore the sqrt is really
1304 # close to abs(z), summing that with z...)
9bc5fa8d
JH
1305 $t = 1/(2 * $z) - 1/(8 * $z**3) + 1/(16 * $z**5) - 5/(128 * $z**7)
1306 if $t == 0;
1fa12f56 1307 return &log($t);
66730be0
RM
1308}
1309
1310#
1311# atanh
1312#
f1e71051 1313# Computes the area/inverse hyperbolic tangent atanh(z) = 1/2 log((1+z) / (1-z)).
66730be0
RM
1314#
1315sub atanh {
1316 my ($z) = @_;
fb73857a 1317 unless (ref $z) {
a8693bd3 1318 return CORE::log((1 + $z)/(1 - $z))/2 if CORE::abs($z) < 1;
fb73857a
PP
1319 $z = cplx($z, 0);
1320 }
1fa12f56
JH
1321 _divbyzero 'atanh(1)', "1 - $z" if (1 - $z == 0);
1322 _logofzero 'atanh(-1)' if (1 + $z == 0);
1323 return 0.5 * &log((1 + $z) / (1 - $z));
66730be0
RM
1324}
1325
1326#
0c721ce2
JH
1327# asech
1328#
f1e71051 1329# Computes the area/inverse hyperbolic secant asech(z) = acosh(1 / z).
0c721ce2
JH
1330#
1331sub asech {
1332 my ($z) = @_;
1fa12f56 1333 _divbyzero 'asech(0)', "$z" if ($z == 0);
0c721ce2
JH
1334 return acosh(1 / $z);
1335}
1336
1337#
1338# acsch
66730be0 1339#
f1e71051 1340# Computes the area/inverse hyperbolic cosecant acsch(z) = asinh(1 / z).
66730be0 1341#
0c721ce2 1342sub acsch {
66730be0 1343 my ($z) = @_;
0e505df1 1344 _divbyzero 'acsch(0)', $z if ($z == 0);
0c721ce2
JH
1345 return asinh(1 / $z);
1346}
1347
1348#
1349# acosech
1350#
1351# Alias for acosh().
1352#
1353sub acosech { Math::Complex::acsch(@_) }
1354
1355#
1356# acoth
1357#
f1e71051 1358# Computes the area/inverse hyperbolic cotangent acoth(z) = 1/2 log((1+z) / (z-1)).
0c721ce2
JH
1359#
1360sub acoth {
1361 my ($z) = @_;
1fa12f56 1362 _divbyzero 'acoth(0)' if ($z == 0);
fb73857a 1363 unless (ref $z) {
a8693bd3 1364 return CORE::log(($z + 1)/($z - 1))/2 if CORE::abs($z) > 1;
fb73857a
PP
1365 $z = cplx($z, 0);
1366 }
1fa12f56
JH
1367 _divbyzero 'acoth(1)', "$z - 1" if ($z - 1 == 0);
1368 _logofzero 'acoth(-1)', "1 + $z" if (1 + $z == 0);
1369 return &log((1 + $z) / ($z - 1)) / 2;
66730be0
RM
1370}
1371
1372#
0c721ce2
JH
1373# acotanh
1374#
1375# Alias for acot().
1376#
1377sub acotanh { Math::Complex::acoth(@_) }
1378
1379#
66730be0
RM
1380# (atan2)
1381#
bf5f1b4c 1382# Compute atan(z1/z2), minding the right quadrant.
66730be0
RM
1383#
1384sub atan2 {
1385 my ($z1, $z2, $inverted) = @_;
fb73857a
PP
1386 my ($re1, $im1, $re2, $im2);
1387 if ($inverted) {
affad850
SP
1388 ($re1, $im1) = ref $z2 ? @{$z2->_cartesian} : ($z2, 0);
1389 ($re2, $im2) = ref $z1 ? @{$z1->_cartesian} : ($z1, 0);
66730be0 1390 } else {
affad850
SP
1391 ($re1, $im1) = ref $z1 ? @{$z1->_cartesian} : ($z1, 0);
1392 ($re2, $im2) = ref $z2 ? @{$z2->_cartesian} : ($z2, 0);
fb73857a 1393 }
bf5f1b4c
JH
1394 if ($im1 || $im2) {
1395 # In MATLAB the imaginary parts are ignored.
1396 # warn "atan2: Imaginary parts ignored";
1397 # http://documents.wolfram.com/mathematica/functions/ArcTan
1398 # NOTE: Mathematica ArcTan[x,y] while atan2(y,x)
1399 my $s = $z1 * $z1 + $z2 * $z2;
1400 _divbyzero("atan2") if $s == 0;
1401 my $i = &i;
1402 my $r = $z2 + $z1 * $i;
1403 return -$i * &log($r / &sqrt( $s ));
66730be0 1404 }
bf5f1b4c 1405 return CORE::atan2($re1, $re2);
66730be0
RM
1406}
1407
1408#
1409# display_format
1410# ->display_format
1411#
16357284 1412# Set (get if no argument) the display format for all complex numbers that
fb73857a 1413# don't happen to have overridden it via ->display_format
66730be0 1414#
16357284 1415# When called as an object method, this actually sets the display format for
66730be0
RM
1416# the current object.
1417#
1418# Valid object formats are 'c' and 'p' for cartesian and polar. The first
1419# letter is used actually, so the type can be fully spelled out for clarity.
1420#
1421sub display_format {
16357284
JH
1422 my $self = shift;
1423 my %display_format = %DISPLAY_FORMAT;
66730be0 1424
16357284
JH
1425 if (ref $self) { # Called as an object method
1426 if (exists $self->{display_format}) {
1427 my %obj = %{$self->{display_format}};
1428 @display_format{keys %obj} = values %obj;
1429 }
476757f7
YN
1430 }
1431 if (@_ == 1) {
1432 $display_format{style} = shift;
1433 } else {
1434 my %new = @_;
1435 @display_format{keys %new} = values %new;
66730be0
RM
1436 }
1437
476757f7 1438 if (ref $self) { # Called as an object method
16357284
JH
1439 $self->{display_format} = { %display_format };
1440 return
1441 wantarray ?
1442 %{$self->{display_format}} :
1443 $self->{display_format}->{style};
66730be0
RM
1444 }
1445
476757f7 1446 # Called as a class method
16357284
JH
1447 %DISPLAY_FORMAT = %display_format;
1448 return
1449 wantarray ?
1450 %DISPLAY_FORMAT :
1451 $DISPLAY_FORMAT{style};
66730be0
RM
1452}
1453
1454#
affad850 1455# (_stringify)
66730be0
RM
1456#
1457# Show nicely formatted complex number under its cartesian or polar form,
1458# depending on the current display format:
1459#
1460# . If a specific display format has been recorded for this object, use it.
1461# . Otherwise, use the generic current default for all complex numbers,
1462# which is a package global variable.
1463#
affad850 1464sub _stringify {
66730be0 1465 my ($z) = shift;
66730be0 1466
16357284
JH
1467 my $style = $z->display_format;
1468
1469 $style = $DISPLAY_FORMAT{style} unless defined $style;
66730be0 1470
affad850
SP
1471 return $z->_stringify_polar if $style =~ /^p/i;
1472 return $z->_stringify_cartesian;
66730be0
RM
1473}
1474
1475#
affad850 1476# ->_stringify_cartesian
66730be0
RM
1477#
1478# Stringify as a cartesian representation 'a+bi'.
1479#
affad850 1480sub _stringify_cartesian {
66730be0 1481 my $z = shift;
affad850 1482 my ($x, $y) = @{$z->_cartesian};
66730be0
RM
1483 my ($re, $im);
1484
16357284
JH
1485 my %format = $z->display_format;
1486 my $format = $format{format};
1487
1fa12f56
JH
1488 if ($x) {
1489 if ($x =~ /^NaN[QS]?$/i) {
1490 $re = $x;
1491 } else {
b57c8994 1492 if ($x =~ /^-?\Q$Inf\E$/oi) {
1fa12f56
JH
1493 $re = $x;
1494 } else {
1495 $re = defined $format ? sprintf($format, $x) : $x;
1496 }
1497 }
1498 } else {
1499 undef $re;
1500 }
1501
1502 if ($y) {
40b904b7 1503 if ($y =~ /^(NaN[QS]?)$/i) {
1fa12f56
JH
1504 $im = $y;
1505 } else {
b57c8994 1506 if ($y =~ /^-?\Q$Inf\E$/oi) {
1fa12f56
JH
1507 $im = $y;
1508 } else {
40b904b7
JH
1509 $im =
1510 defined $format ?
1511 sprintf($format, $y) :
1512 ($y == 1 ? "" : ($y == -1 ? "-" : $y));
1fa12f56
JH
1513 }
1514 }
1515 $im .= "i";
1516 } else {
1517 undef $im;
16357284 1518 }
66730be0 1519
1fa12f56
JH
1520 my $str = $re;
1521
16357284
JH
1522 if (defined $im) {
1523 if ($y < 0) {
1524 $str .= $im;
1fa12f56 1525 } elsif ($y > 0 || $im =~ /^NaN[QS]?i$/i) {
16357284
JH
1526 $str .= "+" if defined $re;
1527 $str .= $im;
1528 }
1fa12f56
JH
1529 } elsif (!defined $re) {
1530 $str = "0";
16357284 1531 }
66730be0
RM
1532
1533 return $str;
1534}
1535
d09ae4e6 1536
66730be0 1537#
affad850 1538# ->_stringify_polar
66730be0
RM
1539#
1540# Stringify as a polar representation '[r,t]'.
1541#
affad850 1542sub _stringify_polar {
66730be0 1543 my $z = shift;
affad850 1544 my ($r, $t) = @{$z->_polar};
66730be0
RM
1545 my $theta;
1546
16357284 1547 my %format = $z->display_format;
1fa12f56 1548 my $format = $format{format};
16357284 1549
b57c8994 1550 if ($t =~ /^NaN[QS]?$/i || $t =~ /^-?\Q$Inf\E$/oi) {
1fa12f56
JH
1551 $theta = $t;
1552 } elsif ($t == pi) {
1553 $theta = "pi";
1554 } elsif ($r == 0 || $t == 0) {
1555 $theta = defined $format ? sprintf($format, $t) : $t;
55497cff 1556 }
66730be0 1557
1fa12f56
JH
1558 return "[$r,$theta]" if defined $theta;
1559
66730be0 1560 #
1fa12f56 1561 # Try to identify pi/n and friends.
66730be0
RM
1562 #
1563
affad850 1564 $t -= int(CORE::abs($t) / pi2) * pi2;
1fa12f56 1565
e97e26fa 1566 if ($format{polar_pretty_print} && $t) {
1fa12f56 1567 my ($a, $b);
9bc5fa8d 1568 for $a (2..9) {
1fa12f56 1569 $b = $t * $a / pi;
e97e26fa 1570 if ($b =~ /^-?\d+$/) {
1fa12f56
JH
1571 $b = $b < 0 ? "-" : "" if CORE::abs($b) == 1;
1572 $theta = "${b}pi/$a";
d09ae4e6 1573 last;
66730be0 1574 }
d09ae4e6 1575 }
66730be0
RM
1576 }
1577
16357284
JH
1578 if (defined $format) {
1579 $r = sprintf($format, $r);
3194a40a 1580 $theta = sprintf($format, $t) unless defined $theta;
1fa12f56
JH
1581 } else {
1582 $theta = $t unless defined $theta;
16357284
JH
1583 }
1584
1fa12f56 1585 return "[$r,$theta]";
a0d0e21e 1586}
a5f75d66 1587
1515bec6
SP
1588sub Inf {
1589 return $Inf;
1590}
1591
a5f75d66
AD
15921;
1593__END__
1594
1cf6bcb8
JH
1595=pod
1596
a5f75d66
AD
1597=head1 NAME
1598
66730be0 1599Math::Complex - complex numbers and associated mathematical functions
a5f75d66
AD
1600
1601=head1 SYNOPSIS
1602
66730be0 1603 use Math::Complex;
fb73857a 1604
66730be0
RM
1605 $z = Math::Complex->make(5, 6);
1606 $t = 4 - 3*i + $z;
1607 $j = cplxe(1, 2*pi/3);
a5f75d66
AD
1608
1609=head1 DESCRIPTION
1610
66730be0
RM
1611This package lets you create and manipulate complex numbers. By default,
1612I<Perl> limits itself to real numbers, but an extra C<use> statement brings
1613full complex support, along with a full set of mathematical functions
1614typically associated with and/or extended to complex numbers.
1615
1616If you wonder what complex numbers are, they were invented to be able to solve
1617the following equation:
1618
1619 x*x = -1
1620
1621and by definition, the solution is noted I<i> (engineers use I<j> instead since
1622I<i> usually denotes an intensity, but the name does not matter). The number
1623I<i> is a pure I<imaginary> number.
1624
1625The arithmetics with pure imaginary numbers works just like you would expect
1626it with real numbers... you just have to remember that
1627
1628 i*i = -1
1629
1630so you have:
1631
1632 5i + 7i = i * (5 + 7) = 12i
1633 4i - 3i = i * (4 - 3) = i
1634 4i * 2i = -8
1635 6i / 2i = 3
1636 1 / i = -i
1637
1638Complex numbers are numbers that have both a real part and an imaginary
1639part, and are usually noted:
1640
1641 a + bi
1642
1643where C<a> is the I<real> part and C<b> is the I<imaginary> part. The
1644arithmetic with complex numbers is straightforward. You have to
1645keep track of the real and the imaginary parts, but otherwise the
1646rules used for real numbers just apply:
1647
1648 (4 + 3i) + (5 - 2i) = (4 + 5) + i(3 - 2) = 9 + i
1649 (2 + i) * (4 - i) = 2*4 + 4i -2i -i*i = 8 + 2i + 1 = 9 + 2i
1650
1651A graphical representation of complex numbers is possible in a plane
1652(also called the I<complex plane>, but it's really a 2D plane).
1653The number
1654
1655 z = a + bi
1656
1657is the point whose coordinates are (a, b). Actually, it would
1658be the vector originating from (0, 0) to (a, b). It follows that the addition
1659of two complex numbers is a vectorial addition.
1660
1661Since there is a bijection between a point in the 2D plane and a complex
1662number (i.e. the mapping is unique and reciprocal), a complex number
1663can also be uniquely identified with polar coordinates:
1664
1665 [rho, theta]
1666
1667where C<rho> is the distance to the origin, and C<theta> the angle between
1668the vector and the I<x> axis. There is a notation for this using the
1669exponential form, which is:
1670
1671 rho * exp(i * theta)
1672
1673where I<i> is the famous imaginary number introduced above. Conversion
1674between this form and the cartesian form C<a + bi> is immediate:
1675
1676 a = rho * cos(theta)
1677 b = rho * sin(theta)
1678
1679which is also expressed by this formula:
1680
fb73857a 1681 z = rho * exp(i * theta) = rho * (cos theta + i * sin theta)
66730be0
RM
1682
1683In other words, it's the projection of the vector onto the I<x> and I<y>
1684axes. Mathematicians call I<rho> the I<norm> or I<modulus> and I<theta>
affad850
SP
1685the I<argument> of the complex number. The I<norm> of C<z> is
1686marked here as C<abs(z)>.
66730be0 1687
affad850
SP
1688The polar notation (also known as the trigonometric representation) is
1689much more handy for performing multiplications and divisions of
1690complex numbers, whilst the cartesian notation is better suited for
1691additions and subtractions. Real numbers are on the I<x> axis, and
1692therefore I<y> or I<theta> is zero or I<pi>.
66730be0
RM
1693
1694All the common operations that can be performed on a real number have
1695been defined to work on complex numbers as well, and are merely
1696I<extensions> of the operations defined on real numbers. This means
1697they keep their natural meaning when there is no imaginary part, provided
1698the number is within their definition set.
1699
1700For instance, the C<sqrt> routine which computes the square root of
fb73857a
PP
1701its argument is only defined for non-negative real numbers and yields a
1702non-negative real number (it is an application from B<R+> to B<R+>).
66730be0
RM
1703If we allow it to return a complex number, then it can be extended to
1704negative real numbers to become an application from B<R> to B<C> (the
1705set of complex numbers):
1706
1707 sqrt(x) = x >= 0 ? sqrt(x) : sqrt(-x)*i
1708
1709It can also be extended to be an application from B<C> to B<C>,
1710whilst its restriction to B<R> behaves as defined above by using
1711the following definition:
1712
1713 sqrt(z = [r,t]) = sqrt(r) * exp(i * t/2)
1714
fb73857a
PP
1715Indeed, a negative real number can be noted C<[x,pi]> (the modulus
1716I<x> is always non-negative, so C<[x,pi]> is really C<-x>, a negative
1717number) and the above definition states that
66730be0
RM
1718
1719 sqrt([x,pi]) = sqrt(x) * exp(i*pi/2) = [sqrt(x),pi/2] = sqrt(x)*i
1720
1721which is exactly what we had defined for negative real numbers above.
b42d0ec9
JH
1722The C<sqrt> returns only one of the solutions: if you want the both,
1723use the C<root> function.
a5f75d66 1724
66730be0
RM
1725All the common mathematical functions defined on real numbers that
1726are extended to complex numbers share that same property of working
1727I<as usual> when the imaginary part is zero (otherwise, it would not
1728be called an extension, would it?).
a5f75d66 1729
66730be0
RM
1730A I<new> operation possible on a complex number that is
1731the identity for real numbers is called the I<conjugate>, and is noted
d1be9408 1732with a horizontal bar above the number, or C<~z> here.
a5f75d66 1733
66730be0
RM
1734 z = a + bi
1735 ~z = a - bi
a5f75d66 1736
66730be0 1737Simple... Now look:
a5f75d66 1738
66730be0 1739 z * ~z = (a + bi) * (a - bi) = a*a + b*b
a5f75d66 1740
66730be0
RM
1741We saw that the norm of C<z> was noted C<abs(z)> and was defined as the
1742distance to the origin, also known as:
a5f75d66 1743
66730be0 1744 rho = abs(z) = sqrt(a*a + b*b)
a5f75d66 1745
66730be0
RM
1746so
1747
1748 z * ~z = abs(z) ** 2
1749
1750If z is a pure real number (i.e. C<b == 0>), then the above yields:
1751
1752 a * a = abs(a) ** 2
1753
1754which is true (C<abs> has the regular meaning for real number, i.e. stands
1755for the absolute value). This example explains why the norm of C<z> is
1756noted C<abs(z)>: it extends the C<abs> function to complex numbers, yet
1757is the regular C<abs> we know when the complex number actually has no
1758imaginary part... This justifies I<a posteriori> our use of the C<abs>
1759notation for the norm.
1760
1761=head1 OPERATIONS
1762
1763Given the following notations:
1764
1765 z1 = a + bi = r1 * exp(i * t1)
1766 z2 = c + di = r2 * exp(i * t2)
1767 z = <any complex or real number>
1768
1769the following (overloaded) operations are supported on complex numbers:
1770
1771 z1 + z2 = (a + c) + i(b + d)
1772 z1 - z2 = (a - c) + i(b - d)
1773 z1 * z2 = (r1 * r2) * exp(i * (t1 + t2))
1774 z1 / z2 = (r1 / r2) * exp(i * (t1 - t2))
1775 z1 ** z2 = exp(z2 * log z1)
b42d0ec9
JH
1776 ~z = a - bi
1777 abs(z) = r1 = sqrt(a*a + b*b)
1778 sqrt(z) = sqrt(r1) * exp(i * t/2)
1779 exp(z) = exp(a) * exp(i * b)
1780 log(z) = log(r1) + i*t
1781 sin(z) = 1/2i (exp(i * z1) - exp(-i * z))
1782 cos(z) = 1/2 (exp(i * z1) + exp(-i * z))
bf5f1b4c
JH
1783 atan2(y, x) = atan(y / x) # Minding the right quadrant, note the order.
1784
1785The definition used for complex arguments of atan2() is
1786
1787 -i log((x + iy)/sqrt(x*x+y*y))
66730be0 1788
affad850
SP
1789Note that atan2(0, 0) is not well-defined.
1790
66730be0
RM
1791The following extra operations are supported on both real and complex
1792numbers:
1793
1794 Re(z) = a
1795 Im(z) = b
1796 arg(z) = t
b42d0ec9 1797 abs(z) = r
66730be0
RM
1798
1799 cbrt(z) = z ** (1/3)
1800 log10(z) = log(z) / log(10)
1801 logn(z, n) = log(z) / log(n)
1802
1803 tan(z) = sin(z) / cos(z)
0c721ce2 1804
5aabfad6
PP
1805 csc(z) = 1 / sin(z)
1806 sec(z) = 1 / cos(z)
0c721ce2 1807 cot(z) = 1 / tan(z)
66730be0
RM
1808
1809 asin(z) = -i * log(i*z + sqrt(1-z*z))
fb73857a 1810 acos(z) = -i * log(z + i*sqrt(1-z*z))
66730be0 1811 atan(z) = i/2 * log((i+z) / (i-z))
0c721ce2 1812
5aabfad6
PP
1813 acsc(z) = asin(1 / z)
1814 asec(z) = acos(1 / z)
8c03c583 1815 acot(z) = atan(1 / z) = -i/2 * log((i+z) / (z-i))
66730be0
RM
1816
1817 sinh(z) = 1/2 (exp(z) - exp(-z))
1818 cosh(z) = 1/2 (exp(z) + exp(-z))
0c721ce2
JH
1819 tanh(z) = sinh(z) / cosh(z) = (exp(z) - exp(-z)) / (exp(z) + exp(-z))
1820
5aabfad6
PP
1821 csch(z) = 1 / sinh(z)
1822 sech(z) = 1 / cosh(z)
0c721ce2 1823 coth(z) = 1 / tanh(z)
fb73857a 1824
66730be0
RM
1825 asinh(z) = log(z + sqrt(z*z+1))
1826 acosh(z) = log(z + sqrt(z*z-1))
1827 atanh(z) = 1/2 * log((1+z) / (1-z))
66730be0 1828
5aabfad6
PP
1829 acsch(z) = asinh(1 / z)
1830 asech(z) = acosh(1 / z)
0c721ce2
JH
1831 acoth(z) = atanh(1 / z) = 1/2 * log((1+z) / (z-1))
1832
b42d0ec9
JH
1833I<arg>, I<abs>, I<log>, I<csc>, I<cot>, I<acsc>, I<acot>, I<csch>,
1834I<coth>, I<acosech>, I<acotanh>, have aliases I<rho>, I<theta>, I<ln>,
1835I<cosec>, I<cotan>, I<acosec>, I<acotan>, I<cosech>, I<cotanh>,
1836I<acosech>, I<acotanh>, respectively. C<Re>, C<Im>, C<arg>, C<abs>,
d1be9408 1837C<rho>, and C<theta> can be used also as mutators. The C<cbrt>
b42d0ec9
JH
1838returns only one of the solutions: if you want all three, use the
1839C<root> function.
0c721ce2
JH
1840
1841The I<root> function is available to compute all the I<n>
66730be0
RM
1842roots of some complex, where I<n> is a strictly positive integer.
1843There are exactly I<n> such roots, returned as a list. Getting the
1844number mathematicians call C<j> such that:
1845
1846 1 + j + j*j = 0;
1847
1848is a simple matter of writing:
1849
1850 $j = ((root(1, 3))[1];
1851
1852The I<k>th root for C<z = [r,t]> is given by:
1853
1854 (root(z, n))[k] = r**(1/n) * exp(i * (t + 2*k*pi)/n)
1855
bf5f1b4c
JH
1856You can return the I<k>th root directly by C<root(z, n, k)>,
1857indexing starting from I<zero> and ending at I<n - 1>.
1858
1515bec6
SP
1859The I<spaceship> numeric comparison operator, E<lt>=E<gt>, is also
1860defined. In order to ensure its restriction to real numbers is conform
1861to what you would expect, the comparison is run on the real part of
1862the complex number first, and imaginary parts are compared only when
1863the real parts match.
66730be0
RM
1864
1865=head1 CREATION
1866
1867To create a complex number, use either:
1868
1869 $z = Math::Complex->make(3, 4);
1870 $z = cplx(3, 4);
1871
1872if you know the cartesian form of the number, or
1873
1874 $z = 3 + 4*i;
1875
fb73857a 1876if you like. To create a number using the polar form, use either:
66730be0
RM
1877
1878 $z = Math::Complex->emake(5, pi/3);
1879 $x = cplxe(5, pi/3);
1880
0c721ce2 1881instead. The first argument is the modulus, the second is the angle
fb73857a
PP
1882(in radians, the full circle is 2*pi). (Mnemonic: C<e> is used as a
1883notation for complex numbers in the polar form).
66730be0
RM
1884
1885It is possible to write:
1886
1887 $x = cplxe(-3, pi/4);
1888
16357284
JH
1889but that will be silently converted into C<[3,-3pi/4]>, since the
1890modulus must be non-negative (it represents the distance to the origin
1891in the complex plane).
66730be0 1892
91cb744f
JH
1893It is also possible to have a complex number as either argument of the
1894C<make>, C<emake>, C<cplx>, and C<cplxe>: the appropriate component of
b42d0ec9
JH
1895the argument will be used.
1896
1897 $z1 = cplx(-2, 1);
1898 $z2 = cplx($z1, 4);
1899
91cb744f
JH
1900The C<new>, C<make>, C<emake>, C<cplx>, and C<cplxe> will also
1901understand a single (string) argument of the forms
1902
1903 2-3i
1904 -3i
1905 [2,3]
bf5f1b4c 1906 [2,-3pi/4]
91cb744f
JH
1907 [2]
1908
1909in which case the appropriate cartesian and exponential components
1910will be parsed from the string and used to create new complex numbers.
1911The imaginary component and the theta, respectively, will default to zero.
1912
bf5f1b4c
JH
1913The C<new>, C<make>, C<emake>, C<cplx>, and C<cplxe> will also
1914understand the case of no arguments: this means plain zero or (0, 0).
1915
1916=head1 DISPLAYING
66730be0
RM
1917
1918When printed, a complex number is usually shown under its cartesian
16357284 1919style I<a+bi>, but there are legitimate cases where the polar style
bf5f1b4c
JH
1920I<[r,t]> is more appropriate. The process of converting the complex
1921number into a string that can be displayed is known as I<stringification>.
66730be0 1922
16357284
JH
1923By calling the class method C<Math::Complex::display_format> and
1924supplying either C<"polar"> or C<"cartesian"> as an argument, you
5287f86b 1925override the default display style, which is C<"cartesian">. Not
16357284 1926supplying any argument returns the current settings.
66730be0
RM
1927
1928This default can be overridden on a per-number basis by calling the
1929C<display_format> method instead. As before, not supplying any argument
5287f86b
JH
1930returns the current display style for this number. Otherwise whatever you
1931specify will be the new display style for I<this> particular number.
66730be0
RM
1932
1933For instance:
1934
1935 use Math::Complex;
1936
1937 Math::Complex::display_format('polar');
16357284
JH
1938 $j = (root(1, 3))[1];
1939 print "j = $j\n"; # Prints "j = [1,2pi/3]"
66730be0
RM
1940 $j->display_format('cartesian');
1941 print "j = $j\n"; # Prints "j = -0.5+0.866025403784439i"
1942
5287f86b 1943The polar style attempts to emphasize arguments like I<k*pi/n>
9bc5fa8d 1944(where I<n> is a positive integer and I<k> an integer within [-9, +9]),
5287f86b 1945this is called I<polar pretty-printing>.
66730be0 1946
bf5f1b4c
JH
1947For the reverse of stringifying, see the C<make> and C<emake>.
1948
16357284
JH
1949=head2 CHANGED IN PERL 5.6
1950
1951The C<display_format> class method and the corresponding
1952C<display_format> object method can now be called using
1953a parameter hash instead of just a one parameter.
1954
1955The old display format style, which can have values C<"cartesian"> or
40b904b7
JH
1956C<"polar">, can be changed using the C<"style"> parameter.
1957
1958 $j->display_format(style => "polar");
1959
1960The one parameter calling convention also still works.
1961
1962 $j->display_format("polar");
16357284
JH
1963
1964There are two new display parameters.
1965
40b904b7
JH
1966The first one is C<"format">, which is a sprintf()-style format string
1967to be used for both numeric parts of the complex number(s). The is
1968somewhat system-dependent but most often it corresponds to C<"%.15g">.
1969You can revert to the default by setting the C<format> to C<undef>.
16357284
JH
1970
1971 # the $j from the above example
1972
1973 $j->display_format('format' => '%.5f');
1974 print "j = $j\n"; # Prints "j = -0.50000+0.86603i"
40b904b7 1975 $j->display_format('format' => undef);
16357284
JH
1976 print "j = $j\n"; # Prints "j = -0.5+0.86603i"
1977
1978Notice that this affects also the return values of the
1979C<display_format> methods: in list context the whole parameter hash
40b904b7
JH
1980will be returned, as opposed to only the style parameter value.
1981This is a potential incompatibility with earlier versions if you
1982have been calling the C<display_format> method in list context.
16357284 1983
5287f86b
JH
1984The second new display parameter is C<"polar_pretty_print">, which can
1985be set to true or false, the default being true. See the previous
1986section for what this means.
16357284 1987
66730be0
RM
1988=head1 USAGE
1989
1990Thanks to overloading, the handling of arithmetics with complex numbers
1991is simple and almost transparent.
1992
1993Here are some examples:
1994
1995 use Math::Complex;
1996
1997 $j = cplxe(1, 2*pi/3); # $j ** 3 == 1
1998 print "j = $j, j**3 = ", $j ** 3, "\n";
1999 print "1 + j + j**2 = ", 1 + $j + $j**2, "\n";
2000
2001 $z = -16 + 0*i; # Force it to be a complex
2002 print "sqrt($z) = ", sqrt($z), "\n";
2003
2004 $k = exp(i * 2*pi/3);
2005 print "$j - $k = ", $j - $k, "\n";
a5f75d66 2006
b42d0ec9
JH
2007 $z->Re(3); # Re, Im, arg, abs,
2008 $j->arg(2); # (the last two aka rho, theta)
2009 # can be used also as mutators.
2010
7637cd07
SP
2011=head1 CONSTANTS
2012
affad850
SP
2013=head2 PI
2014
2015The constant C<pi> and some handy multiples of it (pi2, pi4,
2016and pip2 (pi/2) and pip4 (pi/4)) are also available if separately
2017exported:
2018
2019 use Math::Complex ':pi';
2020 $third_of_circle = pi2 / 3;
2021
1515bec6
SP
2022=head2 Inf
2023
2024The floating point infinity can be exported as a subroutine Inf():
2025
2026 use Math::Complex qw(Inf sinh);
2027 my $AlsoInf = Inf() + 42;
2028 my $AnotherInf = sinh(1e42);
2029 print "$AlsoInf is $AnotherInf\n" if $AlsoInf == $AnotherInf;
2030
2031Note that the stringified form of infinity varies between platforms:
2032it can be for example any of
2033
2034 inf
2035 infinity
2036 INF
2037 1.#INF
2038
2039or it can be something else.
2040
7637cd07
SP
2041Also note that in some platforms trying to use the infinity in
2042arithmetic operations may result in Perl crashing because using
2043an infinity causes SIGFPE or its moral equivalent to be sent.
2044The way to ignore this is
2045
2046 local $SIG{FPE} = sub { };
2047
b42d0ec9 2048=head1 ERRORS DUE TO DIVISION BY ZERO OR LOGARITHM OF ZERO
5aabfad6
PP
2049
2050The division (/) and the following functions
2051
b42d0ec9 2052 log ln log10 logn
2820d885 2053 tan sec csc cot
b42d0ec9
JH
2054 atan asec acsc acot
2055 tanh sech csch coth
2056 atanh asech acsch acoth
5aabfad6
PP
2057
2058cannot be computed for all arguments because that would mean dividing
8c03c583
JH
2059by zero or taking logarithm of zero. These situations cause fatal
2060runtime errors looking like this
5aabfad6
PP
2061
2062 cot(0): Division by zero.
5cd24f17 2063 (Because in the definition of cot(0), the divisor sin(0) is 0)
5aabfad6
PP
2064 Died at ...
2065
8c03c583
JH
2066or
2067
2068 atanh(-1): Logarithm of zero.
2069 Died at...
2070
2071For the C<csc>, C<cot>, C<asec>, C<acsc>, C<acot>, C<csch>, C<coth>,
d1be9408 2072C<asech>, C<acsch>, the argument cannot be C<0> (zero). For the
b42d0ec9
JH
2073logarithmic functions and the C<atanh>, C<acoth>, the argument cannot
2074be C<1> (one). For the C<atanh>, C<acoth>, the argument cannot be
2075C<-1> (minus one). For the C<atan>, C<acot>, the argument cannot be
2076C<i> (the imaginary unit). For the C<atan>, C<acoth>, the argument
2077cannot be C<-i> (the negative imaginary unit). For the C<tan>,
2078C<sec>, C<tanh>, the argument cannot be I<pi/2 + k * pi>, where I<k>
bf5f1b4c
JH
2079is any integer. atan2(0, 0) is undefined, and if the complex arguments
2080are used for atan2(), a division by zero will happen if z1**2+z2**2 == 0.
b42d0ec9
JH
2081
2082Note that because we are operating on approximations of real numbers,
2083these errors can happen when merely `too close' to the singularities
40b904b7 2084listed above.
b42d0ec9
JH
2085
2086=head1 ERRORS DUE TO INDIGESTIBLE ARGUMENTS
2087
2088The C<make> and C<emake> accept both real and complex arguments.
2089When they cannot recognize the arguments they will die with error
2090messages like the following
2091
2092 Math::Complex::make: Cannot take real part of ...
2093 Math::Complex::make: Cannot take real part of ...
2094 Math::Complex::emake: Cannot take rho of ...
2095 Math::Complex::emake: Cannot take theta of ...
5cd24f17 2096
a5f75d66
AD
2097=head1 BUGS
2098
5cd24f17 2099Saying C<use Math::Complex;> exports many mathematical routines in the
bf5f1b4c 2100caller environment and even overrides some (C<sqrt>, C<log>, C<atan2>).
fb73857a 2101This is construed as a feature by the Authors, actually... ;-)
a5f75d66 2102
66730be0
RM
2103All routines expect to be given real or complex numbers. Don't attempt to
2104use BigFloat, since Perl has currently no rule to disambiguate a '+'
2105operation (for instance) between two overloaded entities.
a5f75d66 2106
d09ae4e6
JH
2107In Cray UNICOS there is some strange numerical instability that results
2108in root(), cos(), sin(), cosh(), sinh(), losing accuracy fast. Beware.
2109The bug may be in UNICOS math libs, in UNICOS C compiler, in Math::Complex.
2110Whatever it is, it does not manifest itself anywhere else where Perl runs.
2111
7637cd07
SP
2112=head1 SEE ALSO
2113
2114L<Math::Trig>
2115
0c721ce2 2116=head1 AUTHORS
a5f75d66 2117
3194a40a
Z
2118Daniel S. Lewart <F<lewart!at!uiuc.edu>>,
2119Jarkko Hietaniemi <F<jhi!at!iki.fi>>,
2120Raphael Manfredi <F<Raphael_Manfredi!at!pobox.com>>,
2121Zefram <zefram@fysh.org>
fb73857a 2122
1515bec6
SP
2123=head1 LICENSE
2124
2125This library is free software; you can redistribute it and/or modify
2126it under the same terms as Perl itself.
2127
5cd24f17
PP
2128=cut
2129
b42d0ec9
JH
21301;
2131
5cd24f17 2132# eof