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