Difftime code

[Paul Eggert]

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 at 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 at 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;
}