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