[Bug libc/535] New: difftime double-rounding bug if 64-bit time_t and long double=double

["eggert at gnu dot org" <sourceware-bugzilla at sources dot redhat dot com>]

On hosts with where time_t is 64-bit signed integer, and long double
is IEEE double (i.e., the same type as "double"), glibc difftime
generates answers that are slightly off in some cases.  This bug is
described in the difftime.c source code.  Here's a test case that
illustrates the bug.  This test case is designed to fail on hosts
with the above characteristics, but succeed on other hosts.

  #include <time.h>
  #include <stdlib.h>

  int
  main (void)
  {
    /* Check for double-rounding error when time_t is 64-bit signed int,
       long double is IEEE 754 double with default rounding, time1 =
       9223372036854775807 and time0 = -1536.  Then the true difference
       is 9223372036854777343, which rounds to 9223372036854777856 with
       a total error of 513.  But with the buggy version of difftime,
       time1 - time0 overflows to -9223372036854774273, which rounds to
       -9223372036854774784, and when difftime corrected this by adding
       2**64 = 18446744073709551616 the result was 9223372036854776832,
       which rounds to 9223372036854775808 with a total error of 1535
       instead.  */

    unsigned long long int t1 = 9223372036854775807ull;
    time_t time1 = t1;
    time_t time0 = -1536;
    double dt0 = difftime (time1, time0);
    double dt1 = 9223372036854777343.0;

    if (0 < time1 && t1 == time1 && time0 < 0 && dt0 != dt1)
      exit (1);

    return 0;
  }

This subject came up recently on the time zone mailing list, and Clive
Feather pointed out that the current glibc algorithm has some other
related portability assumptions.  For example, it isn't portable to
some (admittedly weird) C99 hosts that have padding bits (e.g., bits
in integers that do not contribute to their values) in some cases.
Consulting with him on the finer points of C99, I came up with the
following patch.

2004-11-09  Paul Eggert  <eggert@cs.ucla.edu.

	* time/difftime.c: Fix a double-rounding bug on hosts with
	64-bit time_t and long double being IEEE double.  Also, port
	to more valid C99 hosts, even those that have padding bits.
	Don't include <values.h> since it is marked as an obsolescent
	interface.  Include <limits.h>, <float.h>, and <stdint.h> instead.
	(TYPE_BITS, TYPE_FLOATING, TYPE_SIGNED): New macros.
	(subtract): New static function, that works correctly without
	double-rounding, even on hosts with 64-bit time_t.  Also cater
	to hosts with padding bits.
	(__difftime): Use it.  Use DBL_MANT_DIG and LDBL_MANT_DIG to
	determine whether floating types are wide enough: the old
	test (which used sizeof) could in theory report the wrong results
	on hosts with padding bits in floating-point values.

*** libc/time/difftime.c	Thu May  6 21:00:52 2004
--- libc-fixed/time/difftime.c	Mon Nov  8 23:33:46 2004
***************
*** 16,67 ****
     Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
     02111-1307 USA.  */
  
  #include <time.h>
- #include <values.h>
  
  
  /* Return the difference between TIME1 and TIME0.  */
  double
! __difftime (time1, time0)
!      time_t time1;
!      time_t time0;
  {
!   /* Algorithm courtesy Paul Eggert (eggert@twinsun.com).  */
! 
!   time_t delta, hibit;
  
!   if (sizeof (time_t) < sizeof (double))
      return (double) time1 - (double) time0;
!   if (sizeof (time_t) < sizeof (long double))
      return (long double) time1 - (long double) time0;
  
!   if (time1 < time0)
!     return - __difftime (time0, time1);
  
!   /* As much as possible, avoid loss of precision by computing the
!     difference before converting to double.  */
!   delta = time1 - time0;
!   if (delta >= 0)
!     return delta;
! 
!   /* Repair delta overflow.  */
!   hibit = (~ (time_t) 0) << (_TYPEBITS (time_t) - 1);
! 
!   /* The following expression rounds twice, which means the result may not
!      be the closest to the true answer.  For example, suppose time_t is
!      64-bit signed int, long_double is IEEE 754 double with default
!      rounding, time1 = 9223372036854775807 and time0 = -1536.  Then the
!      true difference is 9223372036854777343, which rounds to
!      9223372036854777856 with a total error of 513.  But delta overflows to
!      -9223372036854774273, which rounds to -9223372036854774784, and
!      correcting this by subtracting 2 * (long_double) hibit (i.e. by adding
!      2**64 = 18446744073709551616) yields 9223372036854776832, which rounds
!      to 9223372036854775808 with a total error of 1535 instead.  This
!      problem occurs only with very large differences.  It's too painful to
!      fix this portably.  We are not alone in this problem; many C compilers
!      round twice when converting large unsigned types to small floating
!      types, so if time_t is unsigned the "return delta" above has the same
!      double-rounding problem.  */
!   return delta - 2 * (long double) hibit;
  }
  strong_alias (__difftime, difftime)
--- 16,122 ----
     Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
     02111-1307 USA.  */
  
+ /* Written by Paul Eggert <eggert@cs.ucla.edu>.  */
+ 
  #include <time.h>
  
+ #include <limits.h>
+ #include <float.h>
+ #include <stdint.h>
+ 
+ #define TYPE_BITS(type) (sizeof (type) * CHAR_BIT)
+ #define TYPE_FLOATING(type) ((type) 0.5 == 0.5)
+ #define TYPE_SIGNED(type) ((type) -1 < 0)
+ 
+ /* Return the difference between TIME1 and TIME0, where TIME0 <= TIME1.
+    time_t is known to be an integer type.  */
+ 
+ static double
+ subtract (time_t time1, time_t time0)
+ {
+   if (! TYPE_SIGNED (time_t))
+     return time1 - time0;
+   else
+     {
+       /* Optimize the common special cases where time_t
+ 	 can be converted to uintmax_t without losing information.  */
+       uintmax_t dt = (uintmax_t) time1 - (uintmax_t) time0;
+       double delta = dt;
+ 
+       if (UINTMAX_MAX / 2 < INTMAX_MAX)
+ 	{
+ 	  /* This is a rare host where uintmax_t has padding bits, and possibly
+ 	     information was lost when converting time_t to uintmax_t.
+ 	     Check for overflow by comparing dt/2 to (time1/2 - time0/2).
+ 	     Overflow occurred if they differ by more than a small slop.
+ 	     Thanks to Clive D.W. Feather for detailed technical advice about
+ 	     hosts with padding bits.
+ 
+ 	     In the following code the "h" prefix means half.  By range
+ 	     analysis, we have:
+ 
+                   -0.5 <= ht1 - 0.5*time1 <= 0.5
+                   -0.5 <= ht0 - 0.5*time0 <= 0.5
+                   -1.0 <= dht - 0.5*(time1 - time0) <= 1.0
+ 
+              If overflow has not occurred, we also have:
+ 
+                   -0.5 <= hdt - 0.5*(time1 - time0) <= 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.  */
+ 
+ 	  uintmax_t hdt = dt / 2;
+ 	  time_t ht1 = time1 / 2;
+ 	  time_t ht0 = time0 / 2;
+ 	  time_t dht = ht1 - ht0;
+ 
+ 	  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.  */
+ 
+ 	      delta = dt + 2.0L * (UINTMAX_MAX - UINTMAX_MAX / 2);
+ 	    }
+ 	}
+ 
+       return delta;
+     }
+ }
  
  /* Return the difference between TIME1 and TIME0.  */
  double
! __difftime (time_t time1, time_t time0)
  {
!   /* Convert to double and then subtract if no double-rounding error could
!      result.  */
  
!   if (TYPE_BITS (time_t) <= DBL_MANT_DIG
!       || (TYPE_FLOATING (time_t) && sizeof (time_t) < sizeof (long double)))
      return (double) time1 - (double) time0;
! 
!   /* Likewise for long double.  */
! 
!   if (TYPE_BITS (time_t) <= LDBL_MANT_DIG || TYPE_FLOATING (time_t))
      return (long double) time1 - (long double) time0;
  
!   /* Subtract the smaller integer from the larger, convert the difference to
!      double, and then negate if needed.  */
  
!   return time1 < time0 ? - subtract (time0, time1) : subtract (time1, time0);
  }
  strong_alias (__difftime, difftime)

-- 
           Summary: difftime double-rounding bug if 64-bit time_t and long
                    double=double
           Product: glibc
           Version: 2.3.3
            Status: NEW
          Severity: normal
          Priority: P2
         Component: libc
        AssignedTo: gotom at debian dot or dot jp
        ReportedBy: eggert at gnu dot org
                CC: glibc-bugs at sources dot redhat dot com


http://sources.redhat.com/bugzilla/show_bug.cgi?id=535

------- You are receiving this mail because: -------
You are on the CC list for the bug, or are watching someone who is.