Collect custom C99 math functions in proj_math.h

We are relying more and more on C99 math functions. On C89 systems where
those functions are not available our own custom versions of those
functions are used instead. So far these has been spread across the code
base. This commit gathers them in the same file and introduces the
proj_math.h header. The build system checks for C99 math functions. If
not found the proj_math.h header make sure that C99 functions are
defined as their pj_ equivalent.

Ideally proj_math.h is included instead of math.h. This removes the need
for any checks against HAVE_C99_MATH in the code making it easier to
read.

For this commit the functions hypot, log1p and asinh has been taken care
of.

34 files changed, 155 insertions, 93 deletions

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_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/gie.c b/src/gie.c
index 8c735b88..10b5804e 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/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_math.c b/src/pj_math.c
new file mode 100644
index 00000000..d5242636
--- /dev/null
+++ b/src/pj_math.c
@@ -0,0 +1,66 @@
+/******************************************************************************
+ * 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);
+}
+
+#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..0d28cb30
--- /dev/null
+++ b/src/proj_math.h
@@ -0,0 +1,53 @@
+/******************************************************************************
+ * 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);
+
+#define hypot   pj_hypot
+#define log1p   pj_log1p
+#define asinh   pj_asinh
+
+
+#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.