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