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

Collect custom C99 math functions in proj_math.h

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)

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"])

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 \

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

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,

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.";
 

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=";

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");
 

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
 

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");

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;

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;

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.";
 

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";
 

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";

src/PJ_merc.c

@@ -1,37 +1,18 @@
 #define PJ_LIB__
 #include "proj_internal.h"
 #include "proj.h"
+#include "proj_math.h"
 #include "projects.h"
 #include <float.h>
 
 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));
 }

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)";

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,

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=";
 

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.";

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;

src/PJ_sconics.c

@@ -2,6 +2,7 @@
 #include <errno.h>
 #include "proj.h"
 #include "projects.h"
+#include "proj_math.h"
 
 
 enum Type {

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";

src/PJ_sterea.c

@@ -26,6 +26,7 @@
 #define PJ_LIB__
 #include <errno.h>
 #include "projects.h"
+#include "proj_math.h"
 
 
 struct pj_opaque {

src/PJ_tpeqd.c

@@ -1,6 +1,7 @@
 #define PJ_LIB__
 #include <errno.h>
 #include "proj.h"
+#include "proj_math.h"
 #include "projects.h"
 
 

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

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"

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

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 \

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
 

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;