Merge pull request #2576 from cffk/geod-1.52

Update geodesic routines from GeographicLib 1.52

4 files changed, 39 insertions, 21 deletions

diff --git a/docs/source/geodesic.rst b/docs/source/geodesic.rst
index 11bb5e5c..6d07352c 100644
--- a/docs/source/geodesic.rst
+++ b/docs/source/geodesic.rst
@@ -47,7 +47,7 @@ Traditionally two geodesic problems are considered:
   :math:`\alpha_2`.
 
 PROJ incorporates `C library for Geodesics
-<https://geographiclib.sourceforge.io/1.49/C/>`_ from `GeographicLib
+<https://geographiclib.sourceforge.io/1.52/C/>`_ from `GeographicLib
 <https://geographiclib.sourceforge.io>`_.  This library provides
 routines to solve the direct and inverse geodesic problems.  Full double
 precision accuracy is maintained provided that
@@ -55,7 +55,7 @@ precision accuracy is maintained provided that
 to the
 
     `application programming interface
-    <https://geographiclib.sourceforge.io/1.49/C/geodesic_8h.html>`_
+    <https://geographiclib.sourceforge.io/1.52/C/geodesic_8h.html>`_
 
 for full documentation.  A brief summary of the routines is given in
 geodesic(3).
diff --git a/src/geodesic.c b/src/geodesic.c
index 53ec9ed6..d5751227 100644
--- a/src/geodesic.c
+++ b/src/geodesic.c
@@ -18,7 +18,7 @@
  *
  * See the comments in geodesic.h for documentation.
  *
- * Copyright (c) Charles Karney (2012-2020) <charles@karney.com> and licensed
+ * Copyright (c) Charles Karney (2012-2021) <charles@karney.com> and licensed
  * under the MIT/X11 License.  For more information, see
  * https://geographiclib.sourceforge.io/
  */
@@ -127,7 +127,13 @@ static void swapx(real* x, real* y)
 { real t = *x; *x = *y; *y = t; }
 
 static void norm2(real* sinx, real* cosx) {
+#if defined(_MSC_VER) && _MSC_VER <= 1900
+  /* Visual Studio 2015 (32-bit) has inaccurate hypot, the same as in some
+   * versions of python https://bugs.python.org/issue43088 */
+  real r = sqrt(*sinx * *sinx + *cosx * *cosx);
+#else
   real r = hypot(*sinx, *cosx);
+#endif
   *sinx /= r;
   *cosx /= r;
 }
@@ -164,23 +170,12 @@ static real AngRound(real x) {
 static void sincosdx(real x, real* sinx, real* cosx) {
   /* In order to minimize round-off errors, this function exactly reduces
    * the argument to the range [-45, 45] before converting it to radians. */
-  real r, s, c; int q;
+  real r, s, c; int q = 0;
   r = remquo(x, (real)(90), &q);
   /* now abs(r) <= 45 */
   r *= degree;
   /* Possibly could call the gnu extension sincos */
   s = sin(r); c = cos(r);
-#if defined(_MSC_VER) && _MSC_VER < 1900
-  /*
-   * Before version 14 (2015), Visual Studio had problems dealing
-   * with -0.0.  Specifically
-   *   VC 10,11,12 and 32-bit compile: fmod(-0.0, 360.0) -> +0.0
-   *   VC 12       and 64-bit compile:  sin(-0.0)        -> +0.0
-   * AngNormalize has a similar fix.
-   * python 2.7 on Windows 32-bit machines has the same problem.
-   */
-  if (x == 0) s = x;
-#endif
   switch ((unsigned)q & 3U) {
   case 0U: *sinx =  s; *cosx =  c; break;
   case 1U: *sinx =  c; *cosx = -s; break;
@@ -797,7 +792,9 @@ static real geod_geninverse_int(const struct geod_geodesic* g,
      * not a shortest path. */
     if (sig12 < 1 || m12x >= 0) {
       /* Need at least 2, to handle 90 0 90 180 */
-      if (sig12 < 3 * tiny)
+      if (sig12 < 3 * tiny ||
+          // Prevent negative s12 or m12 for short lines
+          (sig12 < tol0 && (s12x < 0 || m12x < 0)))
         sig12 = m12x = s12x = 0;
       m12x *= g->b;
       s12x *= g->b;
@@ -868,7 +865,6 @@ static real geod_geninverse_int(const struct geod_geodesic* g,
                         slam12, clam12,
                         &salp2, &calp2, &sig12, &ssig1, &csig1, &ssig2, &csig2,
                         &eps, &domg12, numit < maxit1, &dv, Ca);
-        /* 2 * tol0 is approximately 1 ulp for a number in [0, pi]. */
         /* Reversed test to allow escape with NaNs */
         if (tripb || !(fabs(v) >= (tripn ? 8 : 1) * tol0)) break;
         /* Update bracketing values */
diff --git a/src/geodesic.h b/src/geodesic.h
index 548fbf47..b803d2a5 100644
--- a/src/geodesic.h
+++ b/src/geodesic.h
@@ -107,12 +107,12 @@
  * twice about restructuring the internals of the C code since this may make
  * porting fixes from the C++ code more difficult.
  *
- * Copyright (c) Charles Karney (2012-2020) <charles@karney.com> and licensed
+ * Copyright (c) Charles Karney (2012-2021) <charles@karney.com> and licensed
  * under the MIT/X11 License.  For more information, see
  * https://geographiclib.sourceforge.io/
  *
  * This library was distributed with
- * <a href="../index.html">GeographicLib</a> 1.51.
+ * <a href="../index.html">GeographicLib</a> 1.52.
  **********************************************************************/
 
 #if !defined(GEODESIC_H)
@@ -127,7 +127,7 @@
  * The minor version of the geodesic library.  (This tracks the version of
  * GeographicLib.)
  **********************************************************************/
-#define GEODESIC_VERSION_MINOR 51
+#define GEODESIC_VERSION_MINOR 52
 /**
  * The patch level of the geodesic library.  (This tracks the version of
  * GeographicLib.)
diff --git a/src/tests/geodtest.cpp b/src/tests/geodtest.cpp
index 1b1e7aa7..c4d8b098 100644
--- a/src/tests/geodtest.cpp
+++ b/src/tests/geodtest.cpp
@@ -4,7 +4,7 @@
  *
  * Run these tests by configuring with cmake and running "make test".
  *
- * Copyright (c) Charles Karney (2015-2019) <charles@karney.com> and licensed
+ * Copyright (c) Charles Karney (2015-2021) <charles@karney.com> and licensed
  * under the MIT/X11 License.  For more information, see
  * https://geographiclib.sourceforge.io/
  **********************************************************************/
@@ -737,6 +737,9 @@ static int GeodSolve80() {
   result += checkEquals(M12, 1, 1e-15);
   result += checkEquals(M21, 1, 1e-15);
   result += checkEquals(S12, 0, 1e-10);
+  result += 1/a12 > 0 ? 0 : 1;
+  result += 1/s12 > 0 ? 0 : 1;
+  result += 1/m12 > 0 ? 0 : 1;
 
   a12 = geod_geninverse(&g, 90, 0, 90, 180,
                         &s12, &azi1, &azi2, &m12, &M12, &M21, &S12);
@@ -806,6 +809,24 @@ static int GeodSolve84() {
   return result;
 }
 
+static int GeodSolve92() {
+  /* Check fix for inaccurate hypot with python 3.[89].  Problem reported
+   * by agdhruv https://github.com/geopy/geopy/issues/466 ; see
+   * https://bugs.python.org/issue43088 */
+  double azi1, azi2, s12;
+  struct geod_geodesic g;
+  int result = 0;
+  geod_init(&g, wgs84_a, wgs84_f);
+  geod_inverse(&g,
+               37.757540000000006, -122.47018,
+               37.75754,           -122.470177,
+               &s12, &azi1, &azi2);
+  result += checkEquals(azi1, 89.99999923, 1e-7  );
+  result += checkEquals(azi2, 90.00000106, 1e-7  );
+  result += checkEquals(s12,   0.264,      0.5e-3);
+  return result;
+}
+
 static int Planimeter0() {
   /* Check fix for pole-encircling bug found 2011-03-16 */
   double pa[4][2] = {{89, 0}, {89, 90}, {89, 180}, {89, 270}};
@@ -1092,6 +1113,7 @@ int main() {
   if ((i = GeodSolve78())) {++n; printf("GeodSolve78 fail: %d\n", i);}
   if ((i = GeodSolve80())) {++n; printf("GeodSolve80 fail: %d\n", i);}
   if ((i = GeodSolve84())) {++n; printf("GeodSolve84 fail: %d\n", i);}
+  if ((i = GeodSolve92())) {++n; printf("GeodSolve92 fail: %d\n", i);}
   if ((i = Planimeter0())) {++n; printf("Planimeter0 fail: %d\n", i);}
   if ((i = Planimeter5())) {++n; printf("Planimeter5 fail: %d\n", i);}
   if ((i = Planimeter6())) {++n; printf("Planimeter6 fail: %d\n", i);}