X-Git-Url: http://git.shadowcat.co.uk/gitweb/gitweb.cgi?a=blobdiff_plain;f=lib%2FMath%2FComplex.pm;h=5b69039afc07ccb1bdcd2e0c6081fb43444944c5;hb=f4c556ac9d141bf86702c68d95acad2db5ec6874;hp=5ec4a5661ecd7d5620635b9c1b42ec1a44c006a1;hpb=66730be00e20f21554ddafe7899e27bc47edc3cb;p=p5sagit%2Fp5-mst-13.2.git diff --git a/lib/Math/Complex.pm b/lib/Math/Complex.pm index 5ec4a56..5b69039 100644 --- a/lib/Math/Complex.pm +++ b/lib/Math/Complex.pm @@ -1,43 +1,73 @@ -# $RCSFile$ # # Complex numbers and associated mathematical functions -# -- Raphael Manfredi, Sept 1996 +# -- Raphael Manfredi Since Sep 1996 +# -- Jarkko Hietaniemi Since Mar 1997 +# -- Daniel S. Lewart Since Sep 1997 +# require Exporter; -package Math::Complex; @ISA = qw(Exporter); - -@EXPORT = qw( - pi i Re Im arg - log10 logn cbrt root - tan cotan asin acos atan acotan - sinh cosh tanh cotanh asinh acosh atanh acotanh - cplx cplxe +package Math::Complex; + +use strict; + +use vars qw($VERSION @ISA @EXPORT %EXPORT_TAGS); + +my ( $i, $ip2, %logn ); + +$VERSION = sprintf("%s", q$Id: Complex.pm,v 1.26 1998/11/01 00:00:00 dsl Exp $ =~ /(\d+\.\d+)/); + +@ISA = qw(Exporter); + +my @trig = qw( + pi + tan + csc cosec sec cot cotan + asin acos atan + acsc acosec asec acot acotan + sinh cosh tanh + csch cosech sech coth cotanh + asinh acosh atanh + acsch acosech asech acoth acotanh + ); + +@EXPORT = (qw( + i Re Im rho theta arg + sqrt log ln + log10 logn cbrt root + cplx cplxe + ), + @trig); + +%EXPORT_TAGS = ( + 'trig' => [@trig], ); use overload - '+' => \&plus, - '-' => \&minus, - '*' => \&multiply, - '/' => \÷, + '+' => \&plus, + '-' => \&minus, + '*' => \&multiply, + '/' => \÷, '**' => \&power, '<=>' => \&spaceship, 'neg' => \&negate, - '~' => \&conjugate, + '~' => \&conjugate, 'abs' => \&abs, 'sqrt' => \&sqrt, 'exp' => \&exp, 'log' => \&log, 'sin' => \&sin, 'cos' => \&cos, + 'tan' => \&tan, 'atan2' => \&atan2, qw("" stringify); # -# Package globals +# Package "privates" # -$package = 'Math::Complex'; # Package name -$display = 'cartesian'; # Default display format +my $package = 'Math::Complex'; # Package name +my $display = 'cartesian'; # Default display format +my $eps = 1e-14; # Epsilon # # Object attributes (internal): @@ -48,6 +78,12 @@ $display = 'cartesian'; # Default display format # display display format (package's global when not set) # +# Die on bad *make() arguments. + +sub _cannot_make { + die "@{[(caller(1))[3]]}: Cannot take $_[0] of $_[1].\n"; +} + # # ->make # @@ -56,9 +92,26 @@ $display = 'cartesian'; # Default display format sub make { my $self = bless {}, shift; my ($re, $im) = @_; - $self->{cartesian} = [$re, $im]; + my $rre = ref $re; + if ( $rre ) { + if ( $rre eq ref $self ) { + $re = Re($re); + } else { + _cannot_make("real part", $rre); + } + } + my $rim = ref $im; + if ( $rim ) { + if ( $rim eq ref $self ) { + $im = Im($im); + } else { + _cannot_make("imaginary part", $rim); + } + } + $self->{'cartesian'} = [ $re, $im ]; $self->{c_dirty} = 0; $self->{p_dirty} = 1; + $self->display_format('cartesian'); return $self; } @@ -70,10 +123,30 @@ sub make { sub emake { my $self = bless {}, shift; my ($rho, $theta) = @_; - $theta += pi() if $rho < 0; - $self->{polar} = [abs($rho), $theta]; + my $rrh = ref $rho; + if ( $rrh ) { + if ( $rrh eq ref $self ) { + $rho = rho($rho); + } else { + _cannot_make("rho", $rrh); + } + } + my $rth = ref $theta; + if ( $rth ) { + if ( $rth eq ref $self ) { + $theta = theta($theta); + } else { + _cannot_make("theta", $rth); + } + } + if ($rho < 0) { + $rho = -$rho; + $theta = ($theta <= 0) ? $theta + pi() : $theta - pi(); + } + $self->{'polar'} = [$rho, $theta]; $self->{p_dirty} = 0; $self->{c_dirty} = 1; + $self->display_format('polar'); return $self; } @@ -87,7 +160,7 @@ sub new { &make } # For backward compatibility only. # sub cplx { my ($re, $im) = @_; - return $package->make($re, $im); + return $package->make($re, defined $im ? $im : 0); } # @@ -98,18 +171,44 @@ sub cplx { # sub cplxe { my ($rho, $theta) = @_; - return $package->emake($rho, $theta); + return $package->emake($rho, defined $theta ? $theta : 0); } # # pi # -# The number defined as 2 * pi = 360 degrees +# The number defined as pi = 180 degrees # -sub pi () { - $pi = 4 * atan2(1, 1) unless $pi; - return $pi; -} +use constant pi => 4 * CORE::atan2(1, 1); + +# +# pit2 +# +# The full circle +# +use constant pit2 => 2 * pi; + +# +# pip2 +# +# The quarter circle +# +use constant pip2 => pi / 2; + +# +# deg1 +# +# One degree in radians, used in stringify_polar. +# + +use constant deg1 => pi / 180; + +# +# uplog10 +# +# Used in log10(). +# +use constant uplog10 => 1 / CORE::log(10); # # i @@ -117,9 +216,10 @@ sub pi () { # The number defined as i*i = -1; # sub i () { - $i = bless {} unless $i; # There can be only one i - $i->{cartesian} = [0, 1]; - $i->{polar} = [1, pi/2]; + return $i if ($i); + $i = bless {}; + $i->{'cartesian'} = [0, 1]; + $i->{'polar'} = [1, pip2]; $i->{c_dirty} = 0; $i->{p_dirty} = 0; return $i; @@ -129,11 +229,13 @@ sub i () { # Attribute access/set routines # -sub cartesian {$_[0]->{c_dirty} ? $_[0]->update_cartesian : $_[0]->{cartesian}} -sub polar {$_[0]->{p_dirty} ? $_[0]->update_polar : $_[0]->{polar}} +sub cartesian {$_[0]->{c_dirty} ? + $_[0]->update_cartesian : $_[0]->{'cartesian'}} +sub polar {$_[0]->{p_dirty} ? + $_[0]->update_polar : $_[0]->{'polar'}} -sub set_cartesian { $_[0]->{p_dirty}++; $_[0]->{cartesian} = $_[1] } -sub set_polar { $_[0]->{c_dirty}++; $_[0]->{polar} = $_[1] } +sub set_cartesian { $_[0]->{p_dirty}++; $_[0]->{'cartesian'} = $_[1] } +sub set_polar { $_[0]->{c_dirty}++; $_[0]->{'polar'} = $_[1] } # # ->update_cartesian @@ -142,9 +244,9 @@ sub set_polar { $_[0]->{c_dirty}++; $_[0]->{polar} = $_[1] } # sub update_cartesian { my $self = shift; - my ($r, $t) = @{$self->{polar}}; + my ($r, $t) = @{$self->{'polar'}}; $self->{c_dirty} = 0; - return $self->{cartesian} = [$r * cos $t, $r * sin $t]; + return $self->{'cartesian'} = [$r * CORE::cos($t), $r * CORE::sin($t)]; } # @@ -155,10 +257,10 @@ sub update_cartesian { # sub update_polar { my $self = shift; - my ($x, $y) = @{$self->{cartesian}}; + my ($x, $y) = @{$self->{'cartesian'}}; $self->{p_dirty} = 0; - return $self->{polar} = [0, 0] if $x == 0 && $y == 0; - return $self->{polar} = [sqrt($x*$x + $y*$y), atan2($y, $x)]; + return $self->{'polar'} = [0, 0] if $x == 0 && $y == 0; + return $self->{'polar'} = [CORE::sqrt($x*$x + $y*$y), CORE::atan2($y, $x)]; } # @@ -169,7 +271,8 @@ sub update_polar { sub plus { my ($z1, $z2, $regular) = @_; my ($re1, $im1) = @{$z1->cartesian}; - my ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2); + $z2 = cplx($z2) unless ref $z2; + my ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0); unless (defined $regular) { $z1->set_cartesian([$re1 + $re2, $im1 + $im2]); return $z1; @@ -185,7 +288,8 @@ sub plus { sub minus { my ($z1, $z2, $inverted) = @_; my ($re1, $im1) = @{$z1->cartesian}; - my ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2); + $z2 = cplx($z2) unless ref $z2; + my ($re2, $im2) = @{$z2->cartesian}; unless (defined $inverted) { $z1->set_cartesian([$re1 - $re2, $im1 - $im2]); return $z1; @@ -193,6 +297,7 @@ sub minus { return $inverted ? (ref $z1)->make($re2 - $re1, $im2 - $im1) : (ref $z1)->make($re1 - $re2, $im1 - $im2); + } # @@ -201,14 +306,49 @@ sub minus { # Computes z1*z2. # sub multiply { - my ($z1, $z2, $regular) = @_; - my ($r1, $t1) = @{$z1->polar}; - my ($r2, $t2) = ref $z2 ? @{$z2->polar} : (abs($z2), $z2 >= 0 ? 0 : pi); - unless (defined $regular) { - $z1->set_polar([$r1 * $r2, $t1 + $t2]); + my ($z1, $z2, $regular) = @_; + if ($z1->{p_dirty} == 0 and ref $z2 and $z2->{p_dirty} == 0) { + # if both polar better use polar to avoid rounding errors + my ($r1, $t1) = @{$z1->polar}; + my ($r2, $t2) = @{$z2->polar}; + my $t = $t1 + $t2; + if ($t > pi()) { $t -= pit2 } + elsif ($t <= -pi()) { $t += pit2 } + unless (defined $regular) { + $z1->set_polar([$r1 * $r2, $t]); return $z1; + } + return (ref $z1)->emake($r1 * $r2, $t); + } else { + my ($x1, $y1) = @{$z1->cartesian}; + if (ref $z2) { + my ($x2, $y2) = @{$z2->cartesian}; + return (ref $z1)->make($x1*$x2-$y1*$y2, $x1*$y2+$y1*$x2); + } else { + return (ref $z1)->make($x1*$z2, $y1*$z2); + } } - return (ref $z1)->emake($r1 * $r2, $t1 + $t2); +} + +# +# _divbyzero +# +# Die on division by zero. +# +sub _divbyzero { + my $mess = "$_[0]: Division by zero.\n"; + + if (defined $_[1]) { + $mess .= "(Because in the definition of $_[0], the divisor "; + $mess .= "$_[1] " unless ($_[1] eq '0'); + $mess .= "is 0)\n"; + } + + my @up = caller(1); + + $mess .= "Died at $up[1] line $up[2].\n"; + + die $mess; } # @@ -218,15 +358,46 @@ sub multiply { # sub divide { my ($z1, $z2, $inverted) = @_; - my ($r1, $t1) = @{$z1->polar}; - my ($r2, $t2) = ref $z2 ? @{$z2->polar} : (abs($z2), $z2 >= 0 ? 0 : pi); - unless (defined $inverted) { - $z1->set_polar([$r1 / $r2, $t1 - $t2]); - return $z1; + if ($z1->{p_dirty} == 0 and ref $z2 and $z2->{p_dirty} == 0) { + # if both polar better use polar to avoid rounding errors + my ($r1, $t1) = @{$z1->polar}; + my ($r2, $t2) = @{$z2->polar}; + my $t; + if ($inverted) { + _divbyzero "$z2/0" if ($r1 == 0); + $t = $t2 - $t1; + if ($t > pi()) { $t -= pit2 } + elsif ($t <= -pi()) { $t += pit2 } + return (ref $z1)->emake($r2 / $r1, $t); + } else { + _divbyzero "$z1/0" if ($r2 == 0); + $t = $t1 - $t2; + if ($t > pi()) { $t -= pit2 } + elsif ($t <= -pi()) { $t += pit2 } + return (ref $z1)->emake($r1 / $r2, $t); + } + } else { + my ($d, $x2, $y2); + if ($inverted) { + ($x2, $y2) = @{$z1->cartesian}; + $d = $x2*$x2 + $y2*$y2; + _divbyzero "$z2/0" if $d == 0; + return (ref $z1)->make(($x2*$z2)/$d, -($y2*$z2)/$d); + } else { + my ($x1, $y1) = @{$z1->cartesian}; + if (ref $z2) { + ($x2, $y2) = @{$z2->cartesian}; + $d = $x2*$x2 + $y2*$y2; + _divbyzero "$z1/0" if $d == 0; + my $u = ($x1*$x2 + $y1*$y2)/$d; + my $v = ($y1*$x2 - $x1*$y2)/$d; + return (ref $z1)->make($u, $v); + } else { + _divbyzero "$z1/0" if $z2 == 0; + return (ref $z1)->make($x1/$z2, $y1/$z2); + } + } } - return $inverted ? - (ref $z1)->emake($r2 / $r1, $t2 - $t1) : - (ref $z1)->emake($r1 / $r2, $t1 - $t2); } # @@ -236,20 +407,31 @@ sub divide { # sub power { my ($z1, $z2, $inverted) = @_; - return exp($z1 * log $z2) if defined $inverted && $inverted; - return exp($z2 * log $z1); + if ($inverted) { + return 1 if $z1 == 0 || $z2 == 1; + return 0 if $z2 == 0 && Re($z1) > 0; + } else { + return 1 if $z2 == 0 || $z1 == 1; + return 0 if $z1 == 0 && Re($z2) > 0; + } + my $w = $inverted ? CORE::exp($z1 * CORE::log($z2)) + : CORE::exp($z2 * CORE::log($z1)); + # If both arguments cartesian, return cartesian, else polar. + return $z1->{c_dirty} == 0 && + (not ref $z2 or $z2->{c_dirty} == 0) ? + cplx(@{$w->cartesian}) : $w; } # # (spaceship) # # Computes z1 <=> z2. -# Sorts on the real part first, then on the imaginary part. Thus 2-4i > 3+8i. +# Sorts on the real part first, then on the imaginary part. Thus 2-4i < 3+8i. # sub spaceship { my ($z1, $z2, $inverted) = @_; - my ($re1, $im1) = @{$z1->cartesian}; - my ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2); + my ($re1, $im1) = ref $z1 ? @{$z1->cartesian} : ($z1, 0); + my ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0); my $sgn = $inverted ? -1 : 1; return $sgn * ($re1 <=> $re2) if $re1 != $re2; return $sgn * ($im1 <=> $im2); @@ -264,7 +446,8 @@ sub negate { my ($z) = @_; if ($z->{c_dirty}) { my ($r, $t) = @{$z->polar}; - return (ref $z)->emake($r, pi + $t); + $t = ($t <= 0) ? $t + pi : $t - pi; + return (ref $z)->emake($r, $t); } my ($re, $im) = @{$z->cartesian}; return (ref $z)->make(-$re, -$im); @@ -288,47 +471,99 @@ sub conjugate { # # (abs) # -# Compute complex's norm (rho). +# Compute or set complex's norm (rho). # sub abs { - my ($z) = @_; - my ($r, $t) = @{$z->polar}; - return abs($r); + my ($z, $rho) = @_; + return $z unless ref $z; + if (defined $rho) { + $z->{'polar'} = [ $rho, ${$z->polar}[1] ]; + $z->{p_dirty} = 0; + $z->{c_dirty} = 1; + return $rho; + } else { + return ${$z->polar}[0]; + } +} + +sub _theta { + my $theta = $_[0]; + + if ($$theta > pi()) { $$theta -= pit2 } + elsif ($$theta <= -pi()) { $$theta += pit2 } } # # arg # -# Compute complex's argument (theta). +# Compute or set complex's argument (theta). # sub arg { - my ($z) = @_; - return 0 unless ref $z; - my ($r, $t) = @{$z->polar}; - return $t; + my ($z, $theta) = @_; + return $z unless ref $z; + if (defined $theta) { + _theta(\$theta); + $z->{'polar'} = [ ${$z->polar}[0], $theta ]; + $z->{p_dirty} = 0; + $z->{c_dirty} = 1; + } else { + $theta = ${$z->polar}[1]; + _theta(\$theta); + } + return $theta; } # # (sqrt) # -# Compute sqrt(z) (positive only). +# Compute sqrt(z). +# +# It is quite tempting to use wantarray here so that in list context +# sqrt() would return the two solutions. This, however, would +# break things like +# +# print "sqrt(z) = ", sqrt($z), "\n"; +# +# The two values would be printed side by side without no intervening +# whitespace, quite confusing. +# Therefore if you want the two solutions use the root(). # sub sqrt { my ($z) = @_; + my ($re, $im) = ref $z ? @{$z->cartesian} : ($z, 0); + return $re < 0 ? cplx(0, CORE::sqrt(-$re)) : CORE::sqrt($re) if $im == 0; my ($r, $t) = @{$z->polar}; - return (ref $z)->emake(sqrt($r), $t/2); + return (ref $z)->emake(CORE::sqrt($r), $t/2); } # # cbrt # -# Compute cbrt(z) (cubic root, primary only). +# Compute cbrt(z) (cubic root). +# +# Why are we not returning three values? The same answer as for sqrt(). # sub cbrt { my ($z) = @_; - return $z ** (1/3) unless ref $z; + return $z < 0 ? -CORE::exp(CORE::log(-$z)/3) : ($z > 0 ? CORE::exp(CORE::log($z)/3): 0) + unless ref $z; my ($r, $t) = @{$z->polar}; - return (ref $z)->emake($r**(1/3), $t/3); + return (ref $z)->emake(CORE::exp(CORE::log($r)/3), $t/3); +} + +# +# _rootbad +# +# Die on bad root. +# +sub _rootbad { + my $mess = "Root $_[0] not defined, root must be positive integer.\n"; + + my @up = caller(1); + + $mess .= "Died at $up[1] line $up[2].\n"; + + die $mess; } # @@ -343,17 +578,18 @@ sub cbrt { # sub root { my ($z, $n) = @_; - $n = int($n + 0.5); - return undef unless $n > 0; - my ($r, $t) = ref $z ? @{$z->polar} : (abs($z), $z >= 0 ? 0 : pi); + _rootbad($n) if ($n < 1 or int($n) != $n); + my ($r, $t) = ref $z ? @{$z->polar} : (CORE::abs($z), $z >= 0 ? 0 : pi); my @root; my $k; - my $theta_inc = 2 * pi / $n; + my $theta_inc = pit2 / $n; my $rho = $r ** (1/$n); my $theta; - my $complex = ref($z) || $package; + my $cartesian = ref $z && $z->{c_dirty} == 0; for ($k = 0, $theta = $t / $n; $k < $n; $k++, $theta += $theta_inc) { - push(@root, $complex->emake($rho, $theta)); + my $w = cplxe($rho, $theta); + # Yes, $cartesian is loop invariant. + push @root, $cartesian ? cplx(@{$w->cartesian}) : $w; } return @root; } @@ -361,25 +597,53 @@ sub root { # # Re # -# Return Re(z). +# Return or set Re(z). # sub Re { - my ($z) = @_; + my ($z, $Re) = @_; return $z unless ref $z; - my ($re, $im) = @{$z->cartesian}; - return $re; + if (defined $Re) { + $z->{'cartesian'} = [ $Re, ${$z->cartesian}[1] ]; + $z->{c_dirty} = 0; + $z->{p_dirty} = 1; + } else { + return ${$z->cartesian}[0]; + } } # # Im # -# Return Im(z). +# Return or set Im(z). # sub Im { - my ($z) = @_; - return 0 unless ref $z; - my ($re, $im) = @{$z->cartesian}; - return $im; + my ($z, $Im) = @_; + return $z unless ref $z; + if (defined $Im) { + $z->{'cartesian'} = [ ${$z->cartesian}[0], $Im ]; + $z->{c_dirty} = 0; + $z->{p_dirty} = 1; + } else { + return ${$z->cartesian}[1]; + } +} + +# +# rho +# +# Return or set rho(w). +# +sub rho { + Math::Complex::abs(@_); +} + +# +# theta +# +# Return or set theta(w). +# +sub theta { + Math::Complex::arg(@_); } # @@ -390,7 +654,28 @@ sub Im { sub exp { my ($z) = @_; my ($x, $y) = @{$z->cartesian}; - return (ref $z)->emake(exp($x), $y); + return (ref $z)->emake(CORE::exp($x), $y); +} + +# +# _logofzero +# +# Die on logarithm of zero. +# +sub _logofzero { + my $mess = "$_[0]: Logarithm of zero.\n"; + + if (defined $_[1]) { + $mess .= "(Because in the definition of $_[0], the argument "; + $mess .= "$_[1] " unless ($_[1] eq '0'); + $mess .= "is 0)\n"; + } + + my @up = caller(1); + + $mess .= "Died at $up[1] line $up[2].\n"; + + die $mess; } # @@ -400,21 +685,32 @@ sub exp { # sub log { my ($z) = @_; + unless (ref $z) { + _logofzero("log") if $z == 0; + return $z > 0 ? CORE::log($z) : cplx(CORE::log(-$z), pi); + } my ($r, $t) = @{$z->polar}; - return (ref $z)->make(log($r), $t); + _logofzero("log") if $r == 0; + if ($t > pi()) { $t -= pit2 } + elsif ($t <= -pi()) { $t += pit2 } + return (ref $z)->make(CORE::log($r), $t); } # +# ln +# +# Alias for log(). +# +sub ln { Math::Complex::log(@_) } + +# # log10 # # Compute log10(z). # + sub log10 { - my ($z) = @_; - $log10 = log(10) unless defined $log10; - return log($z) / $log10 unless ref $z; - my ($r, $t) = @{$z->polar}; - return (ref $z)->make(log($r) / $log10, $t / $log10); + return Math::Complex::log($_[0]) * uplog10; } # @@ -424,9 +720,10 @@ sub log10 { # sub logn { my ($z, $n) = @_; + $z = cplx($z, 0) unless ref $z; my $logn = $logn{$n}; - $logn = $logn{$n} = log($n) unless defined $logn; # Cache log(n) - return log($z) / log($n); + $logn = $logn{$n} = CORE::log($n) unless defined $logn; # Cache log(n) + return CORE::log($z) / $logn; } # @@ -437,9 +734,10 @@ sub logn { sub cos { my ($z) = @_; my ($x, $y) = @{$z->cartesian}; - my $ey = exp($y); + my $ey = CORE::exp($y); my $ey_1 = 1 / $ey; - return (ref $z)->make(cos($x) * ($ey + $ey_1)/2, sin($x) * ($ey_1 - $ey)/2); + return (ref $z)->make(CORE::cos($x) * ($ey + $ey_1)/2, + CORE::sin($x) * ($ey_1 - $ey)/2); } # @@ -450,9 +748,10 @@ sub cos { sub sin { my ($z) = @_; my ($x, $y) = @{$z->cartesian}; - my $ey = exp($y); + my $ey = CORE::exp($y); my $ey_1 = 1 / $ey; - return (ref $z)->make(sin($x) * ($ey + $ey_1)/2, cos($x) * ($ey - $ey_1)/2); + return (ref $z)->make(CORE::sin($x) * ($ey + $ey_1)/2, + CORE::cos($x) * ($ey - $ey_1)/2); } # @@ -462,29 +761,81 @@ sub sin { # sub tan { my ($z) = @_; - return sin($z) / cos($z); + my $cz = CORE::cos($z); + _divbyzero "tan($z)", "cos($z)" if (CORE::abs($cz) < $eps); + return CORE::sin($z) / $cz; } # -# cotan +# sec +# +# Computes the secant sec(z) = 1 / cos(z). +# +sub sec { + my ($z) = @_; + my $cz = CORE::cos($z); + _divbyzero "sec($z)", "cos($z)" if ($cz == 0); + return 1 / $cz; +} + # -# Computes cotan(z) = 1 / tan(z). +# csc # -sub cotan { +# Computes the cosecant csc(z) = 1 / sin(z). +# +sub csc { my ($z) = @_; - return cos($z) / sin($z); + my $sz = CORE::sin($z); + _divbyzero "csc($z)", "sin($z)" if ($sz == 0); + return 1 / $sz; } # +# cosec +# +# Alias for csc(). +# +sub cosec { Math::Complex::csc(@_) } + +# +# cot +# +# Computes cot(z) = cos(z) / sin(z). +# +sub cot { + my ($z) = @_; + my $sz = CORE::sin($z); + _divbyzero "cot($z)", "sin($z)" if ($sz == 0); + return CORE::cos($z) / $sz; +} + +# +# cotan +# +# Alias for cot(). +# +sub cotan { Math::Complex::cot(@_) } + +# # acos # # Computes the arc cosine acos(z) = -i log(z + sqrt(z*z-1)). # sub acos { - my ($z) = @_; - my $cz = $z*$z - 1; - $cz = cplx($cz, 0) if !ref $cz && $cz < 0; # Force complex if <0 - return ~i * log($z + sqrt $cz); # ~i is -i + my $z = $_[0]; + return CORE::atan2(CORE::sqrt(1-$z*$z), $z) if (! ref $z) && CORE::abs($z) <= 1; + my ($x, $y) = ref $z ? @{$z->cartesian} : ($z, 0); + my $t1 = CORE::sqrt(($x+1)*($x+1) + $y*$y); + my $t2 = CORE::sqrt(($x-1)*($x-1) + $y*$y); + my $alpha = ($t1 + $t2)/2; + my $beta = ($t1 - $t2)/2; + $alpha = 1 if $alpha < 1; + if ($beta > 1) { $beta = 1 } + elsif ($beta < -1) { $beta = -1 } + my $u = CORE::atan2(CORE::sqrt(1-$beta*$beta), $beta); + my $v = CORE::log($alpha + CORE::sqrt($alpha*$alpha-1)); + $v = -$v if $y > 0 || ($y == 0 && $x < -1); + return $package->make($u, $v); } # @@ -493,44 +844,104 @@ sub acos { # Computes the arc sine asin(z) = -i log(iz + sqrt(1-z*z)). # sub asin { - my ($z) = @_; - my $cz = 1 - $z*$z; - $cz = cplx($cz, 0) if !ref $cz && $cz < 0; # Force complex if <0 - return ~i * log(i * $z + sqrt $cz); # ~i is -i + my $z = $_[0]; + return CORE::atan2($z, CORE::sqrt(1-$z*$z)) if (! ref $z) && CORE::abs($z) <= 1; + my ($x, $y) = ref $z ? @{$z->cartesian} : ($z, 0); + my $t1 = CORE::sqrt(($x+1)*($x+1) + $y*$y); + my $t2 = CORE::sqrt(($x-1)*($x-1) + $y*$y); + my $alpha = ($t1 + $t2)/2; + my $beta = ($t1 - $t2)/2; + $alpha = 1 if $alpha < 1; + if ($beta > 1) { $beta = 1 } + elsif ($beta < -1) { $beta = -1 } + my $u = CORE::atan2($beta, CORE::sqrt(1-$beta*$beta)); + my $v = -CORE::log($alpha + CORE::sqrt($alpha*$alpha-1)); + $v = -$v if $y > 0 || ($y == 0 && $x < -1); + return $package->make($u, $v); } # # atan # -# Computes the arc tagent atan(z) = i/2 log((i+z) / (i-z)). +# Computes the arc tangent atan(z) = i/2 log((i+z) / (i-z)). # sub atan { my ($z) = @_; - return i/2 * log((i + $z) / (i - $z)); + return CORE::atan2($z, 1) unless ref $z; + _divbyzero "atan(i)" if ( $z == i); + _divbyzero "atan(-i)" if (-$z == i); + my $log = CORE::log((i + $z) / (i - $z)); + $ip2 = 0.5 * i unless defined $ip2; + return $ip2 * $log; } # -# acotan +# asec # -# Computes the arc cotangent acotan(z) = -i/2 log((i+z) / (z-i)) +# Computes the arc secant asec(z) = acos(1 / z). # -sub acotan { +sub asec { my ($z) = @_; - return i/-2 * log((i + $z) / ($z - i)); + _divbyzero "asec($z)", $z if ($z == 0); + return acos(1 / $z); } # +# acsc +# +# Computes the arc cosecant acsc(z) = asin(1 / z). +# +sub acsc { + my ($z) = @_; + _divbyzero "acsc($z)", $z if ($z == 0); + return asin(1 / $z); +} + +# +# acosec +# +# Alias for acsc(). +# +sub acosec { Math::Complex::acsc(@_) } + +# +# acot +# +# Computes the arc cotangent acot(z) = atan(1 / z) +# +sub acot { + my ($z) = @_; + _divbyzero "acot(0)" if (CORE::abs($z) < $eps); + return ($z >= 0) ? CORE::atan2(1, $z) : CORE::atan2(-1, -$z) unless ref $z; + _divbyzero "acot(i)" if (CORE::abs($z - i) < $eps); + _logofzero "acot(-i)" if (CORE::abs($z + i) < $eps); + return atan(1 / $z); +} + +# +# acotan +# +# Alias for acot(). +# +sub acotan { Math::Complex::acot(@_) } + +# # cosh # # Computes the hyperbolic cosine cosh(z) = (exp(z) + exp(-z))/2. # sub cosh { my ($z) = @_; - my ($x, $y) = ref $z ? @{$z->cartesian} : ($z); - my $ex = exp($x); + my $ex; + unless (ref $z) { + $ex = CORE::exp($z); + return ($ex + 1/$ex)/2; + } + my ($x, $y) = @{$z->cartesian}; + $ex = CORE::exp($x); my $ex_1 = 1 / $ex; - return ($ex + $ex_1)/2 unless ref $z; - return (ref $z)->make(cos($y) * ($ex + $ex_1)/2, sin($y) * ($ex - $ex_1)/2); + return (ref $z)->make(CORE::cos($y) * ($ex + $ex_1)/2, + CORE::sin($y) * ($ex - $ex_1)/2); } # @@ -540,11 +951,16 @@ sub cosh { # sub sinh { my ($z) = @_; - my ($x, $y) = ref $z ? @{$z->cartesian} : ($z); - my $ex = exp($x); + my $ex; + unless (ref $z) { + $ex = CORE::exp($z); + return ($ex - 1/$ex)/2; + } + my ($x, $y) = @{$z->cartesian}; + $ex = CORE::exp($x); my $ex_1 = 1 / $ex; - return ($ex - $ex_1)/2 unless ref $z; - return (ref $z)->make(cos($y) * ($ex - $ex_1)/2, sin($y) * ($ex + $ex_1)/2); + return (ref $z)->make(CORE::cos($y) * ($ex - $ex_1)/2, + CORE::sin($y) * ($ex + $ex_1)/2); } # @@ -554,29 +970,78 @@ sub sinh { # sub tanh { my ($z) = @_; - return sinh($z) / cosh($z); + my $cz = cosh($z); + _divbyzero "tanh($z)", "cosh($z)" if ($cz == 0); + return sinh($z) / $cz; } # -# cotanh +# sech # -# Comptutes the hyperbolic cotangent cotanh(z) = cosh(z) / sinh(z). +# Computes the hyperbolic secant sech(z) = 1 / cosh(z). # -sub cotanh { +sub sech { my ($z) = @_; - return cosh($z) / sinh($z); + my $cz = cosh($z); + _divbyzero "sech($z)", "cosh($z)" if ($cz == 0); + return 1 / $cz; } # +# csch +# +# Computes the hyperbolic cosecant csch(z) = 1 / sinh(z). +# +sub csch { + my ($z) = @_; + my $sz = sinh($z); + _divbyzero "csch($z)", "sinh($z)" if ($sz == 0); + return 1 / $sz; +} + +# +# cosech +# +# Alias for csch(). +# +sub cosech { Math::Complex::csch(@_) } + +# +# coth +# +# Computes the hyperbolic cotangent coth(z) = cosh(z) / sinh(z). +# +sub coth { + my ($z) = @_; + my $sz = sinh($z); + _divbyzero "coth($z)", "sinh($z)" if ($sz == 0); + return cosh($z) / $sz; +} + +# +# cotanh +# +# Alias for coth(). +# +sub cotanh { Math::Complex::coth(@_) } + +# # acosh # # Computes the arc hyperbolic cosine acosh(z) = log(z + sqrt(z*z-1)). # sub acosh { my ($z) = @_; - my $cz = $z*$z - 1; - $cz = cplx($cz, 0) if !ref $cz && $cz < 0; # Force complex if <0 - return log($z + sqrt $cz); + unless (ref $z) { + return CORE::log($z + CORE::sqrt($z*$z-1)) if $z >= 1; + $z = cplx($z, 0); + } + my ($re, $im) = @{$z->cartesian}; + if ($im == 0) { + return cplx(CORE::log($re + CORE::sqrt($re*$re - 1)), 0) if $re >= 1; + return cplx(0, CORE::atan2(CORE::sqrt(1-$re*$re), $re)) if CORE::abs($re) <= 1; + } + return CORE::log($z + CORE::sqrt($z*$z - 1)); } # @@ -586,8 +1051,7 @@ sub acosh { # sub asinh { my ($z) = @_; - my $cz = $z*$z + 1; # Already complex if <0 - return log($z + sqrt $cz); + return CORE::log($z + CORE::sqrt($z*$z + 1)); } # @@ -597,41 +1061,92 @@ sub asinh { # sub atanh { my ($z) = @_; - my $cz = (1 + $z) / (1 - $z); - $cz = cplx($cz, 0) if !ref $cz && $cz < 0; # Force complex if <0 - return log($cz) / 2; + unless (ref $z) { + return CORE::log((1 + $z)/(1 - $z))/2 if CORE::abs($z) < 1; + $z = cplx($z, 0); + } + _divbyzero 'atanh(1)', "1 - $z" if ($z == 1); + _logofzero 'atanh(-1)' if ($z == -1); + return 0.5 * CORE::log((1 + $z) / (1 - $z)); } # -# acotanh +# asech +# +# Computes the hyperbolic arc secant asech(z) = acosh(1 / z). +# +sub asech { + my ($z) = @_; + _divbyzero 'asech(0)', $z if ($z == 0); + return acosh(1 / $z); +} + # -# Computes the arc hyperbolic cotangent acotanh(z) = 1/2 log((1+z) / (z-1)). +# acsch # -sub acotanh { +# Computes the hyperbolic arc cosecant acsch(z) = asinh(1 / z). +# +sub acsch { my ($z) = @_; - my $cz = (1 + $z) / ($z - 1); - $cz = cplx($cz, 0) if !ref $cz && $cz < 0; # Force complex if <0 - return log($cz) / 2; + _divbyzero 'acsch(0)', $z if ($z == 0); + return asinh(1 / $z); } # +# acosech +# +# Alias for acosh(). +# +sub acosech { Math::Complex::acsch(@_) } + +# +# acoth +# +# Computes the arc hyperbolic cotangent acoth(z) = 1/2 log((1+z) / (z-1)). +# +sub acoth { + my ($z) = @_; + _divbyzero 'acoth(0)' if (CORE::abs($z) < $eps); + unless (ref $z) { + return CORE::log(($z + 1)/($z - 1))/2 if CORE::abs($z) > 1; + $z = cplx($z, 0); + } + _divbyzero 'acoth(1)', "$z - 1" if (CORE::abs($z - 1) < $eps); + _logofzero 'acoth(-1)', "1 / $z" if (CORE::abs($z + 1) < $eps); + return CORE::log((1 + $z) / ($z - 1)) / 2; +} + +# +# acotanh +# +# Alias for acot(). +# +sub acotanh { Math::Complex::acoth(@_) } + +# # (atan2) # # Compute atan(z1/z2). # sub atan2 { my ($z1, $z2, $inverted) = @_; - my ($re1, $im1) = @{$z1->cartesian}; - my ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2); - my $tan; - if (defined $inverted && $inverted) { # atan(z2/z1) - return pi * ($re2 > 0 ? 1 : -1) if $re1 == 0 && $im1 == 0; - $tan = $z2 / $z1; + my ($re1, $im1, $re2, $im2); + if ($inverted) { + ($re1, $im1) = ref $z2 ? @{$z2->cartesian} : ($z2, 0); + ($re2, $im2) = @{$z1->cartesian}; } else { - return pi * ($re1 > 0 ? 1 : -1) if $re2 == 0 && $im2 == 0; - $tan = $z1 / $z2; + ($re1, $im1) = @{$z1->cartesian}; + ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0); + } + if ($im2 == 0) { + return cplx(CORE::atan2($re1, $re2), 0) if $im1 == 0; + return cplx(($im1<=>0) * pip2, 0) if $re2 == 0; } - return atan($tan); + my $w = atan($z1/$z2); + my ($u, $v) = ref $w ? @{$w->cartesian} : ($w, 0); + $u += pi if $re2 < 0; + $u -= pit2 if $u > pi; + return cplx($u, $v); } # @@ -639,7 +1154,7 @@ sub atan2 { # ->display_format # # Set (fetch if no argument) display format for all complex numbers that -# don't happen to have overrriden it via ->display_format +# don't happen to have overridden it via ->display_format # # When called as a method, this actually sets the display format for # the current object. @@ -653,7 +1168,7 @@ sub display_format { if (ref $self) { # Called as a method $format = shift; - } else { # Regular procedure call + } else { # Regular procedure call $format = $self; undef $self; } @@ -699,21 +1214,57 @@ sub stringify_cartesian { my ($x, $y) = @{$z->cartesian}; my ($re, $im); - $re = "$x" if abs($x) >= 1e-14; - if ($y == 1) { $im = 'i' } - elsif ($y == -1) { $im = '-i' } - elsif (abs($y) >= 1e-14) { $im = "${y}i" } + $x = int($x + ($x < 0 ? -1 : 1) * $eps) + if int(CORE::abs($x)) != int(CORE::abs($x) + $eps); + $y = int($y + ($y < 0 ? -1 : 1) * $eps) + if int(CORE::abs($y)) != int(CORE::abs($y) + $eps); - my $str; + $re = "$x" if CORE::abs($x) >= $eps; + if ($y == 1) { $im = 'i' } + elsif ($y == -1) { $im = '-i' } + elsif (CORE::abs($y) >= $eps) { $im = $y . "i" } + + my $str = ''; $str = $re if defined $re; $str .= "+$im" if defined $im; $str =~ s/\+-/-/; $str =~ s/^\+//; + $str =~ s/([-+])1i/$1i/; # Not redundant with the above 1/-1 tests. $str = '0' unless $str; return $str; } + +# Helper for stringify_polar, a Greatest Common Divisor with a memory. + +sub _gcd { + my ($a, $b) = @_; + + use integer; + + # Loops forever if given negative inputs. + + if ($b and $a > $b) { return gcd($a % $b, $b) } + elsif ($a and $b > $a) { return gcd($b % $a, $a) } + else { return $a ? $a : $b } +} + +my %gcd; + +sub gcd { + my ($a, $b) = @_; + + my $id = "$a $b"; + + unless (exists $gcd{$id}) { + $gcd{$id} = _gcd($a, $b); + $gcd{"$b $a"} = $gcd{$id}; + } + + return $gcd{$id}; +} + # # ->stringify_polar # @@ -724,35 +1275,59 @@ sub stringify_polar { my ($r, $t) = @{$z->polar}; my $theta; - return '[0,0]' if $r <= 1e-14; + return '[0,0]' if $r <= $eps; - my $tpi = 2 * pi; - my $nt = $t / $tpi; - $nt = ($nt - int($nt)) * $tpi; - $nt += $tpi if $nt < 0; # Range [0, 2pi] + my $nt = $t / pit2; + $nt = ($nt - int($nt)) * pit2; + $nt += pit2 if $nt < 0; # Range [0, 2pi] - if (abs($nt) <= 1e-14) { $theta = 0 } - elsif (abs(pi-$nt) <= 1e-14) { $theta = 'pi' } + if (CORE::abs($nt) <= $eps) { $theta = 0 } + elsif (CORE::abs(pi-$nt) <= $eps) { $theta = 'pi' } - return "\[$r,$theta\]" if defined $theta; + if (defined $theta) { + $r = int($r + ($r < 0 ? -1 : 1) * $eps) + if int(CORE::abs($r)) != int(CORE::abs($r) + $eps); + $theta = int($theta + ($theta < 0 ? -1 : 1) * $eps) + if ($theta ne 'pi' and + int(CORE::abs($theta)) != int(CORE::abs($theta) + $eps)); + return "\[$r,$theta\]"; + } # # Okay, number is not a real. Try to identify pi/n and friends... # - $nt -= $tpi if $nt > pi; - my ($n, $k, $kpi); - - for ($k = 1, $kpi = pi; $k < 10; $k++, $kpi += pi) { + $nt -= pit2 if $nt > pi; + + if (CORE::abs($nt) >= deg1) { + my ($n, $k, $kpi); + + for ($k = 1, $kpi = pi; $k < 10; $k++, $kpi += pi) { $n = int($kpi / $nt + ($nt > 0 ? 1 : -1) * 0.5); - if (abs($kpi/$n - $nt) <= 1e-14) { - $theta = ($nt < 0 ? '-':'').($k == 1 ? 'pi':"${k}pi").'/'.abs($n); - last; + if (CORE::abs($kpi/$n - $nt) <= $eps) { + $n = CORE::abs($n); + my $gcd = gcd($k, $n); + if ($gcd > 1) { + $k /= $gcd; + $n /= $gcd; + } + next if $n > 360; + $theta = ($nt < 0 ? '-':''). + ($k == 1 ? 'pi':"${k}pi"); + $theta .= '/'.$n if $n > 1; + last; } + } } $theta = $nt unless defined $theta; + $r = int($r + ($r < 0 ? -1 : 1) * $eps) + if int(CORE::abs($r)) != int(CORE::abs($r) + $eps); + $theta = int($theta + ($theta < 0 ? -1 : 1) * $eps) + if ($theta !~ m(^-?\d*pi/\d+$) and + int(CORE::abs($theta)) != int(CORE::abs($theta) + $eps)); + return "\[$r,$theta\]"; } @@ -766,6 +1341,7 @@ Math::Complex - complex numbers and associated mathematical functions =head1 SYNOPSIS use Math::Complex; + $z = Math::Complex->make(5, 6); $t = 4 - 3*i + $z; $j = cplxe(1, 2*pi/3); @@ -842,7 +1418,7 @@ between this form and the cartesian form C is immediate: which is also expressed by this formula: - z = rho * exp(i * theta) = rho * (cos theta + i * sin theta) + z = rho * exp(i * theta) = rho * (cos theta + i * sin theta) In other words, it's the projection of the vector onto the I and I axes. Mathematicians call I the I or I and I @@ -852,8 +1428,8 @@ noted C. The polar notation (also known as the trigonometric representation) is much more handy for performing multiplications and divisions of complex numbers, whilst the cartesian notation is better -suited for additions and substractions. Real numbers are on the I -axis, and therefore I is zero. +suited for additions and subtractions. Real numbers are on the I +axis, and therefore I is zero or I. All the common operations that can be performed on a real number have been defined to work on complex numbers as well, and are merely @@ -862,8 +1438,8 @@ they keep their natural meaning when there is no imaginary part, provided the number is within their definition set. For instance, the C routine which computes the square root of -its argument is only defined for positive real numbers and yields a -positive real number (it is an application from B to B). +its argument is only defined for non-negative real numbers and yields a +non-negative real number (it is an application from B to B). If we allow it to return a complex number, then it can be extended to negative real numbers to become an application from B to B (the set of complex numbers): @@ -876,14 +1452,15 @@ the following definition: sqrt(z = [r,t]) = sqrt(r) * exp(i * t/2) -Indeed, a negative real number can be noted C<[x,pi]> -(the modulus I is always positive, so C<[x,pi]> is really C<-x>, a -negative number) -and the above definition states that +Indeed, a negative real number can be noted C<[x,pi]> (the modulus +I is always non-negative, so C<[x,pi]> is really C<-x>, a negative +number) and the above definition states that sqrt([x,pi]) = sqrt(x) * exp(i*pi/2) = [sqrt(x),pi/2] = sqrt(x)*i which is exactly what we had defined for negative real numbers above. +The C returns only one of the solutions: if you want the both, +use the C function. All the common mathematical functions defined on real numbers that are extended to complex numbers share that same property of working @@ -936,14 +1513,13 @@ the following (overloaded) operations are supported on complex numbers: z1 * z2 = (r1 * r2) * exp(i * (t1 + t2)) z1 / z2 = (r1 / r2) * exp(i * (t1 - t2)) z1 ** z2 = exp(z2 * log z1) - ~z1 = a - bi - abs(z1) = r1 = sqrt(a*a + b*b) - sqrt(z1) = sqrt(r1) * exp(i * t1/2) - exp(z1) = exp(a) * exp(i * b) - log(z1) = log(r1) + i*t1 - sin(z1) = 1/2i (exp(i * z1) - exp(-i * z1)) - cos(z1) = 1/2 (exp(i * z1) + exp(-i * z1)) - abs(z1) = r1 + ~z = a - bi + abs(z) = r1 = sqrt(a*a + b*b) + sqrt(z) = sqrt(r1) * exp(i * t/2) + exp(z) = exp(a) * exp(i * b) + log(z) = log(r1) + i*t + sin(z) = 1/2i (exp(i * z1) - exp(-i * z)) + cos(z) = 1/2 (exp(i * z1) + exp(-i * z)) atan2(z1, z2) = atan(z1/z2) The following extra operations are supported on both real and complex @@ -952,30 +1528,51 @@ numbers: Re(z) = a Im(z) = b arg(z) = t + abs(z) = r cbrt(z) = z ** (1/3) log10(z) = log(z) / log(10) logn(z, n) = log(z) / log(n) tan(z) = sin(z) / cos(z) - cotan(z) = 1 / tan(z) + + csc(z) = 1 / sin(z) + sec(z) = 1 / cos(z) + cot(z) = 1 / tan(z) asin(z) = -i * log(i*z + sqrt(1-z*z)) - acos(z) = -i * log(z + sqrt(z*z-1)) + acos(z) = -i * log(z + i*sqrt(1-z*z)) atan(z) = i/2 * log((i+z) / (i-z)) - acotan(z) = -i/2 * log((i+z) / (z-i)) + + acsc(z) = asin(1 / z) + asec(z) = acos(1 / z) + acot(z) = atan(1 / z) = -i/2 * log((i+z) / (z-i)) sinh(z) = 1/2 (exp(z) - exp(-z)) cosh(z) = 1/2 (exp(z) + exp(-z)) - tanh(z) = sinh(z) / cosh(z) - cotanh(z) = 1 / tanh(z) - + tanh(z) = sinh(z) / cosh(z) = (exp(z) - exp(-z)) / (exp(z) + exp(-z)) + + csch(z) = 1 / sinh(z) + sech(z) = 1 / cosh(z) + coth(z) = 1 / tanh(z) + asinh(z) = log(z + sqrt(z*z+1)) acosh(z) = log(z + sqrt(z*z-1)) atanh(z) = 1/2 * log((1+z) / (1-z)) - acotanh(z) = 1/2 * log((1+z) / (z-1)) -The I function is available to compute all the Ith + acsch(z) = asinh(1 / z) + asech(z) = acosh(1 / z) + acoth(z) = atanh(1 / z) = 1/2 * log((1+z) / (z-1)) + +I, I, I, I, I, I, I, I, +I, I, I, have aliases I, I, I, +I, I, I, I, I, I, +I, I, respectively. C, C, C, C, +C, and C can be used also also mutators. The C +returns only one of the solutions: if you want all three, use the +C function. + +The I function is available to compute all the I roots of some complex, where I is a strictly positive integer. There are exactly I such roots, returned as a list. Getting the number mathematicians call C such that: @@ -990,10 +1587,11 @@ The Ith root for C is given by: (root(z, n))[k] = r**(1/n) * exp(i * (t + 2*k*pi)/n) -The I operation is also defined. In order to ensure its -restriction to real numbers is conform to what you would expect, the -comparison is run on the real part of the complex number first, -and imaginary parts are compared only when the real parts match. +The I comparison operator, E=E, is also defined. In +order to ensure its restriction to real numbers is conform to what you +would expect, the comparison is run on the real part of the complex +number first, and imaginary parts are compared only when the real +parts match. =head1 CREATION @@ -1006,23 +1604,30 @@ if you know the cartesian form of the number, or $z = 3 + 4*i; -if you like. To create a number using the trigonometric form, use either: +if you like. To create a number using the polar form, use either: $z = Math::Complex->emake(5, pi/3); $x = cplxe(5, pi/3); -instead. The first argument is the modulus, the second is the angle (in radians). -(Mnmemonic: C is used as a notation for complex numbers in the trigonometric -form). +instead. The first argument is the modulus, the second is the angle +(in radians, the full circle is 2*pi). (Mnemonic: C is used as a +notation for complex numbers in the polar form). It is possible to write: $x = cplxe(-3, pi/4); but that will be silently converted into C<[3,-3pi/4]>, since the modulus -must be positive (it represents the distance to the origin in the complex +must be non-negative (it represents the distance to the origin in the complex plane). +It is also possible to have a complex number as either argument of +either the C or C: the appropriate component of +the argument will be used. + + $z1 = cplx(-2, 1); + $z2 = cplx($z1, 4); + =head1 STRINGIFICATION When printed, a complex number is usually shown under its cartesian @@ -1071,23 +1676,83 @@ Here are some examples: $k = exp(i * 2*pi/3); print "$j - $k = ", $j - $k, "\n"; -=head1 BUGS + $z->Re(3); # Re, Im, arg, abs, + $j->arg(2); # (the last two aka rho, theta) + # can be used also as mutators. + +=head1 ERRORS DUE TO DIVISION BY ZERO OR LOGARITHM OF ZERO + +The division (/) and the following functions -Saying C exports many mathematical routines in the caller -environment. This is construed as a feature by the Author, actually... ;-) + log ln log10 logn + tan sec csc cot + atan asec acsc acot + tanh sech csch coth + atanh asech acsch acoth -The code is not optimized for speed, although we try to use the cartesian -form for addition-like operators and the trigonometric form for all -multiplication-like operators. +cannot be computed for all arguments because that would mean dividing +by zero or taking logarithm of zero. These situations cause fatal +runtime errors looking like this -The arg() routine does not ensure the angle is within the range [-pi,+pi] -(a side effect caused by multiplication and division using the trigonometric -representation). + cot(0): Division by zero. + (Because in the definition of cot(0), the divisor sin(0) is 0) + Died at ... + +or + + atanh(-1): Logarithm of zero. + Died at... + +For the C, C, C, C, C, C, C, +C, C, the argument cannot be C<0> (zero). For the the +logarithmic functions and the C, C, the argument cannot +be C<1> (one). For the C, C, the argument cannot be +C<-1> (minus one). For the C, C, the argument cannot be +C (the imaginary unit). For the C, C, the argument +cannot be C<-i> (the negative imaginary unit). For the C, +C, C, the argument cannot be I, where I +is any integer. + +Note that because we are operating on approximations of real numbers, +these errors can happen when merely `too close' to the singularities +listed above. For example C will die of +division by zero. + +=head1 ERRORS DUE TO INDIGESTIBLE ARGUMENTS + +The C and C accept both real and complex arguments. +When they cannot recognize the arguments they will die with error +messages like the following + + Math::Complex::make: Cannot take real part of ... + Math::Complex::make: Cannot take real part of ... + Math::Complex::emake: Cannot take rho of ... + Math::Complex::emake: Cannot take theta of ... + +=head1 BUGS + +Saying C exports many mathematical routines in the +caller environment and even overrides some (C, C). +This is construed as a feature by the Authors, actually... ;-) All routines expect to be given real or complex numbers. Don't attempt to use BigFloat, since Perl has currently no rule to disambiguate a '+' operation (for instance) between two overloaded entities. -=head1 AUTHOR +In Cray UNICOS there is some strange numerical instability that results +in root(), cos(), sin(), cosh(), sinh(), losing accuracy fast. Beware. +The bug may be in UNICOS math libs, in UNICOS C compiler, in Math::Complex. +Whatever it is, it does not manifest itself anywhere else where Perl runs. + +=head1 AUTHORS + +Raphael Manfredi > and +Jarkko Hietaniemi >. + +Extensive patches by Daniel S. Lewart >. + +=cut + +1; -Raphael Manfredi > +# eof