Re: 961002: problem with ncdump (fwd)

>To: address@hidden
>From: Patrick Guio <address@hidden>
>Subject: problem with ncdump (fwd)
>Organization: .
>Keywords: 199610020843.AA13593

Hi Patrick,

As I suspected, the bug was caused by unusual (but standard) results of
arithmetic and comparisons with an IEEE infinity value (rather than an
IEEE NaN).  I have appended a patch to the ncdump/vardata.c file in the
netCDF 2.4.3 distribution.  At your convenience, please apply it and let
me know whether it fixes the problem you are seeing.  It seems to work
OK here.  If you can't apply a patch in this form, please let me know,
and I'll send you a copy of the patched file instead.

Thanks for reporting the bug!

--Russ

Index: vardata.c
===================================================================
RCS file: /upc/share/CVS/netcdf/ncdump/vardata.c,v
retrieving revision 1.43
diff -c -r1.43 vardata.c
*** 1.43        1996/06/07 14:46:14
--- vardata.c   1996/10/02 19:27:47
***************
*** 128,190 ****
      struct ncvar *varp;               /* variable */
      void *valp;                       /* value, interpreted using varp->type 
*/
  {
!     union {
!       char *cp;
!       short *sp;
!       nclong *lp;
!       float *fp;
!       double *dp;
!     } vp, fillp;
! 
      switch (varp->type) {
      case NC_BYTE:
!       fillp.cp = &varp->fillval.charv;
!       vp.cp = (char *)valp;
!       if (varp->has_fillval && *fillp.cp == *vp.cp) {
            (void) sprintf(sout, FILL_STRING);
        } else {
!           (void) sprintf(sout, fmt, *vp.cp);
        }
        break;
        case NC_SHORT:
!       fillp.sp = &varp->fillval.shortv;
!       vp.sp = (short *)valp;
!       if (varp->has_fillval && *fillp.sp == *vp.sp) {
              (void) sprintf(sout, FILL_STRING);
        } else {
!             (void) sprintf(sout, fmt, *vp.sp);
        }
        break;
        case NC_LONG:
!       fillp.lp = &varp->fillval.longv;
!       vp.lp = (nclong *)valp;
!       if (varp->has_fillval && *fillp.lp == *vp.lp) {
              (void) sprintf(sout, FILL_STRING);
        } else {
!             (void) sprintf(sout, fmt, *vp.lp);
        }
        break;
        case NC_FLOAT:
!       fillp.fp = &varp->fillval.floatv;
!       vp.fp = (float *)valp;
  #define absval(x)  ( (x) < 0 ? -(x) : (x) )
        if(varp->has_fillval &&
!          (*vp.fp > 0) == (*fillp.fp > 0) && /* prevents potential overflow */
!          (absval(*vp.fp - *fillp.fp) <= absval(float_eps * *fillp.fp))) {
              (void) sprintf(sout, FILL_STRING);
        } else {
!             (void) sprintf(sout, fmt, *vp.fp);
        }
        break;
        case NC_DOUBLE:
!       fillp.dp = &varp->fillval.doublev;
!       vp.dp = (double *)valp;
        if(varp->has_fillval &&
!          (*vp.dp > 0) == (*fillp.dp > 0) && /* prevents potential overflow */
!          (absval(*vp.dp - *fillp.dp) <= absval(double_eps * *fillp.dp))) {
              (void) sprintf(sout, FILL_STRING);
        } else {
!             (void) sprintf(sout, fmt, *vp.dp);
        }
        break;
        default:
--- 128,201 ----
      struct ncvar *varp;               /* variable */
      void *valp;                       /* value, interpreted using varp->type 
*/
  {
!     char cval, cfill;
!     short sval, sfill;
!     long lval, lfill;
!     float fval, ffill;
!     double dval, dfill;
!     
      switch (varp->type) {
      case NC_BYTE:
!       cfill = varp->fillval.charv;
!       cval = *(char *)valp;
!       if (varp->has_fillval && cfill == cval) {
            (void) sprintf(sout, FILL_STRING);
        } else {
!           (void) sprintf(sout, fmt, cval);
        }
        break;
        case NC_SHORT:
!       sfill = varp->fillval.shortv;
!       sval = *(short *)valp;
!       if (varp->has_fillval && sfill == sval) {
              (void) sprintf(sout, FILL_STRING);
        } else {
!             (void) sprintf(sout, fmt, sval);
        }
        break;
        case NC_LONG:
!       lfill = varp->fillval.longv;
!       lval = *(nclong *)valp;
!       if (varp->has_fillval && lfill == lval) {
              (void) sprintf(sout, FILL_STRING);
        } else {
!             (void) sprintf(sout, fmt, lval);
        }
        break;
        case NC_FLOAT:
!       ffill = varp->fillval.floatv;
!       fval = *(float *)valp;
  #define absval(x)  ( (x) < 0 ? -(x) : (x) )
+       /* Assuming IEEE arithmetic, the following returns 0 for infinity */
+ #define IS_FINITE(x) (((x)-(x))==((x)-(x)))
+       
        if(varp->has_fillval &&
!          IS_FINITE(ffill) &&  /* needed, because IEEE infinity actually
!                                * passes the next two tests with any
!                                * value! */
!          (fval > 0) == (ffill > 0) && /* prevents potential overflow */
!          (absval(fval - ffill) <= absval(float_eps * ffill))) {
              (void) sprintf(sout, FILL_STRING);
+       } else if (!IS_FINITE(ffill) && !IS_FINITE(fval) && ffill==fval) {
+             (void) sprintf(sout, FILL_STRING);
        } else {
!             (void) sprintf(sout, fmt, fval);
        }
        break;
        case NC_DOUBLE:
!       dfill = varp->fillval.doublev;
!       dval = *(double *)valp;
        if(varp->has_fillval &&
!          IS_FINITE(dfill) &&  /* needed, because IEEE infinity actually
!                                * passes the next two tests with any
!                                * value! */
!          (dval > 0) == (dfill > 0) && /* prevents potential overflow */
!          (absval(dval - dfill) <= absval(double_eps * dfill))) {
!             (void) sprintf(sout, FILL_STRING);
!       } else if (!IS_FINITE(dfill) && !IS_FINITE(dval) && dfill==dval) {
              (void) sprintf(sout, FILL_STRING);
        } else {
!             (void) sprintf(sout, fmt, dval);
        }
        break;
        default: