Patch 1.49.3 for geodesic package.

Set flags for Intel compiler to prevent incorrect optimization of
arithmetic expressions #826.  Guard against nans in sincosdx #834.
Issue #831 is not addressed here (need more information...).

5 files changed, 35 insertions, 11 deletions

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 1d1baef0..8a60b540 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -33,6 +33,18 @@ elseif("${CMAKE_C_COMPILER_ID}" MATCHES "Clang")
     -Wfloat-conversion ${CMAKE_C_FLAGS}")
 endif()
 
+# Tell Intel compiler to do arithmetic accurately.  This is needed to
+# stop the compiler from ignoring parentheses in expressions like
+# (a + b) + c and from simplifying 0.0 + x to x (which is wrong if
+# x = -0.0).
+if (CMAKE_C_COMPILER_ID STREQUAL "Intel")
+  if (MSVC)
+    set (CMAKE_C_FLAGS "${CMAKE_C_FLAGS} /fp:precise")
+  else ()
+    set (CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -fprotect-parens -fsigned-zeros")
+  endif ()
+endif ()
+
 #################################################################################
 # PROJ4 CMake modules
 #################################################################################
diff --git a/configure.ac b/configure.ac
index 898d8b3b..9e4fa3d3 100644
--- a/configure.ac
+++ b/configure.ac
@@ -29,6 +29,12 @@ CFLAGS="$save_CFLAGS"
 dnl We check for headers
 AC_HEADER_STDC
 
+dnl Check flag for accurate arithmetic with Intel compiler.  This is
+dnl needed to stop the compiler from ignoring parentheses in expressions
+dnl like (a + b) + c and from simplifying 0.0 + x to x (which is wrong if
+dnl x = -0.0).
+AX_CHECK_COMPILE_FLAG([-Werror -fprotect-parens -fsigned-zeros],
+        [CFLAGS="$CFLAGS -fprotect-parens -fsigned-zeros"])
 dnl Check for C99 math functions
 save_CFLAGS="$CFLAGS"
 CFLAGS="$CFLAGS -Wall -Werror"
diff --git a/src/geodesic.c b/src/geodesic.c
index 40dffbc1..fc0b5aea 100644
--- a/src/geodesic.c
+++ b/src/geodesic.c
@@ -239,7 +239,8 @@ static void sincosdx(real x, real* sinx, real* cosx) {
   r = remquo(x, (real)(90), &q);
 #else
   r = fmod(x, (real)(360));
-  q = (int)(floor(r / 90 + (real)(0.5)));
+  /* check for NaN */
+  q = r == r ? (int)(floor(r / 90 + (real)(0.5))) : 0;
   r -= 90 * q;
 #endif
   /* now abs(r) <= 45 */
@@ -658,21 +659,21 @@ real geod_genposition(const struct geod_geodesicline* l,
     S12 = l->c2 * atan2(salp12, calp12) + l->A4 * (B42 - l->B41);
   }
 
-  if ((outmask & GEOD_LATITUDE) && plat2)
+  if ((outmask & GEOD_LATITUDE) && plat2) /* plat2 check redundant */
     *plat2 = lat2;
-  if ((outmask & GEOD_LONGITUDE) && plon2)
+  if ((outmask & GEOD_LONGITUDE) && plon2) /* plon2 check redundant */
     *plon2 = lon2;
-  if ((outmask & GEOD_AZIMUTH) && pazi2)
+  if ((outmask & GEOD_AZIMUTH) && pazi2) /* pazi2 check redundant */
     *pazi2 = azi2;
-  if ((outmask & GEOD_DISTANCE) && ps12)
+  if ((outmask & GEOD_DISTANCE) && ps12) /* ps12 check redundant */
     *ps12 = s12;
-  if ((outmask & GEOD_REDUCEDLENGTH) && pm12)
+  if ((outmask & GEOD_REDUCEDLENGTH) && pm12) /* pm12 check redundant */
     *pm12 = m12;
   if (outmask & GEOD_GEODESICSCALE) {
     if (pM12) *pM12 = M12;
     if (pM21) *pM21 = M21;
   }
-  if ((outmask & GEOD_AREA) && pS12)
+  if ((outmask & GEOD_AREA) && pS12) /* pS12 check redundant */
     *pS12 = S12;
 
   return (flags & GEOD_ARCMODE) ? s12_a12 : sig12 / degree;
@@ -1127,7 +1128,7 @@ real SinCosSeries(boolx sinp, real sinx, real cosx, const real c[], int n) {
   real ar, y0, y1;
   c += (n + sinp);              /* Point to one beyond last element */
   ar = 2 * (cosx - sinx) * (cosx + sinx); /* 2 * cos(2 * x) */
-  y0 = (n & 1) ? *--c : 0; y1 = 0;          /* accumulators for sum */
+  y0 = (n & 1) ? *--c : 0; y1 = 0;        /* accumulators for sum */
   /* Now n is even */
   n /= 2;
   while (n--) {
@@ -1894,7 +1895,9 @@ void geod_polygon_addedge(const struct geod_geodesic* g,
                           struct geod_polygon* p,
                           real azi, real s) {
   if (p->num) {              /* Do nothing is num is zero */
-    real lat = 0, lon = 0, S12 = 0;  /* Initialize S12 to stop Visual Studio warning */
+    /* Initialize S12 to stop Visual Studio warning.  Initialization of lat and
+       lon is to make CLang static analyzer happy. */
+    real lat = 0, lon = 0, S12 = 0;
     geod_gendirect(g, p->lat, p->lon, azi, GEOD_LONG_UNROLL, s,
                    &lat, &lon, 0,
                    0, 0, 0, 0, p->polyline ? 0 : &S12);
@@ -2031,6 +2034,8 @@ unsigned geod_polygon_testedge(const struct geod_geodesic* g,
   tempsum = p->A[0];
   crossings = p->crossings;
   {
+    /* Initialization of lat, lon, and S12 is to make CLang static analyzer
+       happy. */
     real lat = 0, lon = 0, s12, S12 = 0;
     geod_gendirect(g, p->lat, p->lon, azi, GEOD_LONG_UNROLL, s,
                    &lat, &lon, 0,
diff --git a/src/geodesic.h b/src/geodesic.h
index cd7c2b70..f9d612f3 100644
--- a/src/geodesic.h
+++ b/src/geodesic.h
@@ -132,7 +132,7 @@
  * The patch level of the geodesic library.  (This tracks the version of
  * GeographicLib.)
  **********************************************************************/
-#define GEODESIC_VERSION_PATCH 2
+#define GEODESIC_VERSION_PATCH 3
 
 /**
  * Pack the version components into a single integer.  Users should not rely on
diff --git a/src/geodtest.c b/src/geodtest.c
index 8cc48134..1025ce08 100644
--- a/src/geodtest.c
+++ b/src/geodtest.c
@@ -616,7 +616,8 @@ static int GeodSolve73() {
               &lat2, &lon2, &azi2);
   result += assertEquals(lat2, 81.04623, 0.5e-5);
   result += assertEquals(lon2, -170, 0.5e-5);
-  result += assertEquals(azi2, 0, 0.5e-5);
+  result += azi2 == 0 ? 0 : 1;
+  result += 1/azi2 > 0 ? 0 : 1; /* Check that azi2 = +0.0 not -0.0 */
   return result;
 }