[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@demon.net> | Tel: +44 20 8495 6138 Internet Expert | Home: <clive@davros.org> | Fax: +44 870 051 9937 Demon Internet | WWW: http://www.davros.org | Mobile: +44 7973 377646 Thus plc | |