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