Difftime code
Clive D.W. Feather
clive at demon.net
Thu Aug 5 20:07:21 UTC 2004
[Hmm, I don't think this got out properly yesterday.]
Okay, I've looked at the problem a bit. I also note that there's some
portability problems in the code in 2004b. In particular, if time_t is a
signed integer type then the calculation of delta can involve undefined
behaviour; you cannot assume twos-complement unflagged wrap-around in
signed integer types.
time_t has to be either an integer type or a real-floating type [C99
allows it to be a complex-floating type, but POSIX doesn't]. We can
easily distinguish these cases:
if ((time_t) 0.2 != (time_t) 0.4)
// It's a real-floating type
else
// It's an integer type
If it's real-floating, then:
return (long double) time1 - time0;
is as accurate as you can get: the subtraction is done in long double
(the conversion to this is lossless) and the rounding to double happens
at the end.
If it's integer, it's a bit harder; I'll assume for simplicity that
time1>time0.
if (time0 >= 0 || time1 < 0)
return time1 - time0;
In this case, the subtraction won't overflow; if the result is
representable, it'll be rounded at the end.
In the remaining case, a simple subtraction could overflow and you're in
trouble. So we need to find the "top type":
* on a C99 implementation, this is uintmax_t;
* on a C89 implementation, this is unsigned long;
* on non-conforming implementations you have a problem, but the most
likely case is unsigned long or unsigned long long.
#if __STDC__ == 1 && __STDC_VERSION__ >= 199901L
#define HAVE_UINTMAX_T
#endif
#ifdef HAVE_UINTMAX_T
#include <stdint.h>
typedef uintmax_t top_type;
#else
#ifdef HAVE_LONG_LONG
typedef unsigned long long top_type;
#else
typedef unsigned long top_type;
#endif
#endif
This code works correctly on C89 and C99 conforming implementations, but
you may have to tune it for others.
If we now cast the two values to (top_type) and subtract, we get a value
which is either:
time1 - time0 OR
time1 - time0 - BASIS
(BASIS is 2 to the power N, where N is the number of bits in top_type).
The latter case can only happen if top_type has the same number of
*value* bits as time_t and so has the same maximum value, *and* the
subtraction overflows.
time_t s4 = time1 / 4 - time0 / 4;
top_type dt = (top_type) time1 - (top_type) time0;
top_type t4 = dt / 4;
At this point, t4 is either approximately equal to s4, in which case dt
is correct, or the right answer is dt + BASIS, in which case t4 is
approximately s4 - BASIS / 4. So:
if (t4 > s4 - (top_type)-1 / 8)
return dt;
So now we can't stay in integer types to calculate the entire value.
Rather, we need to do dt + BASIS in long double.
top_type halfBASIS = (top_type)-1 / 2 + 1;
return 2.0L * halfBASIS + dt;
I'm not in a position to do a full analysis, but let me note that we're
adding two positive values, so you don't have the issues you do with
subtraction and the addition can only lose 1ulp of precision.
Furthermore, on a binary machine the first operand can be represented
exactly.
So, putting that all together and making some cosmetic changes and minor
optimisations:
#if __STDC__ == 1 // and the implementation isn't lying about this
#if __STDC_VERSION >= 199901L
#include <stdint.h>
typedef uintmax_t top_type;
#else
typedef unsigned long top_type;
#endif
typedef long double long_double;
#else
// You'll have to fill this in using your own portability guidelines
typedef ???? top_type; // The widest range unsigned integer type
typedef ???? long_double; // The widest range floating-point type
#endif
double difftime (const time_t time1, const time_t time0)
{
time_t s4;
top_type dt;
if ((time_t) 0.2 != (time_t) 0.4)
return (long_double) time1 - time0;
if (time1 < time0)
return -difftime (time0, time1);
if (time0 >= 0 || time1 < 0)
return time1 - time0;
s4 = (time1 >> 2) - time0 / 4; // time0 < 0, so >> is a bad idea
dt = (top_type) time1 - (top_type) time0;
if ((dt >> 2) > s4 - ((top_type)-1 >> 3))
return dt;
return (long_double) 2 * ~((top_type)-1 >> 1) + dt;
}
If you don't like this approach, then I have one other suggestion.
Rather than using sizeof as a placeholder for the precision of various
types, you can use the values FLT_RADIX and LDBL_MANT_DIG (in <float.h>)
to calculate the largest integer that can be represented exactly in long
double:
static long double maxLDint = 0.0;
if (maxLDint < 100.0)
{
int i;
for (i = 0; i < LDBL_MANT_DIG; i++)
maxLDint = maxLDint * FLT_RADIX + (FLT_RADIX - 1);
}
But you still need to address some of your other bugs in your code.
--
Clive D.W. Feather | Work: <clive at demon.net> | Tel: +44 20 8495 6138
Internet Expert | Home: <clive at davros.org> | Fax: +44 870 051 9937
Demon Internet | WWW: http://www.davros.org | Mobile: +44 7973 377646
Thus plc | |
More information about the tz
mailing list