geodesic.c:sincosdx runtime error: nan is outside the range of representable values of type 'int'

[schwehr]

Original from the fuzzer:

geodesic.c:241:7: runtime error: nan is outside the range of representable values of type 'int'
    #0 sincosdx geodesic.c:241:7
    #1 geod_geninverse_int geodesic.c:773:5
    #2 geod_geninverse geodesic.c:1090:11
    #3 geod_inverse geodesic.c:1117:3
    #4 e_forward PJ_aeqd.c:118:9
    #5 pj_fwd pj_fwd.c:189:18
    #6 pj_transform pj_transform.c:405:33
    #7 TestOneProtoInput(third_party::proj4::tests::Proj4 const&) proj4_fuzzer.cc:59:3
    #8 LLVMFuzzerTestOneInput proj4_fuzzer.cc:48:1

Fuzzer input was:

cat crash-88b2552390e86ab5da34f85a907a3098204ccdf2
src: "+proj=aea"
dst: "+proj=aeqd +ell2=2"
x: nan
y: 0
z: 0

A patch that is unlikely to be good, but prevents the sanitizer failure:

q = isnan(r) ? 0.0 : (int)(floor(r / 90 + (real)(0.5)));

Assuming that I translated this correctly to cs2cs:

echo "nan 0 0" | PROJ_DEBUG=2 cs2cs_nonprod +proj=aea +to +proj=aeqd +ell2=2

pj_open_lib(proj_def.dat): call fopen(proj_def.dat) - succeeded

WARNING: Logging before InitGoogle() is written to STDERR
I0306 14:56:50.668844  199004 google3_proj4_hooks.cc:114] pj_open_lib(proj_def.dat): call fopen(proj_def.dat) - succeeded
pj_open_lib(proj_def.dat): call fopen(proj_def.dat) - succeeded

I0306 14:56:50.669717  199004 google3_proj4_hooks.cc:114] pj_open_lib(proj_def.dat): call fopen(proj_def.dat) - succeeded
pj_ellipsoid - final: a=6378137.000 f=1/298.257, errno=0
I0306 14:56:50.670044  199004 google3_proj4_hooks.cc:114] pj_ellipsoid - final: a=6378137.000 f=1/298.257, errno=0
pj_ellipsoid - final:    ellps=WGS84
I0306 14:56:50.670252  199004 google3_proj4_hooks.cc:114] pj_ellipsoid - final:    ellps=WGS84
pj_open_lib(proj_def.dat): call fopen(proj_def.dat) - succeeded

I0306 14:56:50.670587  199004 google3_proj4_hooks.cc:114] pj_open_lib(proj_def.dat): call fopen(proj_def.dat) - succeeded
pj_ellipsoid - final: a=6378137.000 f=1/298.257, errno=0
I0306 14:56:50.670862  199004 google3_proj4_hooks.cc:114] pj_ellipsoid - final: a=6378137.000 f=1/298.257, errno=0
pj_ellipsoid - final:    ellps=WGS84
I0306 14:56:50.671077  199004 google3_proj4_hooks.cc:114] pj_ellipsoid - final:    ellps=WGS84
nan	nan 0.00

[cffk]

I recommend

  q = r == r ? (int)(floor(r / 90 + (real)(0.5))) : 0;

so that it doesn't require the C99 function isnan. This will be OK, since sincosdx will return nan's in this case (as it should).

[schwehr]

I'm not sure if converting a nan to zero makes sense, but cffk's patch is better than mine.

[cffk]

The nan to zero conversion is safe in this case! (The integer q determines which quadrant the angle falls in. If the angle is a nan, we don't care which quadrant is assigned, since nan's are returned in this case.)

[cffk]

By the way, the current code is "right". The undefined value assigned in

int q = nan;

doesn't affect the result. Ideally, there would be a way to tell fuzzer to ignore this problem. Is there?

[schwehr]

It's not a fun situation. Just because it works now in the compilers that we all are generally seeing doesn't mean that undefined behavior will work for all compilers, especially in the future. I know of one compiler that can trigger SIGILL in this case. In a somewhat worse case, I have found that sometimes overrange narrowing casts don't always produce the same results across compilers and/or compiler options.
For stuff that really is okay or you don't want being fuzzed, there is FUZZING_BUILD_MODE_UNSAFE_FOR_PRODUCTION. But please don't use that in this case.

e.g. https://trac.osgeo.org/gdal/browser/trunk/gdal/gcore/gdaldataset.cpp?rev=41181#L2815

[cffk]

Hmm... I was reading "the behavior is undefined if the value cannot be represented ..." as "the value is undefined if the value cannot be represented ..." . If it's really the case the the program could abort, I'll need to handle this.

[schwehr]

I had to apply the above patch locally as I am indeed seeing SIGILL with the gie tests when a nan is passed.

[cffk]

OK, thanks. I'm looking at similar issues elsewhere in GeographicLib. I want a consistent way of handling this (and possible problems when the value is bigger than the maximum integer). For example in this case, I could use

q = (int)(maxx(-4.0, floor(r / 90 + (real)(0.5))));

[schwehr]

Perhaps a function that does a narrowing cast and clamp? modulo style guide and better naming

int FloatToIntClamp(float val, int min_val, int max_val, int nan_val) {
  if (val != val) return nan_val;
  if (val >= (float)max_val) return max_val;
  if (val <= (float)min_val) return min_val;
  return (int)val;
}

int DoubleToIntClamp(double val, int min_val, int max_val, int nan_val) {
  if (val != val) return nan_val;
  if (val >= (double)max_val) return max_val;
  if (val <= (double)min_val) return min_val;
  return (int)val;
}

// and without range valuesbut use the correct #define limits...
// I miss the C++ std::numeric_limits<int>::max() that is templated to make sure
// you get the value for the correct type.
// These should be what??? http://www.cplusplus.com/reference/climits/

int FloatToInt(float val, int nan_val) {
  if (val != val) return nan_val;
  if (val >= (float)INT_MAX) return max_val;
  if (val <= (float)INT_MIN) return min_val;
  return (int)val;
}

int DoubleToInt(double val, int nan_val) {
  if (val != val) return nan_val;
  if (val >= (double)INT_MAX) return max_val;
  if (val <= (double)INT_MIN) return min_val;
  return (int)val;
}

[kbevers]

@cffk do you have a proper understanding of this issue to supply a patch? If it is possible it would be nice to include a fix in 5.0.1. The same goes for #826 and #831.

[cffk]

Yes, I can do fixes for all three of these. Should I start with the master branch or a 5.0.1-specific branch?

[kbevers]

Just work on the master branch. I'll cherry-pick the relevant commits to the 5.0 branch once committed to master. The 5.0.1 release will then be based on the 5.0 branch that only contains bug fixes cherry-picked from master.

[schwehr]

Other ones I've run into with ASAN and fuzzing. My quick hacks to avoid triggering the issues are not necessarily okay.

--- a/src/PJ_robin.c     2018-03-14 08:30:07.000000000 -0700
+++ b/src/PJ_robin.c   2018-03-14 11:59:47.000000000 -0700
@@ -80,7 +80,7 @@
     double dphi;
     (void) P;
 
-    i = (int)floor((dphi = fabs(lp.phi)) * C1);
+    i = lp.phi == lp.phi ? (int)floor((dphi = fabs(lp.phi)) * C1) : -1;
     if( i < 0 ){
         proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION);
         return xy;
@@ -115,7 +115,7 @@
         }
     } else { /* general problem */
         /* in Y space, reduce to table interval */
-        i = (int)floor(lp.phi * NODES);
+        i = lp.phi == lp.phi ? (int)floor(lp.phi * NODES) : -1;
         if( i < 0 || i >= NODES ) {
             proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION);
             return lp;
--- a/src/nad_intr.c     2018-03-14 08:30:07.000000000 -0700
+++ b/src/nad_intr.c   2018-03-14 11:59:47.000000000 -0700
@@ -10,8 +10,14 @@
        long index;
        int in;
 
-       indx.lam = (int)floor(t.lam /= ct->del.lam);
-       indx.phi = (int)floor(t.phi /= ct->del.phi);
+       indx.lam =
+      t.lam == t.lam && ct->del.lam == ct->del.lam
+      ? (int)floor(t.lam /= ct->del.lam)
+      : 0;
+       indx.phi =
+      t.phi == t.phi && ct->del.phi == ct->del.phi
+      ? (int)floor(t.phi /= ct->del.phi)
+      : 0;
        frct.lam = t.lam - indx.lam;
        frct.phi = t.phi - indx.phi;
        val.lam = val.phi = HUGE_VAL;
--- a/src/pj_gridinfo.c  2018-03-14 08:23:31.000000000 -0700
+++ b/src/pj_gridinfo.c        2018-03-14 11:59:48.000000000 -0700
@@ -916,10 +916,12 @@
         gilist->format = "ctable2";
         gilist->ct = ct;
 
+        if (ct)
         pj_log( ctx, PJ_LOG_DEBUG_MAJOR,
                 "Ctable2 %s %dx%d: LL=(%.9g,%.9g) UR=(%.9g,%.9g)\n",
                 ct->id,
--- a/src/proj_etmerc.c  2018-03-14 08:30:07.000000000 -0700
+++ b/src/proj_etmerc.c        2018-03-14 11:59:48.000000000 -0700
@@ -358,7 +358,10 @@
     }
     else /* nearest central meridian input */
     {
-        zone = (int)(floor ((adjlon (P->lam0) + M_PI) * 30. / M_PI));
+        zone =
+            P->lam0 == P->lam0
+            ? (int)(floor ((adjlon (P->lam0) + M_PI) * 30. / M_PI))
+            : 0;
         if (zone < 0)
             zone = 0;
         else if (zone >= 60)

[cffk]

@schwehr Btw, could you run your tests with HAVE_C99_MATH (e.g., by removing the -std=c89 flag in CMakeLists.txt? I'm curious if remquo (in sincosdx) on a NaN is OK. (The documentation says that a NaN is returned but doesn't say what the integer argument should become.)

[schwehr]

Thanks for pointing out HAVE_C99_MATH. While I have looked at the code many of times, it never registered on me that I should be setting it. Sigh. I am now.

Tried it and it seems to work. Checking the behavior with gunit, it looks safe.

However quo is not always set. And that behavior is not called out in the man page on debian testing linux. I didn't look into math error.

• http://man7.org/linux/man-pages/man3/remquo.3.html
• http://man7.org/linux/man-pages/man7/math_error.7.html

RETURN VALUE
On success, these functions return the same value as the analogous functions described in remainder(3).

   If x or y is a NaN, a NaN is returned.

   If x is an infinity, and y is not a NaN, a domain error occurs, and a NaN is returned.

   If y is zero, and x is not a NaN, a domain error occurs, and a NaN is returned.

ERRORS
See math_error(7) for information on how to determine whether an error has occurred when calling these functions.

   The following errors can occur:

   Domain error: x is an infinity or y is 0, and the other argument is not a NaN
          An invalid floating-point exception (FE_INVALID) is raised.

   These functions do not set errno.

quick gunit test that I just wrote for remquo

// Copyright 2018 Google Inc. All Rights Reserved.
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
//     http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.

#include <cmath>
#include <limits>

#include "testing/base/public/gmock.h"
#include "testing/base/public/gunit.h"

// remquo, remquof, remquol - remainder and part of quotient

namespace {

TEST(RemquoTest, DoesFoo) {
  constexpr double nan = std::numeric_limits<double>::quiet_NaN();
  constexpr double inf = std::numeric_limits<double>::infinity();

  int quo = -999;
  EXPECT_DOUBLE_EQ(0.0, remquo(0.0, 1.0, &quo));
  EXPECT_EQ(0, quo);

  quo = -999;
  EXPECT_TRUE(std::isnan(remquo(0.0, 0.0, &quo)));
  EXPECT_EQ(-999, quo);

  quo = -999;
  EXPECT_TRUE(std::isnan(remquo(nan, nan, &quo)));
  EXPECT_EQ(-999, quo);

  quo = -999;
  EXPECT_TRUE(std::isnan(remquo(0.0, nan, &quo)));
  EXPECT_EQ(-999, quo);

  quo = -999;
  EXPECT_TRUE(std::isnan(remquo(nan, 0.0, &quo)));
  EXPECT_EQ(-999, quo);

  quo = -999;
  EXPECT_TRUE(std::isnan(remquo(inf, inf, &quo)));
  EXPECT_EQ(-999, quo);

  quo = -999;
  EXPECT_DOUBLE_EQ(0.0, remquo(0.0, inf, &quo));
  EXPECT_EQ(0, quo);

  quo = -999;
  EXPECT_TRUE(std::isnan(remquo(inf, 0.0, &quo)));
  EXPECT_EQ(-999, quo);

  quo = -999;
  EXPECT_TRUE(std::isnan(remquo(-inf, -inf, &quo)));
  EXPECT_EQ(-999, quo);

  quo = -999;
  EXPECT_TRUE(std::isnan(remquo(nan, inf, &quo)));
  EXPECT_EQ(-999, quo);

  quo = -999;
  EXPECT_TRUE(std::isnan(remquo(-inf, nan, &quo)));
  EXPECT_EQ(-999, quo);

  quo = -999;
  EXPECT_TRUE(std::isnan(remquo(nan, -inf, &quo)));
  EXPECT_EQ(-999, quo);

  quo = -999;
  EXPECT_TRUE(std::isnan(remquo(-inf, nan, &quo)));
  EXPECT_EQ(-999, quo);
}

}  // namespace

[cffk]

Whoops, I didn't mean to close this issue. I only addressed the one issue of NaN -> int conversion in geodesic.c.

[kbevers]

Fixed in #887.