Add atanh, copysign, cbrt, remainder, remquo to math.cpp.

The supported C99 math functions provided by math.cpp are thus

   hypot
   log1p
   asinh
   atanh
   copysign
   cbrt
   remainder
   remquo
   round
   lround

Make compiler checks in CMakeLists.txt and configure.ac consistent with
this set.

Make geodesic.c use the math.cpp defined (instead of the internally
defined) versions of

   hypot
   atanh
   copysign
   cbrt

This is keyed off the presence of the PROJ_LIB macro.  I had at one
time

    https://github.com/OSGeo/PROJ/pull/1425

suggested supplying an additional macro PROJ_COMPILATION when compiling
geodesic.c.  However, PROJ_LIB seems to fill the bill OK.

The *next* version of geodesic.c (due out in a few weeks) will also use

   remainder
   remquo

All of this is only of concern for C compilers without C99 support.  So
this may become an historical footnote at some point.

5 files changed, 139 insertions, 46 deletions

diff --git a/CMakeLists.txt b/CMakeLists.txt
index d44d20ce..55722f91 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -139,11 +139,11 @@ check_c_source_compiles("
 #include <math.h>
 int main() {
   int q;
-  return (int)(
-    hypot(3.0, 4.0) + atanh(0.8) + cbrt(8.0) +
-    remquo(100.0, 90.0, &q) +
-    remainder(100.0, 90.0) + copysign(1.0, -0.0) +
-    log1p(0.1) + asinh(0.1)) + isnan(0.0);
+  return (int)(hypot(3.0, 4.0) + log1p(2.0) + asinh(10.0) +
+               atanh(0.8) + cbrt(8.0) + remquo(100.0, 90.0, &q) +
+               remainder(100.0, 90.0) + copysign(1.0, -0.0) +
+               round(3.5)) +
+         (int)(lround(-3.5)) + isnan(0.0);
 }
 " C99_MATH)
 if(C99_MATH)
diff --git a/configure.ac b/configure.ac
index a5edd13c..76e45bc0 100644
--- a/configure.ac
+++ b/configure.ac
@@ -204,11 +204,11 @@ AC_MSG_CHECKING([for C99 math functions])
 AC_LINK_IFELSE([AC_LANG_PROGRAM(
         [#include <math.h>],
         [int q;
-         return (int)(hypot(3.0, 4.0) + atanh(0.8) + cbrt(8.0) +
-                      remquo(100.0, 90.0, &q) +
+         return (int)(hypot(3.0, 4.0) + log1p(2.0) + asinh(10.0) +
+                      atanh(0.8) + cbrt(8.0) + remquo(100.0, 90.0, &q) +
                       remainder(100.0, 90.0) + copysign(1.0, -0.0) +
-                      log1p(0.1) + asinh(0.1)) +
-		      isnan(0.0);])],
+                      round(3.5)) +
+                (int)(lround(-3.5)) + isnan(0.0);])],
         [AC_MSG_RESULT([yes]);C99_MATH="-DHAVE_C99_MATH=1"],
         [AC_MSG_RESULT([no]);C99_MATH="-DHAVE_C99_MATH=0"])
 CFLAGS="$save_CFLAGS $C99_MATH"
diff --git a/src/geodesic.c b/src/geodesic.c
index 5504cb3b..887edb42 100644
--- a/src/geodesic.c
+++ b/src/geodesic.c
@@ -23,8 +23,19 @@
  * https://geographiclib.sourceforge.io/
  */
 
+/* The PROJ_COMPILATION flag indicates that this is part of the compilation of
+ * the PROJ library (keyed off the presence of the PROJ_LIB macro which points
+ * to the data directory for PROJ).  If this is set, we use the PROJ supplied
+ * implementations of the C99 math functions instead of the ones defined here.
+ */
+#if defined(PROJ_LIB)
+#define PROJ_COMPILATION 1
+#else
+#define PROJ_COMPILATION 0
+#endif
+
 #include "geodesic.h"
-#ifdef PJ_LIB__
+#if PROJ_COMPILAION
 #include "proj_math.h"
 #else
 #include <math.h>
@@ -122,7 +133,7 @@ enum captype {
 };
 
 static real sq(real x) { return x * x; }
-#if HAVE_C99_MATH
+#if HAVE_C99_MATH || PROJ_COMPILATION
 #define atanhx atanh
 #define copysignx copysign
 #define hypotx hypot
diff --git a/src/math.cpp b/src/math.cpp
index 540ab9eb..0ec4f57f 100644
--- a/src/math.cpp
+++ b/src/math.cpp
@@ -44,64 +44,137 @@ int pj_isnan (double x) {
 
 #if !(defined(HAVE_C99_MATH) && HAVE_C99_MATH)
 
+/* Define C99 compatible versions of
+ *     hypot
+ *     log1p
+ *     asinh
+ *     atanh
+ *     copysign
+ *     cbrt
+ *     remainder
+ *     remquo
+ *     round
+ *     lround
+ */
+
 /* Compute hypotenuse */
 double pj_hypot(double x, double y) {
-    x = fabs(x);
-    y = fabs(y);
-    if ( x < y ) {
-        x /= y;
-        return ( y * sqrt( 1. + x * x ) );
-    } else {
-        y /= (x != 0.0 ? x : 1.0);
-        return ( x * sqrt( 1. + y * y ) );
-    }
+  x = fabs(x);
+  y = fabs(y);
+  if (x < y) {
+    x /= y;                     /* y is nonzero */
+    return y * sqrt(1 + x * x);
+  } else {
+    y /= (x ? x : 1);
+    return x * sqrt(1 + y * y);
+  }
 }
 
 /* Compute log(1+x) accurately */
 double pj_log1p(double x) {
-    volatile double
-      y = 1 + x,
-      z = y - 1;
-    /* Here's the explanation for this magic: y = 1 + z, exactly, and z
-     * approx x, thus log(y)/z (which is nearly constant near z = 0) returns
-     * a good approximation to the true log(1 + x)/x.  The multiplication x *
-     * (log(y)/z) introduces little additional error. */
-    return z == 0 ? x : x * log(y) / z;
+  volatile double
+    y = 1 + x,
+    z = y - 1;
+  /* Here's the explanation for this magic: y = 1 + z, exactly, and z
+   * approx x, thus log(y)/z (which is nearly constant near z = 0) returns
+   * a good approximation to the true log(1 + x)/x.  The multiplication x *
+   * (log(y)/z) introduces little additional error. */
+  return z == 0 ? x : x * log(y) / z;
 }
 
 /* Compute asinh(x) accurately */
 double pj_asinh(double x) {
-    double y = fabs(x); /* Enforce odd parity */
-    y = log1p(y * (1 + y/(hypot(1.0, y) + 1)));
-    return x > 0 ? y : (x < 0 ? -y : x);
+  double y = fabs(x);           /* Enforce odd parity */
+  y = pj_log1p(y * (1 + y/(pj_hypot(1.0, y) + 1)));
+  return x > 0 ? y : (x < 0 ? -y : x); /* asinh(-0.0) = -0.0 */
+}
+
+/* Compute atanh(x) accurately */
+double pj_atanh(double x) {
+  double y = fabs(x);             /* Enforce odd parity */
+  y = pj_log1p(2 * y/(1 - y))/2;
+  return x > 0 ? y : (x < 0 ? -y : x); /* atanh(-0.0) = -0.0 */
+}
+
+/* Implement copysign(x, y)  */
+double pj_copysign(double x, double y) {
+  /* 1/y trick to get the sign of -0.0 */
+  return fabs(x) * (y < 0 || (y == 0 && 1/y < 0) ? -1 : 1);
+}
+
+/* Implement cbrt(x)  */
+double pj_cbrt(double x) {
+  double y = pow(fabs(x), 1/3.0);      /* Return the real cube root */
+  return x > 0 ? y : (x < 0 ? -y : x); /* cbrt(-0.0) = -0.0 */
+}
+
+/* Implement remainder(x, y) with ties to even */
+double pj_remainder(double x, double y) {
+  double z;
+  y = fabs(y);                 /* The result doesn't depend on the sign of y */
+  z = fmod(x, y);
+  if (z == 0)
+    /* This shouldn't be necessary.  However, 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
+     * python 2.7 on Windows 32-bit machines has the same problem. */
+    z = pj_copysign(z, x);
+  else if (2 * fabs(z) == y)
+    z -= fmod(x, 2 * y) - z;    /* Implement ties to even */
+  else if (2 * fabs(z) > y)
+    z += (z < 0 ? y : -y);      /* Fold remaining cases to (-y/2, y/2) */
+  return z;
+}
+
+/* Implement remquo(x, y, n) with n giving low 3 bits + sign of x/y */
+double pj_remquo(double x, double y, int* n) {
+  double z = pj_remainder(x, y);
+  if (n) {
+    double
+      a = pj_remainder(x, 2 * y),
+      b = pj_remainder(x, 4 * y),
+      c = pj_remainder(x, 8 * y);
+    *n  = (a > z ? 1 : (a < z ? -1 : 0));
+    *n += (b > a ? 2 : (b < a ? -2 : 0));
+    *n += (c > b ? 4 : (c < b ? -4 : 0));
+    if (y < 0) *n *= -1;
+    if (y != 0) {
+      if (x/y > 0 && *n <= 0)
+        *n += 8;
+      else if (x/y < 0 && *n >= 0)
+        *n -= 8;
+    }
+  }
+  return z;
 }
 
+/* Implement round(x) */
 double pj_round(double x) {
   /* The handling of corner cases is copied from boost; see
    *   https://github.com/boostorg/math/pull/8
    * with improvements to return -0.0 when appropriate */
   double t;
-  if (x == 0)
-    return x;               /* Retain sign of 0 */
-  else if (0 < x && x <  0.5)
+  if      (0 < x && x <  0.5)
     return +0.0;
   else if (0 > x && x > -0.5)
     return -0.0;
-  else if (x > 0) {
+  else if   (x > 0) {
     t = ceil(x);
     return 0.5 < t - x ? t - 1 : t;
-  } else {                  /* Includes NaN */
+  } else if (x < 0) {
     t = floor(x);
     return 0.5 < x - t ? t + 1 : t;
-  }
+  } else                        /* +/-0 and NaN */
+    return x;                   /* retain sign of 0 */
 }
 
+/* Implement lround(x) */
 long pj_lround(double x) {
   /* Default value for overflow + NaN + (x == LONG_MIN) */
   long r = LONG_MIN;
-  x = round(x);
-  if (fabs(x) < -(double)LONG_MIN) /* Assume (double)LONG_MIN is exact */
-    r = (int)x;
+  x = pj_round(x);
+  if (fabs(x) < -(double)r)     /* Assume (double)LONG_MIN is exact */
+    r = (long)x;
   return r;
 }
 
diff --git a/src/proj_math.h b/src/proj_math.h
index 5aea494d..414df805 100644
--- a/src/proj_math.h
+++ b/src/proj_math.h
@@ -64,16 +64,25 @@ extern "C" {
 double  pj_hypot(double x, double y);
 double  pj_log1p(double x);
 double  pj_asinh(double x);
+double  pj_atanh(double x);
+double  pj_copysign(double x, double y);
+double  pj_cbrt(double x);
+double  pj_remainder(double x, double y);
+double  pj_remquo(double x, double y, int* n);
 double  pj_round(double x);
 long    pj_lround(double x);
 int     PROJ_DLL pj_isnan(double x);
 
-#define hypot   pj_hypot
-#define log1p   pj_log1p
-#define asinh   pj_asinh
-#define round   pj_round
-#define lround  pj_lround
-
+#define hypot     pj_hypot
+#define log1p     pj_log1p
+#define asinh     pj_asinh
+#define atanh     pj_atanh
+#define copysign  pj_copysign
+#define cbrt      pj_cbrt
+#define remainder pj_remainder
+#define remquo    pj_remquo
+#define round     pj_round
+#define lround    pj_lround
 
 #ifdef isnan
 #undef isnan