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