Difftime code

Clive D. W. Feather clive at on-the-train.demon.co.uk
Thu Aug 5 06:03:56 UTC 2004


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     | Internet Expert | Work: <clive at demon.net>
Tel: +44 20 8495 6138  | Demon Internet  | Home: <clive at davros.org>
Fax: +44 870 051 9937  | Thus plc        | Web:  <http://www.davros.org>
Please reply to the Reply-To address, which is:  <clive at demon.net>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 470 bytes
Desc: not available
Url : http://mm.icann.org/pipermail/tz/attachments/20040805/b20f32d1/signature-0001.asc 


More information about the tz mailing list