Merge pull request #937 from kbevers/c99-math-module

Collect custom C99 math functions in proj_math.h

41 files changed, 179 insertions, 119 deletions

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 6d7c4f77..d7afe060 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -96,7 +96,8 @@ 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)) +
+               remainder(100.0, 90.0) + copysign(1.0, -0.0) +
+               log1p(0.1) + asinh(0.1)) +
                isnan(0.0);
 }\n" C99_MATH)
 if (C99_MATH)
diff --git a/configure.ac b/configure.ac
index 269225a5..e89961eb 100644
--- a/configure.ac
+++ b/configure.ac
@@ -45,7 +45,8 @@ AC_LINK_IFELSE([AC_LANG_PROGRAM(
         [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)) +
+                      remainder(100.0, 90.0) + copysign(1.0, -0.0) +
+                      log1p(0.1) + asinh(0.1)) +
 		      isnan(0.0);])],
         [AC_MSG_RESULT([yes]);C99_MATH="-DHAVE_C99_MATH=1"],
         [AC_MSG_RESULT([no]);C99_MATH="-DHAVE_C99_MATH=0"])
diff --git a/src/Makefile.am b/src/Makefile.am
index f94e5f26..78ecfdf8 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -43,7 +43,7 @@ lib_LTLIBRARIES = libproj.la
 libproj_la_LDFLAGS = -no-undefined -version-info 13:0:0
 
 libproj_la_SOURCES = \
-	pj_list.h proj_internal.h\
+	pj_list.h proj_internal.h proj_math.h\
 	PJ_aeqd.c PJ_gnom.c PJ_laea.c PJ_mod_ster.c \
 	PJ_nsper.c PJ_nzmg.c PJ_ortho.c PJ_stere.c PJ_sterea.c \
 	PJ_aea.c PJ_bipc.c PJ_bonne.c PJ_eqdc.c PJ_isea.c PJ_ccon.c\
@@ -83,7 +83,7 @@ libproj_la_SOURCES = \
 	pj_apply_gridshift.c pj_datums.c pj_datum_set.c pj_transform.c \
 	geocent.c geocent.h pj_utils.c pj_gridinfo.c pj_gridlist.c \
 	jniproj.c pj_mutex.c pj_initcache.c pj_apply_vgridshift.c geodesic.c \
-	pj_strtod.c \
+	pj_strtod.c pj_math.c\
 	\
 	proj_4D_api.c PJ_cart.c PJ_pipeline.c PJ_horner.c PJ_helmert.c \
 	PJ_vgridshift.c PJ_hgridshift.c PJ_unitconvert.c PJ_molodensky.c \
diff --git a/src/PJ_aea.c b/src/PJ_aea.c
index 3fe90524..640f1013 100644
--- a/src/PJ_aea.c
+++ b/src/PJ_aea.c
@@ -31,6 +31,8 @@
 #include "proj.h"
 #include <errno.h>
 #include "projects.h"
+#include "proj_math.h"
+
 
 # define EPS10  1.e-10
 # define TOL7   1.e-7
diff --git a/src/PJ_aeqd.c b/src/PJ_aeqd.c
index 84554b91..5d8c3d38 100644
--- a/src/PJ_aeqd.c
+++ b/src/PJ_aeqd.c
@@ -30,6 +30,7 @@
 #include "proj.h"
 #include <errno.h>
 #include "projects.h"
+#include "proj_math.h"
 
 enum Mode {
     N_POLE = 0,
diff --git a/src/PJ_bipc.c b/src/PJ_bipc.c
index 48fc93d0..2f60808d 100644
--- a/src/PJ_bipc.c
+++ b/src/PJ_bipc.c
@@ -2,6 +2,7 @@
 #include "proj.h"
 #include <errno.h>
 #include "projects.h"
+#include "proj_math.h"
 
 PROJ_HEAD(bipc, "Bipolar conic of western hemisphere") "\n\tConic Sph.";
 
diff --git a/src/PJ_bonne.c b/src/PJ_bonne.c
index 40bb28ee..d3d5b757 100644
--- a/src/PJ_bonne.c
+++ b/src/PJ_bonne.c
@@ -2,6 +2,8 @@
 #include <errno.h>
 #include "proj.h"
 #include "projects.h"
+#include "proj_math.h"
+
 
 PROJ_HEAD(bonne, "Bonne (Werner lat_1=90)")
     "\n\tConic Sph&Ell\n\tlat_1=";
diff --git a/src/PJ_cart.c b/src/PJ_cart.c
index ec273f1b..a4fd3254 100644
--- a/src/PJ_cart.c
+++ b/src/PJ_cart.c
@@ -43,8 +43,8 @@
 #define PJ_LIB__
 #include "proj_internal.h"
 #include "projects.h"
+#include "proj_math.h"
 #include <stddef.h>
-#include <math.h>
 #include <errno.h>
 PROJ_HEAD(cart,    "Geodetic/cartesian conversions");
 
diff --git a/src/PJ_ccon.c b/src/PJ_ccon.c
index 02b14388..aed59fb6 100644
--- a/src/PJ_ccon.c
+++ b/src/PJ_ccon.c
@@ -24,6 +24,7 @@
 #include <errno.h>
 #include "proj.h"
 #include "projects.h"
+#include "proj_math.h"
 
 #define EPS10   1e-10
 
diff --git a/src/PJ_deformation.c b/src/PJ_deformation.c
index 58c0e4df..5511eed4 100644
--- a/src/PJ_deformation.c
+++ b/src/PJ_deformation.c
@@ -55,6 +55,7 @@ grid-values in units of mm/year in ENU-space.
 #include <errno.h>
 #include "proj.h"
 #include "proj_internal.h"
+#include "proj_math.h"
 #include "projects.h"
 
 PROJ_HEAD(deformation, "Kinematic grid shift");
diff --git a/src/PJ_eqdc.c b/src/PJ_eqdc.c
index f36e0518..5f65e7d0 100644
--- a/src/PJ_eqdc.c
+++ b/src/PJ_eqdc.c
@@ -2,6 +2,7 @@
 #include <errno.h>
 #include "proj.h"
 #include "projects.h"
+#include "proj_math.h"
 
 struct pj_opaque {
     double phi1;
diff --git a/src/PJ_geos.c b/src/PJ_geos.c
index f0a9bc24..2988e62e 100644
--- a/src/PJ_geos.c
+++ b/src/PJ_geos.c
@@ -31,6 +31,7 @@
 #include <errno.h>
 #include "proj.h"
 #include "projects.h"
+#include "proj_math.h"
 
 struct pj_opaque {
     double h;
diff --git a/src/PJ_gnom.c b/src/PJ_gnom.c
index 8bb24b6a..b4bab0e2 100644
--- a/src/PJ_gnom.c
+++ b/src/PJ_gnom.c
@@ -2,6 +2,7 @@
 #include <errno.h>
 #include "proj.h"
 #include "projects.h"
+#include "proj_math.h"
 
 PROJ_HEAD(gnom, "Gnomonic") "\n\tAzi, Sph.";
 
diff --git a/src/PJ_laea.c b/src/PJ_laea.c
index 0422784b..bcf9c44d 100644
--- a/src/PJ_laea.c
+++ b/src/PJ_laea.c
@@ -2,6 +2,7 @@
 #include <errno.h>
 #include "proj.h"
 #include "projects.h"
+#include "proj_math.h"
 
 PROJ_HEAD(laea, "Lambert Azimuthal Equal Area") "\n\tAzi, Sph&Ell";
 
diff --git a/src/PJ_lcc.c b/src/PJ_lcc.c
index 863f6be9..34ed99cb 100644
--- a/src/PJ_lcc.c
+++ b/src/PJ_lcc.c
@@ -2,6 +2,7 @@
 #include <errno.h>
 #include "proj.h"
 #include "projects.h"
+#include "proj_math.h"
 
 PROJ_HEAD(lcc, "Lambert Conformal Conic")
     "\n\tConic, Sph&Ell\n\tlat_1= and lat_2= or lat_0";
diff --git a/src/PJ_merc.c b/src/PJ_merc.c
index 0bf98625..9123eae1 100644
--- a/src/PJ_merc.c
+++ b/src/PJ_merc.c
@@ -1,6 +1,7 @@
 #define PJ_LIB__
 #include "proj_internal.h"
 #include "proj.h"
+#include "proj_math.h"
 #include "projects.h"
 #include <float.h>
 
@@ -8,30 +9,10 @@ PROJ_HEAD(merc, "Mercator") "\n\tCyl, Sph&Ell\n\tlat_ts=";
 PROJ_HEAD(webmerc, "Web Mercator / Pseudo Mercator") "\n\tCyl, Sph\n\t";
 
 #define EPS10 1.e-10
-
-#if !defined(HAVE_C99_MATH)
-#define HAVE_C99_MATH 0
-#endif
-
-#if HAVE_C99_MATH
-#define log1px log1p
-#else
-static double log1px(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;
-}
-#endif
-
 static double logtanpfpim1(double x) {       /* log(tan(x/2 + M_FORTPI)) */
     if (fabs(x) <= DBL_EPSILON) {
         /* tan(M_FORTPI + .5 * x) can be approximated by  1.0 + x */
-        return log1px(x);
+        return log1p(x);
     }
     return log(tan(M_FORTPI + .5 * x));
 }
diff --git a/src/PJ_mod_ster.c b/src/PJ_mod_ster.c
index d807660c..5e6ce136 100644
--- a/src/PJ_mod_ster.c
+++ b/src/PJ_mod_ster.c
@@ -2,6 +2,7 @@
 #define PJ_LIB__
 #include <errno.h>
 #include "projects.h"
+#include "proj_math.h"
 
 PROJ_HEAD(mil_os, "Miller Oblated Stereographic") "\n\tAzi(mod)";
 PROJ_HEAD(lee_os, "Lee Oblated Stereographic") "\n\tAzi(mod)";
diff --git a/src/PJ_nsper.c b/src/PJ_nsper.c
index 3c0b01f0..223cd75b 100644
--- a/src/PJ_nsper.c
+++ b/src/PJ_nsper.c
@@ -2,6 +2,7 @@
 #include <errno.h>
 #include "proj.h"
 #include "projects.h"
+#include "proj_math.h"
 
 enum Mode {
     N_POLE = 0,
diff --git a/src/PJ_oea.c b/src/PJ_oea.c
index 5112896f..2bfeffa8 100644
--- a/src/PJ_oea.c
+++ b/src/PJ_oea.c
@@ -2,6 +2,7 @@
 #include <errno.h>
 #include "proj.h"
 #include "projects.h"
+#include "proj_math.h"
 
 PROJ_HEAD(oea, "Oblated Equal Area") "\n\tMisc Sph\n\tn= m= theta=";
 
diff --git a/src/PJ_ortho.c b/src/PJ_ortho.c
index 7c27cd56..2b037819 100644
--- a/src/PJ_ortho.c
+++ b/src/PJ_ortho.c
@@ -2,6 +2,7 @@
 #include <errno.h>
 #include "proj.h"
 #include "proj_internal.h"
+#include "proj_math.h"
 #include "projects.h"
 
 PROJ_HEAD(ortho, "Orthographic") "\n\tAzi, Sph.";
diff --git a/src/PJ_robin.c b/src/PJ_robin.c
index 92aebfd3..a1d5ad4d 100644
--- a/src/PJ_robin.c
+++ b/src/PJ_robin.c
@@ -1,4 +1,5 @@
 #define PJ_LIB__
+#include "proj_math.h"
 #include "proj_internal.h"
 #include "proj.h"
 #include "projects.h"
@@ -82,7 +83,7 @@ static XY s_forward (LP lp, PJ *P) {           /* Spheroidal, forward */
     (void) P;
 
     dphi = fabs(lp.phi);
-    i = pj_is_nan(lp.phi) ? -1 : (int)floor(dphi * C1);
+    i = isnan(lp.phi) ? -1 : (int)floor(dphi * C1);
     if( i < 0 ){
         proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION);
         return xy;
@@ -117,7 +118,7 @@ static LP s_inverse (XY xy, PJ *P) {           /* Spheroidal, inverse */
         }
     } else { /* general problem */
         /* in Y space, reduce to table interval */
-        i = pj_is_nan(lp.phi) ? -1 : (int)floor(lp.phi * NODES);
+        i = isnan(lp.phi) ? -1 : (int)floor(lp.phi * NODES);
         if( i < 0 || i >= NODES ) {
             proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION);
             return lp;
diff --git a/src/PJ_sconics.c b/src/PJ_sconics.c
index 34a0f1fd..ce044c24 100644
--- a/src/PJ_sconics.c
+++ b/src/PJ_sconics.c
@@ -2,6 +2,7 @@
 #include <errno.h>
 #include "proj.h"
 #include "projects.h"
+#include "proj_math.h"
 
 
 enum Type {
diff --git a/src/PJ_stere.c b/src/PJ_stere.c
index 62e42460..82fd5c07 100644
--- a/src/PJ_stere.c
+++ b/src/PJ_stere.c
@@ -2,6 +2,7 @@
 #include <errno.h>
 #include "proj.h"
 #include "projects.h"
+#include "proj_math.h"
 
 PROJ_HEAD(stere, "Stereographic") "\n\tAzi, Sph&Ell\n\tlat_ts=";
 PROJ_HEAD(ups, "Universal Polar Stereographic") "\n\tAzi, Sph&Ell\n\tsouth";
diff --git a/src/PJ_sterea.c b/src/PJ_sterea.c
index 96b35609..eb4c9f2c 100644
--- a/src/PJ_sterea.c
+++ b/src/PJ_sterea.c
@@ -26,6 +26,7 @@
 #define PJ_LIB__
 #include <errno.h>
 #include "projects.h"
+#include "proj_math.h"
 
 
 struct pj_opaque {
diff --git a/src/PJ_tpeqd.c b/src/PJ_tpeqd.c
index 21e3c5ed..87877ec1 100644
--- a/src/PJ_tpeqd.c
+++ b/src/PJ_tpeqd.c
@@ -1,6 +1,7 @@
 #define PJ_LIB__
 #include <errno.h>
 #include "proj.h"
+#include "proj_math.h"
 #include "projects.h"
 
 
diff --git a/src/geodesic.c b/src/geodesic.c
index eb956c65..9904c7fa 100644
--- a/src/geodesic.c
+++ b/src/geodesic.c
@@ -24,7 +24,11 @@
  */
 
 #include "geodesic.h"
+#ifdef PJ_LIB__
+#include "proj_math.h"
+#else
 #include <math.h>
+#endif
 
 #if !defined(HAVE_C99_MATH)
 #define HAVE_C99_MATH 0
@@ -239,8 +243,7 @@ static void sincosdx(real x, real* sinx, real* cosx) {
   r = remquo(x, (real)(90), &q);
 #else
   r = fmod(x, (real)(360));
-  /* check for NaN -- do not use pj_is_nan, since we want geodesic.c not to
-   * depend on the rest of proj.4 */
+  /* check for NaN */
   q = r == r ? (int)(floor(r / 90 + (real)(0.5))) : 0;
   r -= 90 * q;
 #endif
diff --git a/src/gie.c b/src/gie.c
index 3831a8f1..b81e6bdc 100644
--- a/src/gie.c
+++ b/src/gie.c
@@ -114,6 +114,7 @@ Thomas Knudsen, thokn@sdfe.dk, 2017-10-01/2017-10-08
 
 #include "proj.h"
 #include "proj_internal.h"
+#include "proj_math.h"
 #include "projects.h"
 
 #include "optargpm.h"
diff --git a/src/hypot.c b/src/hypot.c
deleted file mode 100644
index 822c4595..00000000
--- a/src/hypot.c
+++ /dev/null
@@ -1,36 +0,0 @@
-/* hypot - sqrt(x * x + y * y)
-**
-**	Because this was omitted from the ANSI standards, this version
-** 	is included for those systems that do not include hypot as an
-**	extension to libm.a.  Note: GNU version was not used because it
-**	was not properly coded to minimize potential overflow.
-**
-**	The proper technique for determining hypot is to factor out the
-**	larger of the two terms,  thus leaving a possible case of float
-**	overflow when max(x,y)*sqrt(2) > max machine value.  This allows
-**	a wider range of numbers than the alternative of the sum of the
-**	squares < max machine value.  For an Intel x87 IEEE double of
-**	approximately 1.8e308, only argument values > 1.27e308 are at
-**	risk of causing overflow.  Whereas, not using this method limits
-**	the range to values less that 9.5e153 --- a considerable reduction
-**	in range!
-*/
-extern double sqrt(double);
-	double
-hypot(double x, double y) {
-	if ( x < 0.)
-		x = -x;
-	else if (x == 0.)
-		return (y < 0. ? -y : y);
-	if (y < 0.)
-		y = -y;
-	else if (y == 0.)
-		return (x);
-	if ( x < y ) {
-		x /= y;
-		return ( y * sqrt( 1. + x * x ) );
-	} else {
-		y /= x;
-		return ( x * sqrt( 1. + y * y ) );
-	}
-}
diff --git a/src/lib_proj.cmake b/src/lib_proj.cmake
index c9e4d9e6..c0d0a6a9 100644
--- a/src/lib_proj.cmake
+++ b/src/lib_proj.cmake
@@ -199,6 +199,7 @@ SET(SRC_LIBPROJ_CORE
         pj_list.h
         pj_log.c
         pj_malloc.c
+        pj_math.c
         pj_mlfn.c
         pj_msfn.c
         pj_mutex.c
@@ -218,6 +219,7 @@ SET(SRC_LIBPROJ_CORE
         pj_utils.c
         pj_zpoly1.c
         proj_mdist.c
+        proj_math.h
         proj_rouss.c
         rtodms.c
         vector1.c
diff --git a/src/makefile.vc b/src/makefile.vc
index ee89beb3..ef39b084 100644
--- a/src/makefile.vc
+++ b/src/makefile.vc
@@ -55,7 +55,7 @@ support = \
 	pj_utils.obj pj_gridlist.obj pj_gridinfo.obj \
 	proj_mdist.obj pj_mutex.obj pj_initcache.obj \
 	pj_ctx.obj pj_fileapi.obj pj_log.obj pj_apply_vgridshift.obj \
-	pj_strtod.obj pj_internal.obj
+	pj_strtod.obj pj_internal.obj pj_math.obj
 
 pipeline = \
 	proj_4D_api.obj PJ_cart.obj PJ_pipeline.obj PJ_horner.obj PJ_helmert.obj \
diff --git a/src/nad_cvt.c b/src/nad_cvt.c
index c913511f..d4330d54 100644
--- a/src/nad_cvt.c
+++ b/src/nad_cvt.c
@@ -1,5 +1,6 @@
 #define PJ_LIB__
 #include "projects.h"
+#include "proj_math.h"
 #define MAX_ITERATIONS 10
 #define TOL 1e-12
 
diff --git a/src/nad_intr.c b/src/nad_intr.c
index dc245831..80e428ed 100644
--- a/src/nad_intr.c
+++ b/src/nad_intr.c
@@ -1,6 +1,7 @@
 /* Determine nad table correction value */
 #define PJ_LIB__
 #include "proj_internal.h"
+#include "proj_math.h"
 #include "projects.h"
 
 	LP
@@ -13,9 +14,9 @@ nad_intr(LP t, struct CTABLE *ct) {
 	int in;
 
 	t.lam /= ct->del.lam;
-	indx.lam = pj_is_nan(t.lam) ? 0 : (int)floor(t.lam);
+	indx.lam = isnan(t.lam) ? 0 : (int)floor(t.lam);
 	t.phi /= ct->del.phi;
-	indx.phi = pj_is_nan(t.phi) ? 0 : (int)floor(t.phi);
+	indx.phi = isnan(t.phi) ? 0 : (int)floor(t.phi);
 
 	frct.lam = t.lam - indx.lam;
 	frct.phi = t.phi - indx.phi;
diff --git a/src/pj_apply_vgridshift.c b/src/pj_apply_vgridshift.c
index 05803bc8..e7106b25 100644
--- a/src/pj_apply_vgridshift.c
+++ b/src/pj_apply_vgridshift.c
@@ -29,7 +29,7 @@
 #define PJ_LIB__
 
 #include <string.h>
-#include <math.h>
+#include "proj_math.h"
 #include "proj_internal.h"
 #include "projects.h"
 
@@ -42,7 +42,7 @@ static double read_vgrid_value( PJ *defn, LP input, int *gridlist_count_p, PJ_GR
     float  *cvs;
     /* do not deal with NaN coordinates */
     /* cppcheck-suppress duplicateExpression */
-    if( pj_is_nan(input.phi) || pj_is_nan(input.lam) )
+    if( isnan(input.phi) || isnan(input.lam) )
         itable = *gridlist_count_p;
 
     /* keep trying till we find a table that works */
diff --git a/src/pj_factors.c b/src/pj_factors.c
index e682a12d..373bd743 100644
--- a/src/pj_factors.c
+++ b/src/pj_factors.c
@@ -1,6 +1,7 @@
 /* projection scale factors */
 #define PJ_LIB__
 #include "proj.h"
+#include "proj_math.h"
 #include "projects.h"
 
 #include <errno.h>
diff --git a/src/pj_init.c b/src/pj_init.c
index 9b06ff28..bdaf64f7 100644
--- a/src/pj_init.c
+++ b/src/pj_init.c
@@ -37,6 +37,7 @@
 #include <ctype.h>
 #include "proj.h"
 #include "proj_internal.h"
+#include "proj_math.h"
 #include "projects.h"
 
 
diff --git a/src/pj_internal.c b/src/pj_internal.c
index 891e0b9f..dc528649 100644
--- a/src/pj_internal.c
+++ b/src/pj_internal.c
@@ -438,19 +438,3 @@ void proj_log_func (PJ_CONTEXT *ctx, void *app_data, PJ_LOG_FUNCTION logf) {
     if (0!=logf)
         ctx->logger = logf;
 }
-
-
-#if HAVE_C99_MATH
-/* proj_internal.h defines pj_is_nan as isnan */
-#else
-/*****************************************************************************/
-int pj_is_nan (double val) {
-/******************************************************************************
-    Returns 0 if not a NaN and non-zero if val is a NaN.
-
-    Provides an equivalent to isnan().
-******************************************************************************/
-    /* cppcheck-suppress duplicateExpression */
-    return val != val;
-}
-#endif
diff --git a/src/pj_math.c b/src/pj_math.c
new file mode 100644
index 00000000..751c2750
--- /dev/null
+++ b/src/pj_math.c
@@ -0,0 +1,71 @@
+/******************************************************************************
+ * Project:  PROJ
+ * Purpose:  Make C99 math functions available on C89 systems
+ * Author:   Kristian Evers
+ *
+ ******************************************************************************
+ * Copyright (c) 2018, Kristian Evers
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the "Software"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+ * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+ * DEALINGS IN THE SOFTWARE.
+ *****************************************************************************/
+
+#include "proj_math.h"
+
+#if !(defined(HAVE_C99_MATH) && HAVE_C99_MATH)
+
+/* 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 ) );
+    }
+}
+
+/* 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;
+}
+
+/* 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);
+}
+
+/* Returns 0 if not a NaN and non-zero if val is a NaN */
+int pj_isnan (double x) {
+    /* cppcheck-suppress duplicateExpression */
+    return x != x;
+}
+
+
+#endif /* !(defined(HAVE_C99_MATH) && HAVE_C99_MATH) */
diff --git a/src/proj_4D_api.c b/src/proj_4D_api.c
index aed4685b..70110746 100644
--- a/src/proj_4D_api.c
+++ b/src/proj_4D_api.c
@@ -30,6 +30,7 @@
 #include <errno.h>
 #include "proj.h"
 #include "proj_internal.h"
+#include "proj_math.h"
 #include "projects.h"
 #include "geodesic.h"
 
diff --git a/src/proj_etmerc.c b/src/proj_etmerc.c
index 99646d14..9a89e2b0 100644
--- a/src/proj_etmerc.c
+++ b/src/proj_etmerc.c
@@ -44,6 +44,7 @@
 #include <errno.h>
 #include "proj.h"
 #include "projects.h"
+#include "proj_math.h"
 
 
 struct pj_opaque {
@@ -62,32 +63,6 @@ PROJ_HEAD(utm, "Universal Transverse Mercator (UTM)")
 
 #define PROJ_ETMERC_ORDER 6
 
-
-#ifdef _GNU_SOURCE
-    inline
-#endif
-static double log1py(double x) {              /* Compute log(1+x) accurately */
-    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;
-}
-
-
-#ifdef _GNU_SOURCE
-    inline
-#endif
-static double asinhy(double x) {              /* Compute asinh(x) accurately */
-    double y = fabs(x);         /* Enforce odd parity */
-    y = log1py(y * (1 + y/(hypot(1.0, y) + 1)));
-    return x < 0 ? -y : y;
-}
-
-
 #ifdef _GNU_SOURCE
     inline
 #endif
@@ -182,7 +157,7 @@ static XY e_forward (LP lp, PJ *P) {          /* Ellipsoidal, forward */
     Ce     = atan2 (sin_Ce*cos_Cn,  hypot (sin_Cn, cos_Cn*cos_Ce));
 
     /* compl. sph. N, E -> ell. norm. N, E */
-    Ce  = asinhy ( tan (Ce) );     /* Replaces: Ce  = log(tan(FORTPI + Ce*0.5)); */
+    Ce  = asinh ( tan (Ce) );     /* Replaces: Ce  = log(tan(FORTPI + Ce*0.5)); */
     Cn += clenS (Q->gtu, PROJ_ETMERC_ORDER, 2*Cn, 2*Ce, &dCn, &dCe);
     Ce += dCe;
     if (fabs (Ce) <= 2.623395162778) {
diff --git a/src/proj_math.h b/src/proj_math.h
new file mode 100644
index 00000000..0108a439
--- /dev/null
+++ b/src/proj_math.h
@@ -0,0 +1,55 @@
+/******************************************************************************
+ * Project:  PROJ
+ * Purpose:  Make C99 math functions available on C89 systems
+ * Author:   Kristian Evers
+ *
+ ******************************************************************************
+ * Copyright (c) 2018, Kristian Evers
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the "Software"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+ * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+ * DEALINGS IN THE SOFTWARE.
+ *****************************************************************************/
+
+#include <math.h>
+
+#ifndef PROJ_MATH_H
+#define PROJ_MATH_H
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+#if !(defined(HAVE_C99_MATH) && HAVE_C99_MATH)
+
+
+double  pj_hypot(double x, double y);
+double  pj_log1p(double x);
+double  pj_asinh(double x);
+int     pj_isnan(double x);
+
+#define hypot   pj_hypot
+#define log1p   pj_log1p
+#define asinh   pj_asinh
+#define isnan   pj_isnan
+
+
+#endif /* !(defined(HAVE_C99_MATH) && HAVE_C99_MATH) */
+
+#ifdef __cplusplus
+}
+#endif
+#endif /*PROJ_MATH_H */
diff --git a/src/projects.h b/src/projects.h
index 5c8d1301..d6d6e23e 100644
--- a/src/projects.h
+++ b/src/projects.h
@@ -92,11 +92,6 @@ typedef long pj_int32;
 #define MAX_PATH_FILENAME 1024
 #endif
 
-/* prototype hypot for systems where absent */
-#if !(defined(HAVE_C99_MATH) && HAVE_C99_MATH)
-extern double hypot(double, double);
-#endif
-
 /* If we still haven't got M_PI*, we rely on our own defines.
  * For example, this is necessary when compiling with gcc and
  * the -ansi flag.