Math-Complex: improve numerical accuracy

- Compute the great circle distance with a formula that has better
  numerical properties. The formula used now is accurate for all
  distances.

- Add a few tests to verify.

- This fixes CPAN RT #78938

Improve documentation for great_circle_distance()

For: https://github.com/Perl/perl5/pull/20212

2 files changed, 33 insertions, 25 deletions

diff --git a/dist/Math-Complex/lib/Math/Trig.pm b/dist/Math-Complex/lib/Math/Trig.pm
index 4bd90fbd67..bb96b5e69d 100644
--- a/dist/Math-Complex/lib/Math/Trig.pm
+++ b/dist/Math-Complex/lib/Math/Trig.pm
@@ -154,12 +154,19 @@ sub great_circle_distance {
 
     $rho = 1 unless defined $rho; # Default to the unit sphere.
 
-    my $lat0 = pip2 - $phi0;
-    my $lat1 = pip2 - $phi1;
-
-    return $rho *
-	acos_real( cos( $lat0 ) * cos( $lat1 ) * cos( $theta0 - $theta1 ) +
-		   sin( $lat0 ) * sin( $lat1 ) );
+    my $dphi   = $phi1 - $phi0;
+    my $dtheta = $theta1 - $theta0;
+
+    # A formula that is accurate for all distances is the following special
+    # case of the Vincenty formula for an ellipsoid with equal major and minor
+    # axes.  See
+    # https://en.wikipedia.org/wiki/Great-circle_distance#Computational_formulas
+
+    my $c1 = sin($phi1) * sin($dtheta);
+    my $c2 = sin($phi1) * cos($dtheta);
+    my $c3 = sin($phi0) * cos($phi1) - cos($phi0) * $c2;
+    my $c4 = cos($phi0) * cos($phi1) + sin($phi0) * $c2;
+    return $rho * atan2(sqrt($c1 * $c1 + $c3 * $c3), $c4);
 }
 
 sub great_circle_direction {
@@ -538,26 +545,19 @@ points.
 
 =head2 great_circle_distance
 
-You can compute spherical distances, called B<great circle distances>,
-by importing the great_circle_distance() function:
-
-  use Math::Trig 'great_circle_distance';
+Returns the great circle distance between two points on a sphere.
 
-  $distance = great_circle_distance($theta0, $phi0, $theta1, $phi1, [, $rho]);
+    $distance = great_circle_distance($theta0, $phi0, $theta1, $phi1, [, $rho]);
 
-The I<great circle distance> is the shortest distance between two
-points on a sphere.  The distance is in C<$rho> units.  The C<$rho> is
-optional, it defaults to 1 (the unit sphere), therefore the distance
-defaults to radians.
+Where ($theta0, $phi0) and ($theta1, $phi1) are the spherical coordinates of
+the two points, respectively. The distance is in C<$rho> units. The C<$rho>
+is optional. It defaults to 1 (the unit sphere).
 
-If you think geographically the I<theta> are longitudes: zero at the
-Greenwhich meridian, eastward positive, westward negative -- and the
-I<phi> are latitudes: zero at the North Pole, northward positive,
-southward negative.  B<NOTE>: this formula thinks in mathematics, not
-geographically: the I<phi> zero is at the North Pole, not at the
-Equator on the west coast of Africa (Bay of Guinea).  You need to
-subtract your geographical coordinates from I<pi/2> (also known as 90
-degrees).
+If you are using geographic coordinates, latitude and longitude, you need to
+adjust for the fact that latitude is zero at the equator increasing towards
+the north and decreasing towards the south. Assuming ($lat0, $lon0) and
+($lat1, $lon1) are the geographic coordinates in radians of the two points,
+the distance can be computed with
 
   $distance = great_circle_distance($lon0, pi/2 - $lat0,
                                     $lon1, pi/2 - $lat1, $rho);
diff --git a/dist/Math-Complex/t/Trig.t b/dist/Math-Complex/t/Trig.t
index a9a12556b6..8e2e5d1645 100644
--- a/dist/Math-Complex/t/Trig.t
+++ b/dist/Math-Complex/t/Trig.t
@@ -10,7 +10,7 @@
 
 use strict;
 use warnings;
-use Test::More tests => 153;
+use Test::More tests => 157;
 
 use Math::Trig 1.18;
 use Math::Trig 1.18 qw(:pi Inf);
@@ -363,7 +363,15 @@ 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);
+    cmp_ok(great_circle_distance(0, $e, 0, $e), '<', $e,
+           "great_circle_distance(0, $e, 0, $e) < $e");
+}
+
+for my $e (qw(1e-5 1e-6 1e-7 1e-8)) {
+    # Verify that the distance is positive for points close together. A poor
+    # algorithm is likely to give a distance of zero in some of these cases.
+    cmp_ok(great_circle_distance(2, 2, 2, 2+$e), '>', 0,
+           "great_circle_distance(2, 2, 2, " . (2+$e) . ") > 0");
 }
 
 print "# asin_real, acos_real\n";