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