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