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