This is a live mirror of the Perl 5 development currently hosted at https://github.com/perl/perl5
[inseparable changes from match from perl-5.003_97 to perl-5.003_97a]
[perl5.git] / lib / Math / Complex.pm
CommitLineData
66730be0
RM
1# $RCSFile$
2#
3# Complex numbers and associated mathematical functions
0c721ce2
JH
4# -- Raphael Manfredi, September 1996
5# -- Jarkko Hietaniemi, March 1997
a0d0e21e
LW
6
7require Exporter;
5aabfad6 8package Math::Complex;
a0d0e21e 9
0c721ce2
JH
10use strict;
11
5aabfad6
PP
12use vars qw($VERSION @ISA
13 @EXPORT %EXPORT_TAGS
14 $package $display
0c721ce2
JH
15 $pi $i $ilog10 $logn %logn);
16
5aabfad6
PP
17@ISA = qw(Exporter);
18
19$VERSION = 1.01;
20
21my @trig = qw(
22 pi
23 tan
24 csc cosec sec cot cotan
25 asin acos atan
26 acsc acosec asec acot acotan
27 sinh cosh tanh
28 csch cosech sech coth cotanh
29 asinh acosh atanh
30 acsch acosech asech acoth acotanh
31 );
32
33@EXPORT = (qw(
34 i Re Im arg
35 sqrt exp log ln
36 log10 logn cbrt root
37 cplx cplxe
38 ),
39 @trig);
40
41%EXPORT_TAGS = (
42 'trig' => [@trig],
66730be0 43);
a0d0e21e 44
a5f75d66 45use overload
0c721ce2
JH
46 '+' => \&plus,
47 '-' => \&minus,
48 '*' => \&multiply,
49 '/' => \&divide,
66730be0
RM
50 '**' => \&power,
51 '<=>' => \&spaceship,
52 'neg' => \&negate,
0c721ce2 53 '~' => \&conjugate,
66730be0
RM
54 'abs' => \&abs,
55 'sqrt' => \&sqrt,
56 'exp' => \&exp,
57 'log' => \&log,
58 'sin' => \&sin,
59 'cos' => \&cos,
0c721ce2 60 'tan' => \&tan,
66730be0
RM
61 'atan2' => \&atan2,
62 qw("" stringify);
63
64#
65# Package globals
66#
67
68$package = 'Math::Complex'; # Package name
69$display = 'cartesian'; # Default display format
70
71#
72# Object attributes (internal):
73# cartesian [real, imaginary] -- cartesian form
74# polar [rho, theta] -- polar form
75# c_dirty cartesian form not up-to-date
76# p_dirty polar form not up-to-date
77# display display format (package's global when not set)
78#
79
80#
81# ->make
82#
83# Create a new complex number (cartesian form)
84#
85sub make {
86 my $self = bless {}, shift;
87 my ($re, $im) = @_;
40da2db3 88 $self->{'cartesian'} = [$re, $im];
66730be0
RM
89 $self->{c_dirty} = 0;
90 $self->{p_dirty} = 1;
91 return $self;
92}
93
94#
95# ->emake
96#
97# Create a new complex number (exponential form)
98#
99sub emake {
100 my $self = bless {}, shift;
101 my ($rho, $theta) = @_;
102 $theta += pi() if $rho < 0;
40da2db3 103 $self->{'polar'} = [abs($rho), $theta];
66730be0
RM
104 $self->{p_dirty} = 0;
105 $self->{c_dirty} = 1;
106 return $self;
107}
108
109sub new { &make } # For backward compatibility only.
110
111#
112# cplx
113#
114# Creates a complex number from a (re, im) tuple.
115# This avoids the burden of writing Math::Complex->make(re, im).
116#
117sub cplx {
118 my ($re, $im) = @_;
0c721ce2 119 return $package->make($re, defined $im ? $im : 0);
66730be0
RM
120}
121
122#
123# cplxe
124#
125# Creates a complex number from a (rho, theta) tuple.
126# This avoids the burden of writing Math::Complex->emake(rho, theta).
127#
128sub cplxe {
129 my ($rho, $theta) = @_;
0c721ce2 130 return $package->emake($rho, defined $theta ? $theta : 0);
66730be0
RM
131}
132
133#
134# pi
135#
136# The number defined as 2 * pi = 360 degrees
137#
138sub pi () {
139 $pi = 4 * atan2(1, 1) unless $pi;
140 return $pi;
141}
142
143#
144# i
145#
146# The number defined as i*i = -1;
147#
148sub i () {
149 $i = bless {} unless $i; # There can be only one i
40da2db3
JH
150 $i->{'cartesian'} = [0, 1];
151 $i->{'polar'} = [1, pi/2];
66730be0
RM
152 $i->{c_dirty} = 0;
153 $i->{p_dirty} = 0;
154 return $i;
155}
156
157#
158# Attribute access/set routines
159#
160
0c721ce2
JH
161sub cartesian {$_[0]->{c_dirty} ?
162 $_[0]->update_cartesian : $_[0]->{'cartesian'}}
163sub polar {$_[0]->{p_dirty} ?
164 $_[0]->update_polar : $_[0]->{'polar'}}
66730be0 165
40da2db3
JH
166sub set_cartesian { $_[0]->{p_dirty}++; $_[0]->{'cartesian'} = $_[1] }
167sub set_polar { $_[0]->{c_dirty}++; $_[0]->{'polar'} = $_[1] }
66730be0
RM
168
169#
170# ->update_cartesian
171#
172# Recompute and return the cartesian form, given accurate polar form.
173#
174sub update_cartesian {
175 my $self = shift;
40da2db3 176 my ($r, $t) = @{$self->{'polar'}};
66730be0 177 $self->{c_dirty} = 0;
40da2db3 178 return $self->{'cartesian'} = [$r * cos $t, $r * sin $t];
66730be0
RM
179}
180
181#
182#
183# ->update_polar
184#
185# Recompute and return the polar form, given accurate cartesian form.
186#
187sub update_polar {
188 my $self = shift;
40da2db3 189 my ($x, $y) = @{$self->{'cartesian'}};
66730be0 190 $self->{p_dirty} = 0;
40da2db3
JH
191 return $self->{'polar'} = [0, 0] if $x == 0 && $y == 0;
192 return $self->{'polar'} = [sqrt($x*$x + $y*$y), atan2($y, $x)];
66730be0
RM
193}
194
195#
196# (plus)
197#
198# Computes z1+z2.
199#
200sub plus {
201 my ($z1, $z2, $regular) = @_;
0c721ce2 202 $z2 = cplx($z2, 0) unless ref $z2;
66730be0 203 my ($re1, $im1) = @{$z1->cartesian};
0c721ce2 204 my ($re2, $im2) = @{$z2->cartesian};
66730be0
RM
205 unless (defined $regular) {
206 $z1->set_cartesian([$re1 + $re2, $im1 + $im2]);
207 return $z1;
208 }
209 return (ref $z1)->make($re1 + $re2, $im1 + $im2);
210}
211
212#
213# (minus)
214#
215# Computes z1-z2.
216#
217sub minus {
218 my ($z1, $z2, $inverted) = @_;
0c721ce2 219 $z2 = cplx($z2, 0) unless ref $z2;
66730be0 220 my ($re1, $im1) = @{$z1->cartesian};
0c721ce2 221 my ($re2, $im2) = @{$z2->cartesian};
66730be0
RM
222 unless (defined $inverted) {
223 $z1->set_cartesian([$re1 - $re2, $im1 - $im2]);
224 return $z1;
225 }
226 return $inverted ?
227 (ref $z1)->make($re2 - $re1, $im2 - $im1) :
228 (ref $z1)->make($re1 - $re2, $im1 - $im2);
229}
230
231#
232# (multiply)
233#
234# Computes z1*z2.
235#
236sub multiply {
237 my ($z1, $z2, $regular) = @_;
238 my ($r1, $t1) = @{$z1->polar};
0c721ce2
JH
239 my ($r2, $t2) = ref $z2 ?
240 @{$z2->polar} : (abs($z2), $z2 >= 0 ? 0 : pi);
66730be0
RM
241 unless (defined $regular) {
242 $z1->set_polar([$r1 * $r2, $t1 + $t2]);
243 return $z1;
244 }
245 return (ref $z1)->emake($r1 * $r2, $t1 + $t2);
246}
247
248#
0c721ce2
JH
249# divbyzero
250#
251# Die on division by zero.
252#
253sub divbyzero {
5aabfad6 254 warn "$_[0]: Division by zero.\n";
0c721ce2
JH
255 warn "(Because in the definition of $_[0], $_[1] is 0)\n"
256 if (defined $_[1]);
257 my @up = caller(1);
258 my $dmess = "Died at $up[1] line $up[2].\n";
259 die $dmess;
260}
261
262#
66730be0
RM
263# (divide)
264#
265# Computes z1/z2.
266#
267sub divide {
268 my ($z1, $z2, $inverted) = @_;
269 my ($r1, $t1) = @{$z1->polar};
0c721ce2
JH
270 my ($r2, $t2) = ref $z2 ?
271 @{$z2->polar} : (abs($z2), $z2 >= 0 ? 0 : pi);
66730be0 272 unless (defined $inverted) {
0c721ce2 273 divbyzero "$z1/0" if ($r2 == 0);
66730be0
RM
274 $z1->set_polar([$r1 / $r2, $t1 - $t2]);
275 return $z1;
276 }
0c721ce2
JH
277 if ($inverted) {
278 divbyzero "$z2/0" if ($r1 == 0);
279 return (ref $z1)->emake($r2 / $r1, $t2 - $t1);
280 } else {
281 divbyzero "$z1/0" if ($r2 == 0);
282 return (ref $z1)->emake($r1 / $r2, $t1 - $t2);
283 }
66730be0
RM
284}
285
286#
287# (power)
288#
289# Computes z1**z2 = exp(z2 * log z1)).
290#
291sub power {
292 my ($z1, $z2, $inverted) = @_;
293 return exp($z1 * log $z2) if defined $inverted && $inverted;
294 return exp($z2 * log $z1);
295}
296
297#
298# (spaceship)
299#
300# Computes z1 <=> z2.
301# Sorts on the real part first, then on the imaginary part. Thus 2-4i > 3+8i.
302#
303sub spaceship {
304 my ($z1, $z2, $inverted) = @_;
0c721ce2 305 $z2 = cplx($z2, 0) unless ref $z2;
66730be0 306 my ($re1, $im1) = @{$z1->cartesian};
0c721ce2 307 my ($re2, $im2) = @{$z2->cartesian};
66730be0
RM
308 my $sgn = $inverted ? -1 : 1;
309 return $sgn * ($re1 <=> $re2) if $re1 != $re2;
310 return $sgn * ($im1 <=> $im2);
311}
312
313#
314# (negate)
315#
316# Computes -z.
317#
318sub negate {
319 my ($z) = @_;
320 if ($z->{c_dirty}) {
321 my ($r, $t) = @{$z->polar};
322 return (ref $z)->emake($r, pi + $t);
323 }
324 my ($re, $im) = @{$z->cartesian};
325 return (ref $z)->make(-$re, -$im);
326}
327
328#
329# (conjugate)
330#
331# Compute complex's conjugate.
332#
333sub conjugate {
334 my ($z) = @_;
335 if ($z->{c_dirty}) {
336 my ($r, $t) = @{$z->polar};
337 return (ref $z)->emake($r, -$t);
338 }
339 my ($re, $im) = @{$z->cartesian};
340 return (ref $z)->make($re, -$im);
341}
342
343#
344# (abs)
345#
346# Compute complex's norm (rho).
347#
348sub abs {
349 my ($z) = @_;
0c721ce2 350 return abs($z) unless ref $z;
66730be0
RM
351 my ($r, $t) = @{$z->polar};
352 return abs($r);
353}
354
355#
356# arg
357#
358# Compute complex's argument (theta).
359#
360sub arg {
361 my ($z) = @_;
0c721ce2 362 return ($z < 0 ? pi : 0) unless ref $z;
66730be0
RM
363 my ($r, $t) = @{$z->polar};
364 return $t;
365}
366
367#
368# (sqrt)
369#
0c721ce2 370# Compute sqrt(z).
66730be0
RM
371#
372sub sqrt {
373 my ($z) = @_;
0c721ce2 374 $z = cplx($z, 0) unless ref $z;
66730be0
RM
375 my ($r, $t) = @{$z->polar};
376 return (ref $z)->emake(sqrt($r), $t/2);
377}
378
379#
380# cbrt
381#
0c721ce2 382# Compute cbrt(z) (cubic root).
66730be0
RM
383#
384sub cbrt {
385 my ($z) = @_;
0c721ce2 386 return cplx($z, 0) ** (1/3) unless ref $z;
66730be0
RM
387 my ($r, $t) = @{$z->polar};
388 return (ref $z)->emake($r**(1/3), $t/3);
389}
390
391#
392# root
393#
394# Computes all nth root for z, returning an array whose size is n.
395# `n' must be a positive integer.
396#
397# The roots are given by (for k = 0..n-1):
398#
399# z^(1/n) = r^(1/n) (cos ((t+2 k pi)/n) + i sin ((t+2 k pi)/n))
400#
401sub root {
402 my ($z, $n) = @_;
403 $n = int($n + 0.5);
404 return undef unless $n > 0;
405 my ($r, $t) = ref $z ? @{$z->polar} : (abs($z), $z >= 0 ? 0 : pi);
406 my @root;
407 my $k;
408 my $theta_inc = 2 * pi / $n;
409 my $rho = $r ** (1/$n);
410 my $theta;
411 my $complex = ref($z) || $package;
412 for ($k = 0, $theta = $t / $n; $k < $n; $k++, $theta += $theta_inc) {
413 push(@root, $complex->emake($rho, $theta));
a0d0e21e 414 }
66730be0 415 return @root;
a0d0e21e
LW
416}
417
66730be0
RM
418#
419# Re
420#
421# Return Re(z).
422#
a0d0e21e 423sub Re {
66730be0
RM
424 my ($z) = @_;
425 return $z unless ref $z;
426 my ($re, $im) = @{$z->cartesian};
427 return $re;
a0d0e21e
LW
428}
429
66730be0
RM
430#
431# Im
432#
433# Return Im(z).
434#
a0d0e21e 435sub Im {
66730be0
RM
436 my ($z) = @_;
437 return 0 unless ref $z;
438 my ($re, $im) = @{$z->cartesian};
439 return $im;
a0d0e21e
LW
440}
441
66730be0
RM
442#
443# (exp)
444#
445# Computes exp(z).
446#
447sub exp {
448 my ($z) = @_;
0c721ce2 449 $z = cplx($z, 0) unless ref $z;
66730be0
RM
450 my ($x, $y) = @{$z->cartesian};
451 return (ref $z)->emake(exp($x), $y);
452}
453
454#
455# (log)
456#
457# Compute log(z).
458#
459sub log {
460 my ($z) = @_;
0c721ce2 461 $z = cplx($z, 0) unless ref $z;
66730be0 462 my ($r, $t) = @{$z->polar};
0c721ce2
JH
463 my ($x, $y) = @{$z->cartesian};
464 $t -= 2 * pi if ($t > pi() and $x < 0);
465 $t += 2 * pi if ($t < -pi() and $x < 0);
66730be0
RM
466 return (ref $z)->make(log($r), $t);
467}
468
469#
0c721ce2
JH
470# ln
471#
472# Alias for log().
473#
474sub ln { Math::Complex::log(@_) }
475
476#
66730be0
RM
477# log10
478#
479# Compute log10(z).
480#
481sub log10 {
482 my ($z) = @_;
0c721ce2
JH
483 my $ilog10 = 1 / log(10) unless defined $ilog10;
484 return log(cplx($z, 0)) * $ilog10 unless ref $z;
66730be0 485 my ($r, $t) = @{$z->polar};
0c721ce2 486 return (ref $z)->make(log($r) * $ilog10, $t * $ilog10);
66730be0
RM
487}
488
489#
490# logn
491#
492# Compute logn(z,n) = log(z) / log(n)
493#
494sub logn {
495 my ($z, $n) = @_;
0c721ce2 496 $z = cplx($z, 0) unless ref $z;
66730be0
RM
497 my $logn = $logn{$n};
498 $logn = $logn{$n} = log($n) unless defined $logn; # Cache log(n)
0c721ce2 499 return log($z) / $logn;
66730be0
RM
500}
501
502#
503# (cos)
504#
505# Compute cos(z) = (exp(iz) + exp(-iz))/2.
506#
507sub cos {
508 my ($z) = @_;
509 my ($x, $y) = @{$z->cartesian};
510 my $ey = exp($y);
511 my $ey_1 = 1 / $ey;
0c721ce2
JH
512 return (ref $z)->make(cos($x) * ($ey + $ey_1)/2,
513 sin($x) * ($ey_1 - $ey)/2);
66730be0
RM
514}
515
516#
517# (sin)
518#
519# Compute sin(z) = (exp(iz) - exp(-iz))/2.
520#
521sub sin {
522 my ($z) = @_;
523 my ($x, $y) = @{$z->cartesian};
524 my $ey = exp($y);
525 my $ey_1 = 1 / $ey;
0c721ce2
JH
526 return (ref $z)->make(sin($x) * ($ey + $ey_1)/2,
527 cos($x) * ($ey - $ey_1)/2);
66730be0
RM
528}
529
530#
531# tan
532#
533# Compute tan(z) = sin(z) / cos(z).
534#
535sub tan {
536 my ($z) = @_;
0c721ce2
JH
537 my $cz = cos($z);
538 divbyzero "tan($z)", "cos($z)" if ($cz == 0);
539 return sin($z) / $cz;
66730be0
RM
540}
541
542#
0c721ce2
JH
543# sec
544#
545# Computes the secant sec(z) = 1 / cos(z).
546#
547sub sec {
548 my ($z) = @_;
549 my $cz = cos($z);
550 divbyzero "sec($z)", "cos($z)" if ($cz == 0);
551 return 1 / $cz;
552}
553
554#
555# csc
556#
557# Computes the cosecant csc(z) = 1 / sin(z).
558#
559sub csc {
560 my ($z) = @_;
561 my $sz = sin($z);
562 divbyzero "csc($z)", "sin($z)" if ($sz == 0);
563 return 1 / $sz;
564}
565
66730be0 566#
0c721ce2 567# cosec
66730be0 568#
0c721ce2
JH
569# Alias for csc().
570#
571sub cosec { Math::Complex::csc(@_) }
572
573#
574# cot
575#
576# Computes cot(z) = 1 / tan(z).
577#
578sub cot {
66730be0 579 my ($z) = @_;
0c721ce2
JH
580 my $sz = sin($z);
581 divbyzero "cot($z)", "sin($z)" if ($sz == 0);
582 return cos($z) / $sz;
66730be0
RM
583}
584
585#
0c721ce2
JH
586# cotan
587#
588# Alias for cot().
589#
590sub cotan { Math::Complex::cot(@_) }
591
592#
66730be0
RM
593# acos
594#
595# Computes the arc cosine acos(z) = -i log(z + sqrt(z*z-1)).
596#
597sub acos {
598 my ($z) = @_;
0c721ce2
JH
599 $z = cplx($z, 0) unless ref $z;
600 return ~i * log($z + (Re($z) * Im($z) > 0 ? 1 : -1) * sqrt($z*$z - 1));
66730be0
RM
601}
602
603#
604# asin
605#
606# Computes the arc sine asin(z) = -i log(iz + sqrt(1-z*z)).
607#
608sub asin {
609 my ($z) = @_;
0c721ce2
JH
610 $z = cplx($z, 0) unless ref $z;
611 return ~i * log(i * $z + sqrt(1 - $z*$z));
66730be0
RM
612}
613
614#
615# atan
616#
0c721ce2 617# Computes the arc tangent atan(z) = i/2 log((i+z) / (i-z)).
66730be0
RM
618#
619sub atan {
620 my ($z) = @_;
0c721ce2
JH
621 divbyzero "atan($z)", "i - $z" if ($z == i);
622 return i/2*log((i + $z) / (i - $z));
a0d0e21e
LW
623}
624
66730be0 625#
0c721ce2
JH
626# asec
627#
628# Computes the arc secant asec(z) = acos(1 / z).
629#
630sub asec {
631 my ($z) = @_;
632 return acos(1 / $z);
633}
634
635#
636# acosec
637#
638# Computes the arc cosecant sec(z) = asin(1 / z).
639#
640sub acosec {
641 my ($z) = @_;
642 return asin(1 / $z);
643}
644
645#
646# acsc
66730be0 647#
0c721ce2
JH
648# Alias for acosec().
649#
650sub acsc { Math::Complex::acosec(@_) }
651
66730be0 652#
0c721ce2
JH
653# acot
654#
655# Computes the arc cotangent acot(z) = -i/2 log((i+z) / (z-i))
656#
657sub acot {
66730be0 658 my ($z) = @_;
0c721ce2 659 divbyzero "acot($z)", "$z - i" if ($z == i);
66730be0
RM
660 return i/-2 * log((i + $z) / ($z - i));
661}
662
663#
0c721ce2
JH
664# acotan
665#
666# Alias for acot().
667#
668sub acotan { Math::Complex::acot(@_) }
669
670#
66730be0
RM
671# cosh
672#
673# Computes the hyperbolic cosine cosh(z) = (exp(z) + exp(-z))/2.
674#
675sub cosh {
676 my ($z) = @_;
0c721ce2
JH
677 $z = cplx($z, 0) unless ref $z;
678 my ($x, $y) = @{$z->cartesian};
66730be0
RM
679 my $ex = exp($x);
680 my $ex_1 = 1 / $ex;
681 return ($ex + $ex_1)/2 unless ref $z;
0c721ce2
JH
682 return (ref $z)->make(cos($y) * ($ex + $ex_1)/2,
683 sin($y) * ($ex - $ex_1)/2);
66730be0
RM
684}
685
686#
687# sinh
688#
689# Computes the hyperbolic sine sinh(z) = (exp(z) - exp(-z))/2.
690#
691sub sinh {
692 my ($z) = @_;
0c721ce2
JH
693 $z = cplx($z, 0) unless ref $z;
694 my ($x, $y) = @{$z->cartesian};
66730be0
RM
695 my $ex = exp($x);
696 my $ex_1 = 1 / $ex;
697 return ($ex - $ex_1)/2 unless ref $z;
0c721ce2
JH
698 return (ref $z)->make(cos($y) * ($ex - $ex_1)/2,
699 sin($y) * ($ex + $ex_1)/2);
66730be0
RM
700}
701
702#
703# tanh
704#
705# Computes the hyperbolic tangent tanh(z) = sinh(z) / cosh(z).
706#
707sub tanh {
708 my ($z) = @_;
0c721ce2
JH
709 my $cz = cosh($z);
710 divbyzero "tanh($z)", "cosh($z)" if ($cz == 0);
711 return sinh($z) / $cz;
66730be0
RM
712}
713
714#
0c721ce2
JH
715# sech
716#
717# Computes the hyperbolic secant sech(z) = 1 / cosh(z).
718#
719sub sech {
720 my ($z) = @_;
721 my $cz = cosh($z);
722 divbyzero "sech($z)", "cosh($z)" if ($cz == 0);
723 return 1 / $cz;
724}
725
726#
727# csch
728#
729# Computes the hyperbolic cosecant csch(z) = 1 / sinh(z).
66730be0 730#
0c721ce2
JH
731sub csch {
732 my ($z) = @_;
733 my $sz = sinh($z);
734 divbyzero "csch($z)", "sinh($z)" if ($sz == 0);
735 return 1 / $sz;
736}
737
738#
739# cosech
740#
741# Alias for csch().
742#
743sub cosech { Math::Complex::csch(@_) }
744
66730be0 745#
0c721ce2
JH
746# coth
747#
748# Computes the hyperbolic cotangent coth(z) = cosh(z) / sinh(z).
749#
750sub coth {
66730be0 751 my ($z) = @_;
0c721ce2
JH
752 my $sz = sinh($z);
753 divbyzero "coth($z)", "sinh($z)" if ($sz == 0);
754 return cosh($z) / $sz;
66730be0
RM
755}
756
757#
0c721ce2
JH
758# cotanh
759#
760# Alias for coth().
761#
762sub cotanh { Math::Complex::coth(@_) }
763
764#
66730be0
RM
765# acosh
766#
767# Computes the arc hyperbolic cosine acosh(z) = log(z + sqrt(z*z-1)).
768#
769sub acosh {
770 my ($z) = @_;
0c721ce2
JH
771 $z = cplx($z, 0) unless ref $z; # asinh(-2)
772 return log($z + sqrt($z*$z - 1));
66730be0
RM
773}
774
775#
776# asinh
777#
778# Computes the arc hyperbolic sine asinh(z) = log(z + sqrt(z*z-1))
779#
780sub asinh {
781 my ($z) = @_;
0c721ce2
JH
782 $z = cplx($z, 0) unless ref $z; # asinh(-2)
783 return log($z + sqrt($z*$z + 1));
66730be0
RM
784}
785
786#
787# atanh
788#
789# Computes the arc hyperbolic tangent atanh(z) = 1/2 log((1+z) / (1-z)).
790#
791sub atanh {
792 my ($z) = @_;
0c721ce2
JH
793 $z = cplx($z, 0) unless ref $z; # atanh(-2)
794 divbyzero 'atanh(1)', "1 - $z" if ($z == 1);
66730be0 795 my $cz = (1 + $z) / (1 - $z);
66730be0
RM
796 return log($cz) / 2;
797}
798
799#
0c721ce2
JH
800# asech
801#
802# Computes the hyperbolic arc secant asech(z) = acosh(1 / z).
803#
804sub asech {
805 my ($z) = @_;
806 divbyzero 'asech(0)', $z if ($z == 0);
807 return acosh(1 / $z);
808}
809
810#
811# acsch
66730be0 812#
0c721ce2 813# Computes the hyperbolic arc cosecant acsch(z) = asinh(1 / z).
66730be0 814#
0c721ce2 815sub acsch {
66730be0 816 my ($z) = @_;
0c721ce2
JH
817 divbyzero 'acsch(0)', $z if ($z == 0);
818 return asinh(1 / $z);
819}
820
821#
822# acosech
823#
824# Alias for acosh().
825#
826sub acosech { Math::Complex::acsch(@_) }
827
828#
829# acoth
830#
831# Computes the arc hyperbolic cotangent acoth(z) = 1/2 log((1+z) / (z-1)).
832#
833sub acoth {
834 my ($z) = @_;
835 $z = cplx($z, 0) unless ref $z; # acoth(-2)
836 divbyzero 'acoth(1)', "$z - 1" if ($z == 1);
66730be0 837 my $cz = (1 + $z) / ($z - 1);
66730be0
RM
838 return log($cz) / 2;
839}
840
841#
0c721ce2
JH
842# acotanh
843#
844# Alias for acot().
845#
846sub acotanh { Math::Complex::acoth(@_) }
847
848#
66730be0
RM
849# (atan2)
850#
851# Compute atan(z1/z2).
852#
853sub atan2 {
854 my ($z1, $z2, $inverted) = @_;
855 my ($re1, $im1) = @{$z1->cartesian};
0c721ce2 856 my ($re2, $im2) = @{$z2->cartesian};
66730be0
RM
857 my $tan;
858 if (defined $inverted && $inverted) { # atan(z2/z1)
859 return pi * ($re2 > 0 ? 1 : -1) if $re1 == 0 && $im1 == 0;
860 $tan = $z2 / $z1;
861 } else {
862 return pi * ($re1 > 0 ? 1 : -1) if $re2 == 0 && $im2 == 0;
863 $tan = $z1 / $z2;
864 }
865 return atan($tan);
866}
867
868#
869# display_format
870# ->display_format
871#
872# Set (fetch if no argument) display format for all complex numbers that
873# don't happen to have overrriden it via ->display_format
874#
875# When called as a method, this actually sets the display format for
876# the current object.
877#
878# Valid object formats are 'c' and 'p' for cartesian and polar. The first
879# letter is used actually, so the type can be fully spelled out for clarity.
880#
881sub display_format {
882 my $self = shift;
883 my $format = undef;
884
885 if (ref $self) { # Called as a method
886 $format = shift;
0c721ce2 887 } else { # Regular procedure call
66730be0
RM
888 $format = $self;
889 undef $self;
890 }
891
892 if (defined $self) {
893 return defined $self->{display} ? $self->{display} : $display
894 unless defined $format;
895 return $self->{display} = $format;
896 }
897
898 return $display unless defined $format;
899 return $display = $format;
900}
901
902#
903# (stringify)
904#
905# Show nicely formatted complex number under its cartesian or polar form,
906# depending on the current display format:
907#
908# . If a specific display format has been recorded for this object, use it.
909# . Otherwise, use the generic current default for all complex numbers,
910# which is a package global variable.
911#
a0d0e21e 912sub stringify {
66730be0
RM
913 my ($z) = shift;
914 my $format;
915
916 $format = $display;
917 $format = $z->{display} if defined $z->{display};
918
919 return $z->stringify_polar if $format =~ /^p/i;
920 return $z->stringify_cartesian;
921}
922
923#
924# ->stringify_cartesian
925#
926# Stringify as a cartesian representation 'a+bi'.
927#
928sub stringify_cartesian {
929 my $z = shift;
930 my ($x, $y) = @{$z->cartesian};
931 my ($re, $im);
932
55497cff
PP
933 $x = int($x + ($x < 0 ? -1 : 1) * 1e-14)
934 if int(abs($x)) != int(abs($x) + 1e-14);
935 $y = int($y + ($y < 0 ? -1 : 1) * 1e-14)
936 if int(abs($y)) != int(abs($y) + 1e-14);
937
66730be0
RM
938 $re = "$x" if abs($x) >= 1e-14;
939 if ($y == 1) { $im = 'i' }
940 elsif ($y == -1) { $im = '-i' }
40da2db3 941 elsif (abs($y) >= 1e-14) { $im = $y . "i" }
66730be0 942
0c721ce2 943 my $str = '';
66730be0
RM
944 $str = $re if defined $re;
945 $str .= "+$im" if defined $im;
946 $str =~ s/\+-/-/;
947 $str =~ s/^\+//;
948 $str = '0' unless $str;
949
950 return $str;
951}
952
953#
954# ->stringify_polar
955#
956# Stringify as a polar representation '[r,t]'.
957#
958sub stringify_polar {
959 my $z = shift;
960 my ($r, $t) = @{$z->polar};
961 my $theta;
0c721ce2 962 my $eps = 1e-14;
66730be0 963
0c721ce2 964 return '[0,0]' if $r <= $eps;
a0d0e21e 965
66730be0
RM
966 my $tpi = 2 * pi;
967 my $nt = $t / $tpi;
968 $nt = ($nt - int($nt)) * $tpi;
969 $nt += $tpi if $nt < 0; # Range [0, 2pi]
a0d0e21e 970
0c721ce2
JH
971 if (abs($nt) <= $eps) { $theta = 0 }
972 elsif (abs(pi-$nt) <= $eps) { $theta = 'pi' }
66730be0 973
55497cff 974 if (defined $theta) {
0c721ce2
JH
975 $r = int($r + ($r < 0 ? -1 : 1) * $eps)
976 if int(abs($r)) != int(abs($r) + $eps);
977 $theta = int($theta + ($theta < 0 ? -1 : 1) * $eps)
978 if ($theta ne 'pi' and
979 int(abs($theta)) != int(abs($theta) + $eps));
55497cff
PP
980 return "\[$r,$theta\]";
981 }
66730be0
RM
982
983 #
984 # Okay, number is not a real. Try to identify pi/n and friends...
985 #
986
987 $nt -= $tpi if $nt > pi;
988 my ($n, $k, $kpi);
989
990 for ($k = 1, $kpi = pi; $k < 10; $k++, $kpi += pi) {
991 $n = int($kpi / $nt + ($nt > 0 ? 1 : -1) * 0.5);
0c721ce2
JH
992 if (abs($kpi/$n - $nt) <= $eps) {
993 $theta = ($nt < 0 ? '-':'').
994 ($k == 1 ? 'pi':"${k}pi").'/'.abs($n);
66730be0
RM
995 last;
996 }
997 }
998
999 $theta = $nt unless defined $theta;
1000
0c721ce2
JH
1001 $r = int($r + ($r < 0 ? -1 : 1) * $eps)
1002 if int(abs($r)) != int(abs($r) + $eps);
1003 $theta = int($theta + ($theta < 0 ? -1 : 1) * $eps)
1004 if ($theta !~ m(^-?\d*pi/\d+$) and
1005 int(abs($theta)) != int(abs($theta) + $eps));
55497cff 1006
66730be0 1007 return "\[$r,$theta\]";
a0d0e21e 1008}
a5f75d66
AD
1009
10101;
1011__END__
1012
1013=head1 NAME
1014
66730be0 1015Math::Complex - complex numbers and associated mathematical functions
a5f75d66
AD
1016
1017=head1 SYNOPSIS
1018
66730be0 1019 use Math::Complex;
5aabfad6 1020
66730be0
RM
1021 $z = Math::Complex->make(5, 6);
1022 $t = 4 - 3*i + $z;
1023 $j = cplxe(1, 2*pi/3);
a5f75d66
AD
1024
1025=head1 DESCRIPTION
1026
66730be0
RM
1027This package lets you create and manipulate complex numbers. By default,
1028I<Perl> limits itself to real numbers, but an extra C<use> statement brings
1029full complex support, along with a full set of mathematical functions
1030typically associated with and/or extended to complex numbers.
1031
1032If you wonder what complex numbers are, they were invented to be able to solve
1033the following equation:
1034
1035 x*x = -1
1036
1037and by definition, the solution is noted I<i> (engineers use I<j> instead since
1038I<i> usually denotes an intensity, but the name does not matter). The number
1039I<i> is a pure I<imaginary> number.
1040
1041The arithmetics with pure imaginary numbers works just like you would expect
1042it with real numbers... you just have to remember that
1043
1044 i*i = -1
1045
1046so you have:
1047
1048 5i + 7i = i * (5 + 7) = 12i
1049 4i - 3i = i * (4 - 3) = i
1050 4i * 2i = -8
1051 6i / 2i = 3
1052 1 / i = -i
1053
1054Complex numbers are numbers that have both a real part and an imaginary
1055part, and are usually noted:
1056
1057 a + bi
1058
1059where C<a> is the I<real> part and C<b> is the I<imaginary> part. The
1060arithmetic with complex numbers is straightforward. You have to
1061keep track of the real and the imaginary parts, but otherwise the
1062rules used for real numbers just apply:
1063
1064 (4 + 3i) + (5 - 2i) = (4 + 5) + i(3 - 2) = 9 + i
1065 (2 + i) * (4 - i) = 2*4 + 4i -2i -i*i = 8 + 2i + 1 = 9 + 2i
1066
1067A graphical representation of complex numbers is possible in a plane
1068(also called the I<complex plane>, but it's really a 2D plane).
1069The number
1070
1071 z = a + bi
1072
1073is the point whose coordinates are (a, b). Actually, it would
1074be the vector originating from (0, 0) to (a, b). It follows that the addition
1075of two complex numbers is a vectorial addition.
1076
1077Since there is a bijection between a point in the 2D plane and a complex
1078number (i.e. the mapping is unique and reciprocal), a complex number
1079can also be uniquely identified with polar coordinates:
1080
1081 [rho, theta]
1082
1083where C<rho> is the distance to the origin, and C<theta> the angle between
1084the vector and the I<x> axis. There is a notation for this using the
1085exponential form, which is:
1086
1087 rho * exp(i * theta)
1088
1089where I<i> is the famous imaginary number introduced above. Conversion
1090between this form and the cartesian form C<a + bi> is immediate:
1091
1092 a = rho * cos(theta)
1093 b = rho * sin(theta)
1094
1095which is also expressed by this formula:
1096
1097 z = rho * exp(i * theta) = rho * (cos theta + i * sin theta)
1098
1099In other words, it's the projection of the vector onto the I<x> and I<y>
1100axes. Mathematicians call I<rho> the I<norm> or I<modulus> and I<theta>
1101the I<argument> of the complex number. The I<norm> of C<z> will be
1102noted C<abs(z)>.
1103
1104The polar notation (also known as the trigonometric
1105representation) is much more handy for performing multiplications and
1106divisions of complex numbers, whilst the cartesian notation is better
1107suited for additions and substractions. Real numbers are on the I<x>
1108axis, and therefore I<theta> is zero.
1109
1110All the common operations that can be performed on a real number have
1111been defined to work on complex numbers as well, and are merely
1112I<extensions> of the operations defined on real numbers. This means
1113they keep their natural meaning when there is no imaginary part, provided
1114the number is within their definition set.
1115
1116For instance, the C<sqrt> routine which computes the square root of
1117its argument is only defined for positive real numbers and yields a
1118positive real number (it is an application from B<R+> to B<R+>).
1119If we allow it to return a complex number, then it can be extended to
1120negative real numbers to become an application from B<R> to B<C> (the
1121set of complex numbers):
1122
1123 sqrt(x) = x >= 0 ? sqrt(x) : sqrt(-x)*i
1124
1125It can also be extended to be an application from B<C> to B<C>,
1126whilst its restriction to B<R> behaves as defined above by using
1127the following definition:
1128
1129 sqrt(z = [r,t]) = sqrt(r) * exp(i * t/2)
1130
1131Indeed, a negative real number can be noted C<[x,pi]>
1132(the modulus I<x> is always positive, so C<[x,pi]> is really C<-x>, a
1133negative number)
1134and the above definition states that
1135
1136 sqrt([x,pi]) = sqrt(x) * exp(i*pi/2) = [sqrt(x),pi/2] = sqrt(x)*i
1137
1138which is exactly what we had defined for negative real numbers above.
a5f75d66 1139
66730be0
RM
1140All the common mathematical functions defined on real numbers that
1141are extended to complex numbers share that same property of working
1142I<as usual> when the imaginary part is zero (otherwise, it would not
1143be called an extension, would it?).
a5f75d66 1144
66730be0
RM
1145A I<new> operation possible on a complex number that is
1146the identity for real numbers is called the I<conjugate>, and is noted
1147with an horizontal bar above the number, or C<~z> here.
a5f75d66 1148
66730be0
RM
1149 z = a + bi
1150 ~z = a - bi
a5f75d66 1151
66730be0 1152Simple... Now look:
a5f75d66 1153
66730be0 1154 z * ~z = (a + bi) * (a - bi) = a*a + b*b
a5f75d66 1155
66730be0
RM
1156We saw that the norm of C<z> was noted C<abs(z)> and was defined as the
1157distance to the origin, also known as:
a5f75d66 1158
66730be0 1159 rho = abs(z) = sqrt(a*a + b*b)
a5f75d66 1160
66730be0
RM
1161so
1162
1163 z * ~z = abs(z) ** 2
1164
1165If z is a pure real number (i.e. C<b == 0>), then the above yields:
1166
1167 a * a = abs(a) ** 2
1168
1169which is true (C<abs> has the regular meaning for real number, i.e. stands
1170for the absolute value). This example explains why the norm of C<z> is
1171noted C<abs(z)>: it extends the C<abs> function to complex numbers, yet
1172is the regular C<abs> we know when the complex number actually has no
1173imaginary part... This justifies I<a posteriori> our use of the C<abs>
1174notation for the norm.
1175
1176=head1 OPERATIONS
1177
1178Given the following notations:
1179
1180 z1 = a + bi = r1 * exp(i * t1)
1181 z2 = c + di = r2 * exp(i * t2)
1182 z = <any complex or real number>
1183
1184the following (overloaded) operations are supported on complex numbers:
1185
1186 z1 + z2 = (a + c) + i(b + d)
1187 z1 - z2 = (a - c) + i(b - d)
1188 z1 * z2 = (r1 * r2) * exp(i * (t1 + t2))
1189 z1 / z2 = (r1 / r2) * exp(i * (t1 - t2))
1190 z1 ** z2 = exp(z2 * log z1)
1191 ~z1 = a - bi
1192 abs(z1) = r1 = sqrt(a*a + b*b)
1193 sqrt(z1) = sqrt(r1) * exp(i * t1/2)
1194 exp(z1) = exp(a) * exp(i * b)
1195 log(z1) = log(r1) + i*t1
1196 sin(z1) = 1/2i (exp(i * z1) - exp(-i * z1))
1197 cos(z1) = 1/2 (exp(i * z1) + exp(-i * z1))
1198 abs(z1) = r1
1199 atan2(z1, z2) = atan(z1/z2)
1200
1201The following extra operations are supported on both real and complex
1202numbers:
1203
1204 Re(z) = a
1205 Im(z) = b
1206 arg(z) = t
1207
1208 cbrt(z) = z ** (1/3)
1209 log10(z) = log(z) / log(10)
1210 logn(z, n) = log(z) / log(n)
1211
1212 tan(z) = sin(z) / cos(z)
0c721ce2 1213
5aabfad6
PP
1214 csc(z) = 1 / sin(z)
1215 sec(z) = 1 / cos(z)
0c721ce2 1216 cot(z) = 1 / tan(z)
66730be0
RM
1217
1218 asin(z) = -i * log(i*z + sqrt(1-z*z))
1219 acos(z) = -i * log(z + sqrt(z*z-1))
1220 atan(z) = i/2 * log((i+z) / (i-z))
0c721ce2 1221
5aabfad6
PP
1222 acsc(z) = asin(1 / z)
1223 asec(z) = acos(1 / z)
0c721ce2 1224 acot(z) = -i/2 * log((i+z) / (z-i))
66730be0
RM
1225
1226 sinh(z) = 1/2 (exp(z) - exp(-z))
1227 cosh(z) = 1/2 (exp(z) + exp(-z))
0c721ce2
JH
1228 tanh(z) = sinh(z) / cosh(z) = (exp(z) - exp(-z)) / (exp(z) + exp(-z))
1229
5aabfad6
PP
1230 csch(z) = 1 / sinh(z)
1231 sech(z) = 1 / cosh(z)
0c721ce2 1232 coth(z) = 1 / tanh(z)
66730be0
RM
1233
1234 asinh(z) = log(z + sqrt(z*z+1))
1235 acosh(z) = log(z + sqrt(z*z-1))
1236 atanh(z) = 1/2 * log((1+z) / (1-z))
66730be0 1237
5aabfad6
PP
1238 acsch(z) = asinh(1 / z)
1239 asech(z) = acosh(1 / z)
0c721ce2
JH
1240 acoth(z) = atanh(1 / z) = 1/2 * log((1+z) / (z-1))
1241
1242I<log>, I<csc>, I<cot>, I<acsc>, I<acot>, I<csch>, I<coth>,
1243I<acosech>, I<acotanh>, have aliases I<ln>, I<cosec>, I<cotan>,
1244I<acosec>, I<acotan>, I<cosech>, I<cotanh>, I<acosech>, I<acotanh>,
1245respectively.
1246
1247The I<root> function is available to compute all the I<n>
66730be0
RM
1248roots of some complex, where I<n> is a strictly positive integer.
1249There are exactly I<n> such roots, returned as a list. Getting the
1250number mathematicians call C<j> such that:
1251
1252 1 + j + j*j = 0;
1253
1254is a simple matter of writing:
1255
1256 $j = ((root(1, 3))[1];
1257
1258The I<k>th root for C<z = [r,t]> is given by:
1259
1260 (root(z, n))[k] = r**(1/n) * exp(i * (t + 2*k*pi)/n)
1261
0c721ce2
JH
1262The I<spaceship> comparison operator is also defined. In order to
1263ensure its restriction to real numbers is conform to what you would
1264expect, the comparison is run on the real part of the complex number
1265first, and imaginary parts are compared only when the real parts
1266match.
66730be0
RM
1267
1268=head1 CREATION
1269
1270To create a complex number, use either:
1271
1272 $z = Math::Complex->make(3, 4);
1273 $z = cplx(3, 4);
1274
1275if you know the cartesian form of the number, or
1276
1277 $z = 3 + 4*i;
1278
1279if you like. To create a number using the trigonometric form, use either:
1280
1281 $z = Math::Complex->emake(5, pi/3);
1282 $x = cplxe(5, pi/3);
1283
0c721ce2
JH
1284instead. The first argument is the modulus, the second is the angle
1285(in radians, the full circle is 2*pi). (Mnmemonic: C<e> is used as a
1286notation for complex numbers in the trigonometric form).
66730be0
RM
1287
1288It is possible to write:
1289
1290 $x = cplxe(-3, pi/4);
1291
1292but that will be silently converted into C<[3,-3pi/4]>, since the modulus
1293must be positive (it represents the distance to the origin in the complex
1294plane).
1295
1296=head1 STRINGIFICATION
1297
1298When printed, a complex number is usually shown under its cartesian
1299form I<a+bi>, but there are legitimate cases where the polar format
1300I<[r,t]> is more appropriate.
1301
1302By calling the routine C<Math::Complex::display_format> and supplying either
1303C<"polar"> or C<"cartesian">, you override the default display format,
1304which is C<"cartesian">. Not supplying any argument returns the current
1305setting.
1306
1307This default can be overridden on a per-number basis by calling the
1308C<display_format> method instead. As before, not supplying any argument
1309returns the current display format for this number. Otherwise whatever you
1310specify will be the new display format for I<this> particular number.
1311
1312For instance:
1313
1314 use Math::Complex;
1315
1316 Math::Complex::display_format('polar');
1317 $j = ((root(1, 3))[1];
1318 print "j = $j\n"; # Prints "j = [1,2pi/3]
1319 $j->display_format('cartesian');
1320 print "j = $j\n"; # Prints "j = -0.5+0.866025403784439i"
1321
1322The polar format attempts to emphasize arguments like I<k*pi/n>
1323(where I<n> is a positive integer and I<k> an integer within [-9,+9]).
1324
1325=head1 USAGE
1326
1327Thanks to overloading, the handling of arithmetics with complex numbers
1328is simple and almost transparent.
1329
1330Here are some examples:
1331
1332 use Math::Complex;
1333
1334 $j = cplxe(1, 2*pi/3); # $j ** 3 == 1
1335 print "j = $j, j**3 = ", $j ** 3, "\n";
1336 print "1 + j + j**2 = ", 1 + $j + $j**2, "\n";
1337
1338 $z = -16 + 0*i; # Force it to be a complex
1339 print "sqrt($z) = ", sqrt($z), "\n";
1340
1341 $k = exp(i * 2*pi/3);
1342 print "$j - $k = ", $j - $k, "\n";
a5f75d66 1343
5aabfad6
PP
1344=head1 CAVEATS
1345
1346The division (/) and the following functions
1347
1348 tan
1349 sec
1350 csc
1351 cot
1352 atan
1353 acot
1354 tanh
1355 sech
1356 csch
1357 coth
1358 atanh
1359 asech
1360 acsch
1361 acoth
1362
1363cannot be computed for all arguments because that would mean dividing
1364by zero. These situations cause fatal runtime errors looking like this
1365
1366 cot(0): Division by zero.
1367 (Because in the definition of cot(0), sin(0) is 0)
1368 Died at ...
1369
a5f75d66
AD
1370=head1 BUGS
1371
66730be0
RM
1372Saying C<use Math::Complex;> exports many mathematical routines in the caller
1373environment. This is construed as a feature by the Author, actually... ;-)
1374
1375The code is not optimized for speed, although we try to use the cartesian
1376form for addition-like operators and the trigonometric form for all
1377multiplication-like operators.
1378
1379The arg() routine does not ensure the angle is within the range [-pi,+pi]
1380(a side effect caused by multiplication and division using the trigonometric
1381representation).
a5f75d66 1382
66730be0
RM
1383All routines expect to be given real or complex numbers. Don't attempt to
1384use BigFloat, since Perl has currently no rule to disambiguate a '+'
1385operation (for instance) between two overloaded entities.
a5f75d66 1386
0c721ce2 1387=head1 AUTHORS
a5f75d66 1388
0c721ce2
JH
1389 Raphael Manfredi <F<Raphael_Manfredi@grenoble.hp.com>>
1390 Jarkko Hietaniemi <F<jhi@iki.fi>>