double rounding in xtime_diff()?

[Markus Kuhn]

Paul Eggert wrote on 1998-10-12 19:57 UTC:
> >     ``non-pathological'' = year range -9999 .. +9999
> 
> OK, then here's a non-pathological example of a double-rounding error:
> 
>   1970-01-01 00:00:01.851401618 UTC
> - 1970-01-01 00:00:00.000000000 UTC
> 
> On the host that I'm talking about (sparc, IEEE, round-to-even)
> your xtime_diff implementation yields
> 1.851401618000000137698179969447664916515350341796875 exactly, whereas
> 1.8514016179999999156535750444163568317890167236328125, which is also
> exactly representable as an IEEE double, is closer to the correct answer.

That's very interesting. Under Linux/Pentium_II (cpu=686, model=3,
vendor_id=GenuineIntel, stepping=4, gcc 2.7.2.3) the effect does not
show up immediately. I had to enforce using the software FP-library with
gcc -msoft-float to reproduce this effect, because without this option
with the hardware FPU I get the correct result (and this is apparently
not thanks to Intel's internal 80-bit arithmetic, because I did use gcc
-ffloat-store to enforce 64-bit rounding).

There seem to be interesting differences between IEEE and IEEE.

-----------------------------------------------------------------------
#include <stdio.h>

typedef long long int_fast64_t;
typedef long      int_fast32_t;

struct xtime {
  int_fast64_t sec;
  int_fast32_t nsec;
};

double xtime_diff(struct xtime t1, struct xtime t2)
{
  if (t1.nsec > 1000000000) t1.nsec = 1000000000;
  if (t2.nsec > 1000000000) t2.nsec = 1000000000;

  return (t1.sec - t2.sec) + (t1.nsec - t2.nsec) / 1.0e9;
}

int main()
{
  struct xtime t1 = { 1, 851401618 }, t2 = { 0, 0 };
  double d;

  d = xtime_diff(t1, t2);
  printf("%.60f\n", d);

  return 0;
}
-----------------------------------------------------------------------

Suggestions for a better portable implementation?

A somewhat related question: Given a floating-point value X,
is there in C9X a portable way to get the next smaller and the
next larger floating point value? Something like the Double'Succ(X)
and Double'Pred(X) built-in function attributes available in Ada
<http://wuarchive.wustl.edu/languages/ada/userdocs/docadalt/rm95/0305.htm>.
This might be very useful to write portable C code that is guaranteed to
deliver the closest value.

Markus

-- 
Markus G. Kuhn, Security Group, Computer Lab, Cambridge University, UK
email: mkuhn at acm.org,  home page: <http://www.cl.cam.ac.uk/~mgk25/>