X-Git-Url: http://git.shadowcat.co.uk/gitweb/gitweb.cgi?a=blobdiff_plain;f=lib%2FMath%2FComplex.pm;h=5b69039afc07ccb1bdcd2e0c6081fb43444944c5;hb=f4c556ac9d141bf86702c68d95acad2db5ec6874;hp=a501b99ddd29a0627115a1c2debb8bec946a38b0;hpb=0e505df1d74bd65a42a69dbec1b44eebf2bd0b93;p=p5sagit%2Fp5-mst-13.2.git diff --git a/lib/Math/Complex.pm b/lib/Math/Complex.pm index a501b99..5b69039 100644 --- a/lib/Math/Complex.pm +++ b/lib/Math/Complex.pm @@ -1,26 +1,26 @@ -# $RCSFile$ # # Complex numbers and associated mathematical functions -# -- Raphael Manfredi, September 1996 -# -- Jarkko Hietaniemi, March-April 1997 +# -- Raphael Manfredi Since Sep 1996 +# -- Jarkko Hietaniemi Since Mar 1997 +# -- Daniel S. Lewart Since Sep 1997 +# require Exporter; package Math::Complex; use strict; -use vars qw($VERSION @ISA - @EXPORT %EXPORT_TAGS - $package $display - $i $logn %logn); +use vars qw($VERSION @ISA @EXPORT %EXPORT_TAGS); -@ISA = qw(Exporter); +my ( $i, $ip2, %logn ); + +$VERSION = sprintf("%s", q$Id: Complex.pm,v 1.26 1998/11/01 00:00:00 dsl Exp $ =~ /(\d+\.\d+)/); -$VERSION = 1.01; +@ISA = qw(Exporter); my @trig = qw( pi - sin cos tan + tan csc cosec sec cot cotan asin acos atan acsc acosec asec acot acotan @@ -31,8 +31,8 @@ my @trig = qw( ); @EXPORT = (qw( - i Re Im arg - sqrt exp log ln + i Re Im rho theta arg + sqrt log ln log10 logn cbrt root cplx cplxe ), @@ -62,11 +62,12 @@ use overload 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): @@ -77,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 # @@ -85,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; } @@ -99,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; } @@ -133,18 +177,38 @@ sub cplxe { # # pi # -# The number defined as 2 * pi = 360 degrees +# The number defined as pi = 180 degrees # +use constant pi => 4 * CORE::atan2(1, 1); -use constant pi => 4 * atan2(1, 1); +# +# pit2 +# +# The full circle +# +use constant pit2 => 2 * pi; # -# log2inv +# pip2 # -# Used in log10(). +# The quarter circle # +use constant pip2 => pi / 2; -use constant log10inv => 1 / log(10); +# +# 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 @@ -155,7 +219,7 @@ sub i () { return $i if ($i); $i = bless {}; $i->{'cartesian'} = [0, 1]; - $i->{'polar'} = [1, pi/2]; + $i->{'polar'} = [1, pip2]; $i->{c_dirty} = 0; $i->{p_dirty} = 0; return $i; @@ -182,7 +246,7 @@ sub update_cartesian { my $self = shift; 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)]; } # @@ -196,7 +260,7 @@ sub update_polar { 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'} = [CORE::sqrt($x*$x + $y*$y), CORE::atan2($y, $x)]; } # @@ -242,15 +306,28 @@ sub minus { # Computes z1*z2. # sub multiply { - my ($z1, $z2, $regular) = @_; - my ($r1, $t1) = @{$z1->polar}; - $z2 = cplxe(abs($z2), $z2 >= 0 ? 0 : pi) unless ref $z2; - my ($r2, $t2) = @{$z2->polar}; - 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); } # @@ -268,7 +345,7 @@ sub _divbyzero { } my @up = caller(1); - + $mess .= "Died at $up[1] line $up[2].\n"; die $mess; @@ -281,63 +358,75 @@ sub _divbyzero { # sub divide { my ($z1, $z2, $inverted) = @_; - my ($r1, $t1) = @{$z1->polar}; - $z2 = cplxe(abs($z2), $z2 >= 0 ? 0 : pi) unless ref $z2; - my ($r2, $t2) = @{$z2->polar}; - unless (defined $inverted) { - _divbyzero "$z1/0" if ($r2 == 0); - $z1->set_polar([$r1 / $r2, $t1 - $t2]); - return $z1; - } - if ($inverted) { + 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); - return (ref $z1)->emake($r2 / $r1, $t2 - $t1); - } else { + $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); - return (ref $z1)->emake($r1 / $r2, $t1 - $t2); + $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); + } + } } } # -# _zerotozero -# -# Die on zero raised to the zeroth. -# -sub _zerotozero { - my $mess = "The zero raised to the zeroth power is not defined.\n"; - - my @up = caller(1); - - $mess .= "Died at $up[1] line $up[2].\n"; - - die $mess; -} - -# # (power) # # Computes z1**z2 = exp(z2 * log z1)). # sub power { my ($z1, $z2, $inverted) = @_; - _zerotozero if ($z1 == 0 and $z2 == 0); - return 1 if ($z1 == 1); - return 0 if ($z2 == 0); - $z2 = cplx($z2) unless ref $z2; - unless (defined $inverted) { - my $z3 = exp($z2 * log $z1); - $z1->set_cartesian([@{$z3->cartesian}]); - return $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; } - return exp($z2 * log $z1) unless $inverted; - return exp($z1 * log $z2); + 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) = @_; @@ -357,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); @@ -381,25 +471,46 @@ sub conjugate { # # (abs) # -# Compute complex's norm (rho). +# Compute or set complex's norm (rho). # sub abs { - my ($z) = @_; - return abs($z) unless ref $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 ($z < 0 ? pi : 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; } # @@ -407,11 +518,22 @@ sub arg { # # 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) = @_; - $z = cplx($z, 0) unless ref $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); } # @@ -419,11 +541,14 @@ sub sqrt { # # Compute cbrt(z) (cubic root). # +# Why are we not returning three values? The same answer as for sqrt(). +# sub cbrt { my ($z) = @_; - return cplx($z, 0) ** (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); } # @@ -435,7 +560,7 @@ 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; @@ -454,15 +579,17 @@ sub _rootbad { sub root { my ($z, $n) = @_; _rootbad($n) if ($n < 1 or int($n) != $n); - my ($r, $t) = ref $z ? @{$z->polar} : (abs($z), $z >= 0 ? 0 : pi); + 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; } @@ -470,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(@_); } # @@ -498,9 +653,29 @@ sub Im { # sub exp { my ($z) = @_; - $z = cplx($z, 0) unless ref $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; } # @@ -510,12 +685,15 @@ sub exp { # sub log { my ($z) = @_; - $z = cplx($z, 0) unless ref $z; - my ($x, $y) = @{$z->cartesian}; + unless (ref $z) { + _logofzero("log") if $z == 0; + return $z > 0 ? CORE::log($z) : cplx(CORE::log(-$z), pi); + } my ($r, $t) = @{$z->polar}; - $t -= 2 * pi if ($t > pi() and $x < 0); - $t += 2 * pi if ($t < -pi() and $x < 0); - 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); } # @@ -532,11 +710,7 @@ sub ln { Math::Complex::log(@_) } # sub log10 { - my ($z) = @_; - - return log(cplx($z, 0)) * log10inv unless ref $z; - my ($r, $t) = @{$z->polar}; - return (ref $z)->make(log($r) * log10inv, $t * log10inv); + return Math::Complex::log($_[0]) * uplog10; } # @@ -548,8 +722,8 @@ 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) / $logn; + $logn = $logn{$n} = CORE::log($n) unless defined $logn; # Cache log(n) + return CORE::log($z) / $logn; } # @@ -559,12 +733,11 @@ sub logn { # sub cos { my ($z) = @_; - $z = cplx($z, 0) unless ref $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); } # @@ -574,12 +747,11 @@ sub cos { # sub sin { my ($z) = @_; - $z = cplx($z, 0) unless ref $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); } # @@ -589,9 +761,9 @@ sub sin { # sub tan { my ($z) = @_; - my $cz = cos($z); - _divbyzero "tan($z)", "cos($z)" if ($cz == 0); - return sin($z) / $cz; + my $cz = CORE::cos($z); + _divbyzero "tan($z)", "cos($z)" if (CORE::abs($cz) < $eps); + return CORE::sin($z) / $cz; } # @@ -601,7 +773,7 @@ sub tan { # sub sec { my ($z) = @_; - my $cz = cos($z); + my $cz = CORE::cos($z); _divbyzero "sec($z)", "cos($z)" if ($cz == 0); return 1 / $cz; } @@ -613,7 +785,7 @@ sub sec { # sub csc { my ($z) = @_; - my $sz = sin($z); + my $sz = CORE::sin($z); _divbyzero "csc($z)", "sin($z)" if ($sz == 0); return 1 / $sz; } @@ -628,13 +800,13 @@ sub cosec { Math::Complex::csc(@_) } # # cot # -# Computes cot(z) = 1 / tan(z). +# Computes cot(z) = cos(z) / sin(z). # sub cot { my ($z) = @_; - my $sz = sin($z); + my $sz = CORE::sin($z); _divbyzero "cot($z)", "sin($z)" if ($sz == 0); - return cos($z) / $sz; + return CORE::cos($z) / $sz; } # @@ -650,9 +822,20 @@ sub cotan { Math::Complex::cot(@_) } # Computes the arc cosine acos(z) = -i log(z + sqrt(z*z-1)). # sub acos { - my ($z) = @_; - $z = cplx($z, 0) unless ref $z; - return ~i * log($z + (Re($z) * Im($z) > 0 ? 1 : -1) * sqrt($z*$z - 1)); + 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); } # @@ -661,9 +844,20 @@ sub acos { # Computes the arc sine asin(z) = -i log(iz + sqrt(1-z*z)). # sub asin { - my ($z) = @_; - $z = cplx($z, 0) unless ref $z; - return ~i * log(i * $z + sqrt(1 - $z*$z)); + 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); } # @@ -673,9 +867,12 @@ sub asin { # sub atan { my ($z) = @_; - $z = cplx($z, 0) unless ref $z; - _divbyzero "atan($z)", "i - $z" if ($z == i); - 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; } # @@ -692,7 +889,7 @@ sub asec { # # acsc # -# Computes the arc cosecant sec(z) = asin(1 / z). +# Computes the arc cosecant acsc(z) = asin(1 / z). # sub acsc { my ($z) = @_; @@ -710,13 +907,15 @@ sub acosec { Math::Complex::acsc(@_) } # # acot # -# Computes the arc cotangent acot(z) = -i/2 log((i+z) / (z-i)) +# Computes the arc cotangent acot(z) = atan(1 / z) # sub acot { my ($z) = @_; - $z = cplx($z, 0) unless ref $z; - _divbyzero "acot($z)", "$z - i" if ($z == i); - return i/-2 * log((i + $z) / ($z - i)); + _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); } # @@ -733,17 +932,16 @@ sub acotan { Math::Complex::acot(@_) } # sub cosh { my ($z) = @_; - my $real; + my $ex; unless (ref $z) { - $z = cplx($z, 0); - $real = 1; + $ex = CORE::exp($z); + return ($ex + 1/$ex)/2; } my ($x, $y) = @{$z->cartesian}; - my $ex = exp($x); + $ex = CORE::exp($x); my $ex_1 = 1 / $ex; - return ($ex + $ex_1)/2 if $real; - 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); } # @@ -753,17 +951,16 @@ sub cosh { # sub sinh { my ($z) = @_; - my $real; + my $ex; unless (ref $z) { - $z = cplx($z, 0); - $real = 1; + $ex = CORE::exp($z); + return ($ex - 1/$ex)/2; } my ($x, $y) = @{$z->cartesian}; - my $ex = exp($x); + $ex = CORE::exp($x); my $ex_1 = 1 / $ex; - return ($ex - $ex_1)/2 if $real; - 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); } # @@ -835,8 +1032,16 @@ sub cotanh { Math::Complex::coth(@_) } # sub acosh { my ($z) = @_; - $z = cplx($z, 0) unless ref $z; - return log($z + sqrt($z*$z - 1)); + 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)); } # @@ -846,8 +1051,7 @@ sub acosh { # sub asinh { my ($z) = @_; - $z = cplx($z, 0) unless ref $z; - return log($z + sqrt($z*$z + 1)); + return CORE::log($z + CORE::sqrt($z*$z + 1)); } # @@ -857,10 +1061,13 @@ sub asinh { # sub atanh { my ($z) = @_; - _divbyzero 'atanh(1)', "1 - $z" if ($z == 1); - $z = cplx($z, 0) unless ref $z; - my $cz = (1 + $z) / (1 - $z); - 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)); } # @@ -899,10 +1106,14 @@ sub acosech { Math::Complex::acsch(@_) } # sub acoth { my ($z) = @_; - _divbyzero 'acoth(1)', "$z - 1" if ($z == 1); - $z = cplx($z, 0) unless ref $z; - my $cz = (1 + $z) / ($z - 1); - return log($cz) / 2; + _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; } # @@ -919,17 +1130,23 @@ sub acotanh { Math::Complex::acoth(@_) } # sub atan2 { my ($z1, $z2, $inverted) = @_; - my ($re1, $im1) = ref $z1 ? @{$z1->cartesian} : ($z1, 0); - my ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0); - 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); } - return atan($tan); + if ($im2 == 0) { + return cplx(CORE::atan2($re1, $re2), 0) if $im1 == 0; + return cplx(($im1<=>0) * pip2, 0) if $re2 == 0; + } + 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); } # @@ -937,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. @@ -997,26 +1214,57 @@ sub stringify_cartesian { my ($x, $y) = @{$z->cartesian}; my ($re, $im); - $x = int($x + ($x < 0 ? -1 : 1) * 1e-14) - if int(abs($x)) != int(abs($x) + 1e-14); - $y = int($y + ($y < 0 ? -1 : 1) * 1e-14) - if int(abs($y)) != int(abs($y) + 1e-14); + $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); - $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" } + $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 # @@ -1026,24 +1274,22 @@ sub stringify_polar { my $z = shift; my ($r, $t) = @{$z->polar}; my $theta; - my $eps = 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) <= $eps) { $theta = 0 } - elsif (abs(pi-$nt) <= $eps) { $theta = 'pi' } + if (CORE::abs($nt) <= $eps) { $theta = 0 } + elsif (CORE::abs(pi-$nt) <= $eps) { $theta = 'pi' } if (defined $theta) { $r = int($r + ($r < 0 ? -1 : 1) * $eps) - if int(abs($r)) != int(abs($r) + $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(abs($theta)) != int(abs($theta) + $eps)); + int(CORE::abs($theta)) != int(CORE::abs($theta) + $eps)); return "\[$r,$theta\]"; } @@ -1051,25 +1297,36 @@ sub stringify_polar { # 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) <= $eps) { - $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(abs($r)) != int(abs($r) + $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(abs($theta)) != int(abs($theta) + $eps)); + int(CORE::abs($theta)) != int(CORE::abs($theta) + $eps)); return "\[$r,$theta\]"; } @@ -1084,7 +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); @@ -1161,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 @@ -1171,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 @@ -1181,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): @@ -1195,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 @@ -1255,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 @@ -1271,6 +1528,7 @@ numbers: Re(z) = a Im(z) = b arg(z) = t + abs(z) = r cbrt(z) = z ** (1/3) log10(z) = log(z) / log(10) @@ -1283,12 +1541,12 @@ numbers: 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)) acsc(z) = asin(1 / z) asec(z) = acos(1 / z) - acot(z) = -i/2 * log((i+z) / (z-i)) + 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)) @@ -1297,7 +1555,7 @@ numbers: 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)) @@ -1306,10 +1564,13 @@ numbers: 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, have aliases I, I, I, -I, I, I, I, I, I, -respectively. +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. @@ -1326,11 +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 comparison operator 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 @@ -1343,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, the full circle is 2*pi). (Mnmemonic: C is used as a -notation for complex numbers in the trigonometric form). +(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 @@ -1408,65 +1676,83 @@ Here are some examples: $k = exp(i * 2*pi/3); print "$j - $k = ", $j - $k, "\n"; -=head1 ERRORS DUE TO DIVISION BY ZERO + $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 - tan - sec - csc - cot - asec - acsc - atan - acot - tanh - sech - csch - coth - atanh - asech - acsch - acoth + log ln log10 logn + tan sec csc cot + atan asec acsc acot + tanh sech csch coth + atanh asech acsch acoth cannot be computed for all arguments because that would mean dividing -by zero. These situations cause fatal runtime errors looking like this +by zero or taking logarithm of zero. These situations cause fatal +runtime errors looking like this cot(0): Division by zero. (Because in the definition of cot(0), the divisor sin(0) is 0) Died at ... -For the C, C, C, C, C, C, C, -C, the argument cannot be C<0> (zero). For the C, -C, the argument cannot be C<1> (one). For the C, C, -the argument cannot be C (the imaginary unit). For the C, -C, C, C, the argument cannot be I, where -I is any integer. +or -=head1 BUGS + atanh(-1): Logarithm of zero. + Died at... -Saying C exports many mathematical routines in the -caller environment and even overrides some (C, C, C, -C, C). This is construed as a feature by the Authors, -actually... ;-) +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 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. +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 -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). + 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. +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 > - Jarkko Hietaniemi > +Raphael Manfredi > and +Jarkko Hietaniemi >. + +Extensive patches by Daniel S. Lewart >. =cut +1; + # eof