From: Steve Peters Date: Thu, 17 Jan 2008 12:24:21 +0000 (+0000) Subject: Upgrade to Math-Complex-1.47 X-Git-Url: http://git.shadowcat.co.uk/gitweb/gitweb.cgi?a=commitdiff_plain;h=f3828575f121139b95363dc2ba6aab15ad201a85;p=p5sagit%2Fp5-mst-13.2.git Upgrade to Math-Complex-1.47 p4raw-id: //depot/perl@32989 --- diff --git a/lib/Math/Complex.pm b/lib/Math/Complex.pm index c9af42a..36dc1ea 100644 --- a/lib/Math/Complex.pm +++ b/lib/Math/Complex.pm @@ -9,7 +9,7 @@ package Math::Complex; use vars qw($VERSION @ISA @EXPORT @EXPORT_OK %EXPORT_TAGS $Inf); -$VERSION = 1.44; +$VERSION = 1.47; BEGIN { # For 64-bit doubles, anyway. diff --git a/lib/Math/Complex.t b/lib/Math/Complex.t index 1abb5b5..cbd5ed3 100755 --- a/lib/Math/Complex.t +++ b/lib/Math/Complex.t @@ -13,7 +13,7 @@ BEGIN { } } -use Math::Complex 1.44; +use Math::Complex 1.47; use vars qw($VERSION); diff --git a/lib/Math/Trig.pm b/lib/Math/Trig.pm index 3c864b1..d8ecf69 100644 --- a/lib/Math/Trig.pm +++ b/lib/Math/Trig.pm @@ -10,21 +10,23 @@ package Math::Trig; use 5.005; use strict; -use Math::Complex 1.44; +use Math::Complex 1.47; use Math::Complex qw(:trig :pi); use vars qw($VERSION $PACKAGE @ISA @EXPORT @EXPORT_OK %EXPORT_TAGS); @ISA = qw(Exporter); -$VERSION = 1.09; +$VERSION = 1.12; my @angcnv = qw(rad2deg rad2grad deg2rad deg2grad grad2rad grad2deg); +my @areal = qw(asin_real acos_real); + @EXPORT = (@{$Math::Complex::EXPORT_TAGS{'trig'}}, - @angcnv); + @angcnv, @areal); my @rdlcnv = qw(cartesian_to_cylindrical cartesian_to_spherical @@ -92,6 +94,22 @@ sub rad2grad ($;$) { my $d = _RG * $_[0]; $_[1] ? $d : grad2grad($d) } sub grad2rad ($;$) { my $d = _GR * $_[0]; $_[1] ? $d : rad2rad($d) } +# +# acos and asin functions which always return a real number +# + +sub acos_real { + return 0 if $_[0] >= 1; + return pi if $_[0] <= -1; + return acos($_[0]); +} + +sub asin_real { + return &pip2 if $_[0] >= 1; + return -&pip2 if $_[0] <= -1; + return asin($_[0]); +} + sub cartesian_to_spherical { my ( $x, $y, $z ) = @_; @@ -99,7 +117,7 @@ sub cartesian_to_spherical { return ( $rho, atan2( $y, $x ), - $rho ? acos( $z / $rho ) : 0 ); + $rho ? acos_real( $z / $rho ) : 0 ); } sub spherical_to_cartesian { @@ -141,8 +159,8 @@ sub great_circle_distance { my $lat1 = pip2 - $phi1; return $rho * - acos(cos( $lat0 ) * cos( $lat1 ) * cos( $theta0 - $theta1 ) + - sin( $lat0 ) * sin( $lat1 ) ); + acos_real( cos( $lat0 ) * cos( $lat1 ) * cos( $theta0 - $theta1 ) + + sin( $lat0 ) * sin( $lat1 ) ); } sub great_circle_direction { @@ -154,9 +172,9 @@ sub great_circle_direction { my $lat1 = pip2 - $phi1; my $direction = - acos((sin($lat1) - sin($lat0) * cos($distance)) / - (cos($lat0) * sin($distance))); - + acos_real((sin($lat1) - sin($lat0) * cos($distance)) / + (cos($lat0) * sin($distance))); + $direction = pi2 - $direction if sin($theta1 - $theta0) < 0; @@ -189,7 +207,7 @@ sub great_circle_waypoint { my $z = $A * sin($lat0) + $B * sin($lat1); my $theta = atan2($y, $x); - my $phi = acos($z); + my $phi = acos_real($z); return ($theta, $phi); } @@ -203,7 +221,8 @@ sub great_circle_destination { my $lat0 = pip2 - $phi0; - my $phi1 = asin(sin($lat0)*cos($dst)+cos($lat0)*sin($dst)*cos($dir0)); + my $phi1 = asin_real(sin($lat0)*cos($dst) + + cos($lat0)*sin($dst)*cos($dir0)); my $theta1 = $theta0 + atan2(sin($dir0)*sin($dst)*cos($lat0), cos($dst)-sin($lat0)*sin($phi1)); @@ -538,7 +557,7 @@ optional, it defaults to 1 (the unit sphere), therefore the distance defaults to radians. If you think geographically the I are longitudes: zero at the -Greenwhich meridian, eastward positive, westward negative--and the +Greenwhich meridian, eastward positive, westward negative -- and the I are latitudes: zero at the North Pole, northward positive, southward negative. B: this formula thinks in mathematics, not geographically: the I zero is at the North Pole, not at the @@ -617,9 +636,11 @@ You can import all the great circle formulas by Notice that the resulting directions might be somewhat surprising if you are looking at a flat worldmap: in such map projections the great -circles quite often do not look like the shortest routes-- but for +circles quite often do not look like the shortest routes -- but for example the shortest possible routes from Europe or North America to -Asia do often cross the polar regions. +Asia do often cross the polar regions. (The common Mercator projection +does B show great circles as straight lines: straight lines in the +Mercator projection are lines of constant bearing.) =head1 EXAMPLES @@ -655,6 +676,51 @@ The answers may be off by few percentages because of the irregular (slightly aspherical) form of the Earth. The errors are at worst about 0.55%, but generally below 0.3%. +=head2 Real-valued asin and acos + +For small inputs asin() and acos() may return complex numbers even +when real numbers would be enough and correct, this happens because of +floating-point inaccuracies. You can see these inaccuracies for +example by trying theses: + + print cos(1e-6)**2+sin(1e-6)**2 - 1,"\n"; + printf "%.20f", cos(1e-6)**2+sin(1e-6)**2,"\n"; + +which will print something like this + + -1.11022302462516e-16 + 0.99999999999999988898 + +even though the expected results are of course exactly zero and one. +The formulas used to compute asin() and acos() are quite sensitive to +this, and therefore they might accidentally slip into the complex +plane even when they should not. To counter this there are two +interfaces that are guaranteed to return a real-valued output. + +=over 4 + +=item asin_real + + use Math::Trig qw(asin_real); + + $real_angle = asin_real($input_sin); + +Return a real-valued arcus sine if the input is between [-1, 1], +B the endpoints. For inputs greater than one, pi/2 +is returned. For inputs less than minus one, -pi/2 is returned. + +=item acos_real + + use Math::Trig qw(acos_real); + + $real_angle = acos_real($input_cos); + +Return a real-valued arcus cosine if the input is between [-1, 1], +B the endpoints. For inputs greater than one, zero +is returned. For inputs less than minus one, pi is returned. + +=back + =head1 BUGS Saying C exports many mathematical routines in the diff --git a/lib/Math/Trig.t b/lib/Math/Trig.t index 9febacc..7a88fe9 100755 --- a/lib/Math/Trig.t +++ b/lib/Math/Trig.t @@ -26,10 +26,10 @@ BEGIN { } } -plan(tests => 135); +plan(tests => 153); -use Math::Trig 1.09; -use Math::Trig 1.09 qw(Inf); +use Math::Trig 1.12; +use Math::Trig 1.12 qw(:pi Inf); my $pip2 = pi / 2; @@ -363,4 +363,30 @@ cmp_ok(csch(-1e5), '==', 0); cmp_ok(tanh(-1e5), '==', -1); cmp_ok(coth(-1e5), '==', -1); +print "# great_circle_distance with small angles\n"; + +for my $e (qw(1e-2 1e-3 1e-4 1e-5)) { + # Can't assume == 0 because of floating point fuzz, + # but let's hope for at least < $e. + cmp_ok(great_circle_distance(0, $e, 0, $e), '<', $e); +} + +print "# asin_real, acos_real\n"; + +is(acos_real(-2.0), pi); +is(acos_real(-1.0), pi); +is(acos_real(-0.5), acos(-0.5)); +is(acos_real( 0.0), acos( 0.0)); +is(acos_real( 0.5), acos( 0.5)); +is(acos_real( 1.0), 0); +is(acos_real( 2.0), 0); + +is(asin_real(-2.0), -&pip2); +is(asin_real(-1.0), -&pip2); +is(asin_real(-0.5), asin(-0.5)); +is(asin_real( 0.0), asin( 0.0)); +is(asin_real( 0.5), asin( 0.5)); +is(asin_real( 1.0), pip2); +is(asin_real( 2.0), pip2); + # eof