
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@demon.net> Tel: +44 20 8495 6138 | Demon Internet | Home: <clive@davros.org> Fax: +44 870 051 9937 | Thus plc | Web: <http://www.davros.org> Please reply to the Reply-To address, which is: <clive@demon.net>

Wow, that's pretty a complicated way to subtract two numbers. :-) If we can assume uintmax_t (by typedefing it on older hosts), isn't there a much simpler approach entirely? Something like the following. The idea behind the sizeof test is to avoid "long double" if it's safe to do so, since long double is expensive on some hosts. This use of sizeof conforms even to the weird C99 rules. (Tricky, huh?) #define TYPE_FLOATING(type) ((type) 0.5 != 0) #define TYPE_SIGNED(type) (((type) -1) < 0) #define TYPE_BIT(type) (sizeof (type) * CHAR_BIT) static double simple_difftime (time_t time1, time_t time0) { if (TYPE_SIGNED (time_t)) return (uintmax_t) time1 - (uintmax_t) time0; else return time1 - time0; } double difftime (time_t time1, time_t time0) { if (TYPE_BIT (time_t) <= DBL_MANT_DIG || (TYPE_FLOATING (time_t) && sizeof (time_t) != sizeof (long double))) return (double) time1 - (double) time0; if (TYPE_FLOATING (time_t)) return (long_double) time1 - (long_double) time0; if (time1 < time0) return -simple_difftime (time0, time1); else return simple_difftime (time1, time0); }

Paul Eggert said:
Wow, that's pretty a complicated way to subtract two numbers. :-)
Only if you want the right answer :-)
If we can assume uintmax_t (by typedefing it on older hosts), isn't there a much simpler approach entirely?
No. [Note: my top_type is effectively that typedef.]
The idea behind the sizeof test is to avoid "long double" if it's safe to do so, since long double is expensive on some hosts.
Is it *that* expensive? I thought you were trying to avoid rounding errors.
#define TYPE_FLOATING(type) ((type) 0.5 != 0) #define TYPE_SIGNED(type) (((type) -1) < 0) #define TYPE_BIT(type) (sizeof (type) * CHAR_BIT) [...] double difftime (time_t time1, time_t time0) { if (TYPE_BIT (time_t) <= DBL_MANT_DIG
That's a conservative test I have no objection to. It's *very* conservative, since if FLT_RADIX isn't 2 it will seriously underestimate how big double is. You might be better off comparing time1 and time0 with DBL_EPSILON or use the maxLDint trick I described. But if you're worried about efficiency, why are you doing this in floating point when they're integers?
|| (TYPE_FLOATING (time_t) && sizeof (time_t) != sizeof (long double)))
Okay, this proves that time_t isn't long double. Clever.
return (double) time1 - (double) time0; if (TYPE_FLOATING (time_t)) return (long_double) time1 - (long_double) time0;
if (time1 < time0) return -simple_difftime (time0, time1); else return simple_difftime (time1, time0); }
Moved here for expositional purposes:
static double simple_difftime (time_t time1, time_t time0) { if (TYPE_SIGNED (time_t)) return (uintmax_t) time1 - (uintmax_t) time0; else return time1 - time0; }
That will sometimes get the answer badly wrong. The problem occurs when time_t is signed and the maximum value of time_t is the same as the maximum value of uintmax_t. For example, a C89 system where long is 60 bits including sign and unsigned long is 59 bits. On such systems, the maximum possible difference is greater than the maximum value of uintmax_t, and your subtract will get it wrong. I gate this case by looking for time1 >=0 and time0 < 0. In fact, you can be safer than that: #define HALFMAX ((uintmax_t)-1 >> 1) if (time1 <= HALFMAX && (time0 >= 0 || (uintmax_t) time0 >= -HALFMAX)) return (uintmax_t) time1 - (uintmax_t) time0; However, the remaining cases have to allow for overflow in the subtraction, and that's the complicated bit. -- 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 | |

Being in a hurry this morning, I wrote:
|| (TYPE_FLOATING (time_t) && sizeof (time_t) != sizeof (long double)))
Okay, this proves that time_t isn't long double. Clever.
That is, if this is true then time_t isn't long double. Still clever.
The problem occurs when time_t is signed and the maximum value of time_t is the same as the maximum value of uintmax_t.
On C99, this can only happen if UINTMAX_MAX == INTMAX_MAX. On C89 it can only happen if ULONG_MAX == LONG_MAX. If you're willing to ensure these symbols are present on non-Standard systems, you can filter out many cases at compile time.
For example, a C89 system where long is 60 bits including sign and unsigned long is 59 bits. On such systems, the maximum possible difference is greater than the maximum value of uintmax_t, and your subtract will get it wrong.
The worst case on such a system (using ** for powers) is: time1 == 2**59 - 1 time0 == -2**59 UINTMAX_MAX == 2**59 - 1 time1-time0 == 2**60 - 1, which will overflow and come out as 2**59 - 1
I gate this case by looking for time1 >=0 and time0 < 0. In fact, you can be safer than that:
#define HALFMAX ((uintmax_t)-1 >> 1)
if (time1 <= HALFMAX && (time0 >= 0 || (uintmax_t) time0 >= -HALFMAX)) return (uintmax_t) time1 - (uintmax_t) time0;
There are in fact many ways to filter out easy cases: * If time1 and time0 have the same sign (treating 0 as positive) then you can't have an overflow (original posting). * If time1 and -time0 both have the top bit zero, you can't have an overflow (posting quoted just above). * If time1/2 - time0/2 is less than UINTMAX_MAX/2, you can't have an overflow. Take your pick. Attempting to summarize yet again: * If time_t is a floating point type, we can just do the calculation in long double. Furthermore, unless it's long double we can do it in double. [Note, however, that C99 left open the possibility of more floating point types in future standards.] * If time_t is an integer type, we can nearly always convert to uintmax_t, subtract, and convert to double. * There are a few implementations where this subtraction could overflow. There are compile- and run-time tests which can identify whether a particular call to difftime is "at risk"; if it is, it then becomes necessary to either: - determine whether the simple subtraction is correct or is too low by UINTMAX_MAX + 1, and adjust the result in the latter case; - do the calculation in floating point, risking a loss of precision. Any code which fails to do all of this is sometimes going to get the wrong answer. Any code which does all of this is going to look "pretty complicated". -- 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 | |

"Clive D.W. Feather" <clive@demon.net> writes:
The idea behind the sizeof test is to avoid "long double" if it's safe to do so, since long double is expensive on some hosts.
Is it *that* expensive? I thought you were trying to avoid rounding errors.
Avoiding rounding errors is the primary goal, but efficiency is also important.
You might be better off comparing time1 and time0 with DBL_EPSILON or use the maxLDint trick I described.
I'd rather avoid run-time tests. That is, I would rather have the normal case execute code without any conditional branches.
But if you're worried about efficiency, why are you doing this in floating point when they're integers?
On the modern processors where this code typically executes, floating point is pretty fast. (Except maybe long double, I suppose.) Conditional branches hurt performance, though.
The problem occurs when time_t is signed and the maximum value of time_t is the same as the maximum value of uintmax_t.
I see. But this will happen only on weird hosts with padding bits, right? On other hosts (i.e. the vast majority) we don't need to check for overflow.
s4 = (time1 >> 2) - time0 / 4; // time0 < 0, so >> is a bad idea
We can't use >> here, since time_t might be a floating type. I looked at your code more carefully and have the following thoughts. * The idea of adding 2 * ~((top_type)-1 >> 1) is basically the same idea as the current difftime's approach of subtracting 2 * (long_double) hibit, except it's being done with unsigned arithmetic and therefore avoids integer-overflow problems. * I suspect the add-the-high-bit approach will suffer from double-rounding on some theoretical platforms (given that it clearly suffers from similar problems in the existing difftime.c), but I know of no actual platforms where double-rounding will occur so it's hard to test whether it will introduce an error. * I didn't understand why you computed time1 / 4 - time0 / 4. Presumably this was to avoid overflow, but surely it suffices to divide by 2, and compute time1 / 2 - time0 / 2. * I got confused by the "t4 > s4 - (top_type)-1 / 8" business; I don't really know what that code is doing. Sorry. * Even on weird hosts I'd rather avoid runtime tests like time1 >=0 and time0 < 0. * As before, the code will have more rounding errors on hsots where FLT_RADIX is not a power of 2. But I don't think we need to worry about this, from a practical viewpoint. Here's the result of my cogitating on this. From a practical point of view, the main improvement over the current difftime is that it avoids double-rounding on platforms with 64-bit time_t and 64-bit long double; this is a somewhat-orthogonal change but it is a real win. But I hope it also addresses your concerns about correctness on hosts that have padding bits. /* This part will be configured differently for pre-C99 hosts. */ #include <stdint.h> #include <float.h> #define uintmax uintmax_t #define long_double long double #include <time.h> #include <limits.h> #define TYPE_FLOATING(type) ((type) 0.4 != 0) #define TYPE_SIGNED(type) (((type) -1) < 0) #define TYPE_BIT(type) (sizeof (type) * CHAR_BIT) /* ** Return TIME1 - TIME0, where TIME0 <= TIME1. ** This function is executed only if time_t is an integer type. */ static double simple_difftime (time_t time1, time_t time0) { /* ** Optimize the common special cases where time_t is unsigned, ** or can be converted to uintmax without losing information. */ if (!TYPE_SIGNED (time_t)) return time1 - time0; else if (INTMAX_MAX <= UINTMAX_MAX / 2) return (uintmax) time1 - (uintmax) time0; else { /* ** Weird. uintmax has padding bits, and time_t is signed. ** Check for overflow: compare delta/2 to (time1/2 - time0/2). ** Overflow occurred if they differ by more than a small slop. */ uintmax delta = (uintmax) time1 - (uintmax) time0; uintmax hdelta = delta / 2; time_t ht1 = time1 / 2; time_t ht0 = time0 / 2; time_t dht = ht1 - ht0; /* ** "h" here means half. By simple range analysis, we have: ** -0.5 <= ht1 - time1/2 <= 0.5 ** -0.5 <= ht0 - time0/2 <= 0.5 ** -1.0 <= dht - (time1 - time0)/2 <= 1.0 ** If overflow has not occurred, we also have: ** -0.5 <= hdelta - (time1 - time0)/2 <= 0 ** -1.0 <= dht - hdelta <= 1.5 ** and since dht - hdelta is an integer, we also have: ** -1 <= dht - hdelta <= 1 ** or equivalently: ** 0 <= dht - hdelta + 1 <= 2 ** In the above analysis, all the operators have ** their exact mathematical semantics, not C semantics. ** However, dht - hdelta + 1 is unsigned in C, ** so it need not be compared to zero. */ if (dht - hdelta + 1 <= 2) return delta; else { /* ** Repair delta overflow. ** ** The following expression contains a second rounding, ** so the result may not be the closest to the true answer. ** This problem occurs only with very large differences, ** It's too painful to fix this portably. ** We are not alone in this problem; ** some C compilers round twice when converting ** large unsigned types to small floating types, ** so if time_t is unsigned the "return time1 - time0" above ** has the same double-rounding problem with those compilers. */ long_double hibit = ~(UINTMAX_MAX / 2); return delta + 2 * hibit; } } } double difftime (time_t time1, time_t time0) { if (TYPE_BIT (time_t) <= DBL_MANT_DIG || (TYPE_FLOATING (time_t) && sizeof (time_t) != sizeof (long_double))) return (double) time1 - (double) time0; if (TYPE_BIT (time_t) <= LDBL_MANT_DIG || TYPE_FLOATING (time_t)) return (long_double) time1 - (long_double) time0; if (time1 < time0) return -simple_difftime (time0, time1); else return simple_difftime (time1, time0); }

Paul Eggert said:
The idea behind the sizeof test is to avoid "long double" if it's safe to do so, since long double is expensive on some hosts. Is it *that* expensive? I thought you were trying to avoid rounding errors. Avoiding rounding errors is the primary goal, but efficiency is also important.
Well, your choice.
You might be better off comparing time1 and time0 with DBL_EPSILON or use the maxLDint trick I described. I'd rather avoid run-time tests. That is, I would rather have the normal case execute code without any conditional branches.
Hmm.
But if you're worried about efficiency, why are you doing this in floating point when they're integers? On the modern processors where this code typically executes, floating point is pretty fast. (Except maybe long double, I suppose.)
I'm surprised that long double takes more than twice a normal double or Much more than the function call, particularly if:
Conditional branches hurt performance, though.
I'm obviously out of date.
The problem occurs when time_t is signed and the maximum value of time_t is the same as the maximum value of uintmax_t. I see. But this will happen only on weird hosts with padding bits, right?
Not that weird.
On other hosts (i.e. the vast majority) we don't need to check for overflow.
In particular, it only happens if: #if UINTMAX_MAX == INTMAX_MAX (so if you have those symbols, or ULONG_MAX == LONG_MAX on C89, you can test at compile time.
s4 = (time1 >> 2) - time0 / 4; // time0 < 0, so >> is a bad idea We can't use >> here, since time_t might be a floating type.
Oops. The code will only get executed if time_t is integer, but you're right. Silly me. A test for time1 >= 0 would be useful here; firstly, if it's false, overflow can't happen. Secondly, it it's true, the compiler has enough information to know that /4 can be converted to >>2.
I looked at your code more carefully and have the following thoughts. * The idea of adding 2 * ~((top_type)-1 >> 1) is basically the same idea as the current difftime's approach of subtracting 2 * (long_double) hibit, except it's being done with unsigned arithmetic and therefore avoids integer-overflow problems.
Right. Remember that overflow in signed integers is allowed to produce *any* answer or even raise a signal.
* I suspect the add-the-high-bit approach will suffer from double-rounding on some theoretical platforms (given that it clearly suffers from similar problems in the existing difftime.c), but I know of no actual platforms where double-rounding will occur so it's hard to test whether it will introduce an error.
Not my area of expertise.
* I didn't understand why you computed time1 / 4 - time0 / 4. Presumably this was to avoid overflow, but surely it suffices to divide by 2, and compute time1 / 2 - time0 / 2.
At the time I wrote it I'd convinced myself there was a corner case that still went wrong. Dividing by any number greater than 2 avoided it.
* I got confused by the "t4 > s4 - (top_type)-1 / 8" business; I don't really know what that code is doing. Sorry.
s4 is 1/4 of a value within +/- 4 of the true answer. t4 is 1/4 of dt, which is what you get if you subtract in top_type. dt is either correct or it's wrong by UINTMAX_MAX+1. So t4 - s4 is either close to zero or is about -UINTMAX_MAX/4. That test determines which.
* Even on weird hosts I'd rather avoid runtime tests like time1 >=0 and time0 < 0.
Yet you're happy with an extra function call?
Here's the result of my cogitating on this. From a practical point of view, the main improvement over the current difftime is that it avoids double-rounding on platforms with 64-bit time_t and 64-bit long double; this is a somewhat-orthogonal change but it is a real win. But I hope it also addresses your concerns about correctness on hosts that have padding bits.
/* This part will be configured differently for pre-C99 hosts. */ #include <stdint.h> #include <float.h> #define uintmax uintmax_t #define long_double long double
#include <time.h> #include <limits.h>
#define TYPE_FLOATING(type) ((type) 0.4 != 0) #define TYPE_SIGNED(type) (((type) -1) < 0) #define TYPE_BIT(type) (sizeof (type) * CHAR_BIT)
/* ** Return TIME1 - TIME0, where TIME0 <= TIME1. ** This function is executed only if time_t is an integer type. */ static double simple_difftime (time_t time1, time_t time0) { /* ** Optimize the common special cases where time_t is unsigned, ** or can be converted to uintmax without losing information. */ if (!TYPE_SIGNED (time_t)) return time1 - time0; else if (INTMAX_MAX <= UINTMAX_MAX / 2) return (uintmax) time1 - (uintmax) time0;
This second test can be "INTMAX_MAX != UINTMAX_MAX", and can be a #if.
else { /* ** Weird. uintmax has padding bits, and time_t is signed. ** Check for overflow: compare delta/2 to (time1/2 - time0/2). ** Overflow occurred if they differ by more than a small slop. */ uintmax delta = (uintmax) time1 - (uintmax) time0; uintmax hdelta = delta / 2; time_t ht1 = time1 / 2; time_t ht0 = time0 / 2; time_t dht = ht1 - ht0; /* ** "h" here means half. By simple range analysis, we have: ** -0.5 <= ht1 - time1/2 <= 0.5 ** -0.5 <= ht0 - time0/2 <= 0.5 ** -1.0 <= dht - (time1 - time0)/2 <= 1.0 ** If overflow has not occurred, we also have: ** -0.5 <= hdelta - (time1 - time0)/2 <= 0 ** -1.0 <= dht - hdelta <= 1.5 ** and since dht - hdelta is an integer, we also have: ** -1 <= dht - hdelta <= 1 ** or equivalently: ** 0 <= dht - hdelta + 1 <= 2 ** In the above analysis, all the operators have ** their exact mathematical semantics, not C semantics. ** However, dht - hdelta + 1 is unsigned in C, ** so it need not be compared to zero. */ if (dht - hdelta + 1 <= 2) return delta;
You've basically done the same as the code you said you didn't understand! * I used /4 instead of /2 because I wanted to ensure we never went near the maximum values of types. * If there *is* an overflow, dht - hdelta is going to be around UINTMAX_MAX/2. So I used the eqivalent of UINTMAX_MAX/4 as the gateway; I'd certainly use a value bigger than 2 just in case I'd made an error in the analysis. In fact, it might be a good idea to insert, at this point: assert (dht - hdelta + 1 >= UINTMAX_MAX/2 - 2 && dht - hdelta + 1 <= UINTMAX_MAX/2 + 2); to ensure nothing's gone wrong (I *think* +/-2 suffices).
else { /* ** Repair delta overflow. ** ** The following expression contains a second rounding, ** so the result may not be the closest to the true answer. ** This problem occurs only with very large differences, ** It's too painful to fix this portably. ** We are not alone in this problem; ** some C compilers round twice when converting ** large unsigned types to small floating types, ** so if time_t is unsigned the "return time1 - time0" above ** has the same double-rounding problem with those compilers. */ long_double hibit = ~(UINTMAX_MAX / 2); return delta + 2 * hibit; } } }
double difftime (time_t time1, time_t time0) { if (TYPE_BIT (time_t) <= DBL_MANT_DIG || (TYPE_FLOATING (time_t) && sizeof (time_t) != sizeof (long_double))) return (double) time1 - (double) time0; if (TYPE_BIT (time_t) <= LDBL_MANT_DIG || TYPE_FLOATING (time_t)) return (long_double) time1 - (long_double) time0;
if (time1 < time0) return -simple_difftime (time0, time1); else return simple_difftime (time1, time0); }
-- 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 | |

Thanks for the tour of some of the darker corners of C; I didn't realize all the ins and outs of hosts with padding bits. I responded to some of your points (e.g., avoiding function calls), and compiled and looked at the assembly code output on both sparc and x86 and tuned it a bit more, and came up with the following proposed rewrite of difftime.c. It does fix some problems on fairly-popular hosts (e.g., it avoids double-rounding on 64-bit sparc and probably several other 64-bit hosts) and it fixes the padding-bit bugs you noted (at least in theory: I can't easily test this), so it looks like a win overall. /* ** This file is in the public domain, so clarified as of ** June 5, 1996 by Arthur David Olson (arthur_david_olson@nih.gov). */ #ifndef lint #ifndef NOID static char elsieid[] = "@(#)difftime.c 7.9"; #endif /* !defined NOID */ #endif /* !defined lint */ /*LINTLIBRARY*/ #include "private.h" /* ** Algorithm courtesy Paul Eggert (eggert@cs.ucla.edu). ** ** Most other code assumes that time_t is an integer type without ** padding bits, and that integer arithmetic is modular two's ** complement without overflow traps, but (just for fun) this works ** even if time_t is an integer type with padding bits, or a real ** floating type, and it works even if signed integer overflow ** has undefined behavior. */ #include <float.h> #define TYPE_FLOATING(type) ((type) 0.4 != 0) #if !defined HAVE_LONG_DOUBLE && defined __STDC__ #define HAVE_LONG_DOUBLE 1 #endif /* !defined HAVE_LONG_DOUBLE && defined __STDC__ */ #ifndef HAVE_LONG_DOUBLE #define HAVE_LONG_DOUBLE 0 #endif /* !defined HAVE_LONG_DOUBLE */ #if HAVE_LONG_DOUBLE #define long_double long double #endif /* HAVE_LONG_DOUBLE */ #if !HAVE_LONG_DOUBLE #define long_double double #endif /* !HAVE_LONG_DOUBLE */ #ifndef HAVE_STDINT_H #define HAVE_STDINT_H (199901 <= __STDC_VERSION__) #endif /* !defined HAVE_STDINT_H */ #if HAVE_STDINT_H #include <stdint.h> #define uintmax uintmax_t #endif /* HAVE_STDINT_H */ #if !defined uintmax && defined ULLONG_MAX #define uintmax unsigned long long int #endif /* !defined uintmax && defined ULLONG_MAX */ #ifndef uintmax #define uintmax unsigned long int #endif /* defined uintmax */ #ifndef UINTMAX_MAX #define UINTMAX_MAX ((uintmax) -1) #endif /* !defined UINTMAX_MAX */ #if !defined INTMAX_MAX && defined LLONG_MAX #define INTMAX_MAX LLONG_MAX #endif /* !defined INTMAX_MAX && defined LLONG_MAX */ #ifndef INTMAX_MAX #define INTMAX_MAX LONG_MAX #endif /* !defined INTMAX_MAX */ double difftime(time1, time0) time_t time1; time_t time0; { int time1_is_smaller; double delta; /* ** Use floating point if there should be no double-rounding error. ** However, avoid long double if it must be wider than needed, ** as it's sometimes much more expensive in these cases ** (e.g., 64-bit sparc). */ if (TYPE_BIT(time_t) <= DBL_MANT_DIG || (TYPE_FLOATING(time_t) && sizeof(time_t) < sizeof(long_double))) { double t1 = time1; double t0 = time0; return t1 - t0; } if ((TYPE_BIT(time_t) <= LDBL_MANT_DIG && (TYPE_BIT(time_t) == LDBL_MANT_DIG || (TYPE_SIGNED(time_t) && UINTMAX_MAX / 2 < INTMAX_MAX))) || TYPE_FLOATING(time_t)) { long_double t1 = time1; long_double t0 = time0; return t1 - t0; } time1_is_smaller = time1 < time0; if (time1_is_smaller) { time_t t = time1; time0 = time1; time1 = t; } /* ** Now time0 <= time1, and time_t is an integer type. ** Optimize the common special cases where time_t is unsigned, ** or can be converted to uintmax without losing information. */ if (! TYPE_SIGNED(time_t)) delta = time1 - time0; else { uintmax t1 = time1; uintmax t0 = time0; uintmax dt = t1 - t0; delta = dt; if (UINTMAX_MAX / 2 < INTMAX_MAX) { /* ** uintmax has padding bits, and time_t is signed. ** Check for overflow: compare dt/2 to (time1/2 - ** time0/2). Overflow occurred if they differ by ** more than a small slop. */ uintmax hdt = dt / 2; time_t ht1 = time1 / 2; time_t ht0 = time0 / 2; time_t dht = ht1 - ht0; /* ** "h" here means half. By range analysis, we have: ** -0.5 <= ht1 - time1/2 <= 0.5 ** -0.5 <= ht0 - time0/2 <= 0.5 ** -1.0 <= dht - (time1 - time0)/2 <= 1.0 ** If overflow has not occurred, we also have: ** -0.5 <= hdt - (time1 - time0)/2 <= 0 ** -1.0 <= dht - hdt <= 1.5 ** and since dht - hdt is an integer, we also have: ** -1 <= dht - hdt <= 1 ** or equivalently: ** 0 <= dht - hdt + 1 <= 2 ** In the above analysis, all the operators have ** their exact mathematical semantics, not C semantics. ** However, dht - hdt + 1 is unsigned in C, ** so it need not be compared to zero. */ if (2 < dht - hdt + 1) { /* ** Repair delta overflow. ** ** The following expression contains a second ** rounding, so the result may not be the ** closest to the true answer. This problem ** occurs only with very large differences, ** It's too painful to fix this portably. ** We are not alone in this problem; some C ** compilers round twice when converting ** large unsigned types to small floating ** types, so if time_t is unsigned the ** "delta = dt" above has the same ** double-rounding problem with those ** compilers. */ long_double hibit = ~(UINTMAX_MAX / 2); delta = dt + 2 * hibit; } } } return time1_is_smaller ? -delta : delta; }

Paul Eggert said:
Thanks for the tour of some of the darker corners of C; I didn't realize all the ins and outs of hosts with padding bits.
You're welcome.
/* ** Algorithm courtesy Paul Eggert (eggert@cs.ucla.edu).
No credit?
#ifndef HAVE_STDINT_H #define HAVE_STDINT_H (199901 <= __STDC_VERSION__)
19910L (can you say "16 bits"?)
#if HAVE_STDINT_H #include <stdint.h> #define uintmax uintmax_t
I always feel typedef is safer for this sort of thing.
/* ** Use floating point if there should be no double-rounding error. ** However, avoid long double if it must be wider than needed, ** as it's sometimes much more expensive in these cases ** (e.g., 64-bit sparc). */ if (TYPE_BIT(time_t) <= DBL_MANT_DIG || (TYPE_FLOATING(time_t) && sizeof(time_t) < sizeof(long_double))) {
Using != allows for certain perverse implementations. Not important, though.
if ((TYPE_BIT(time_t) <= LDBL_MANT_DIG && (TYPE_BIT(time_t) == LDBL_MANT_DIG || (TYPE_SIGNED(time_t) && UINTMAX_MAX / 2 < INTMAX_MAX)))
Why write this last test like this? UINTMAX_MAX and INTMAX_MAX are both going to be numbers of the form 2**N-1. Furthermore, UINTMAX_MAX >= INTMAX_MAX. So that last one can be done as UINTMAX_MAX == INTMAX_MAX. It happens again later. -- 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 | |

"Clive D.W. Feather" <clive@demon.net> writes:
No credit?
Sure. How about adding this comment? ** Thanks to Clive D.W. Feather for detailed technical ** advice about hosts with padding bits.
#ifndef HAVE_STDINT_H #define HAVE_STDINT_H (199901 <= __STDC_VERSION__)
19910L (can you say "16 bits"?)
No, 16 bits is a red herring. The expression is used only in preprocessor contexts, so arithmetic must be at least 32 bits (64 with C99) and the "L" is unnecessary. And even if the expression were used outside of preprocessor contexts, __STDC_VERSION__ can't be negative so the expression would still work. Anyway, I'll put the "L" in (if only to forestall pedantic bug reports :-).
if ((TYPE_BIT(time_t) <= LDBL_MANT_DIG && (TYPE_BIT(time_t) == LDBL_MANT_DIG || (TYPE_SIGNED(time_t) && UINTMAX_MAX / 2 < INTMAX_MAX)))
Why write this last test like this? UINTMAX_MAX and INTMAX_MAX are both going to be numbers of the form 2**N-1. Furthermore, UINTMAX_MAX >= INTMAX_MAX. So that last one can be done as UINTMAX_MAX == INTMAX_MAX.
It's just a style thing, as the two expressions are logically equivalent under the constraints that you mention. "UINTMAX_MAX / 2 < INTMAX_MAX" is more closely related to the constraints that were in my head when I wrote the code for the case where uintmax has padding bits and time_t is signed.
participants (3)
-
Clive D. W. Feather
-
Clive D.W. Feather
-
Paul Eggert