[PATCH] Fix undefined behaviour in lround,llround.

[Dave Korn <dave dot korn dot cygwin at gmail dot com>]

    Hi all,

  Back at the start of June we had a report(*) of a bug in the gfortran nint()
function, which I tracked down to a problem in the underlying implementation
of lround() in newlib.  The bug occurs here:

>   else if (exponent_less_1023 < (8 * sizeof (long int)) - 1)
>     {
>       if (exponent_less_1023 >= 52)
>         result = ((long int) msw << (exponent_less_1023 - 20)) | (lsw << (exponent_less_1023 - 52));
>       else
>         {
>           unsigned int tmp = lsw + (0x80000000 >> (exponent_less_1023 - 20));
>           if (tmp < lsw)
>             ++msw;
>           result = ((long int) msw << (exponent_less_1023 - 20)) | (tmp >> (52 - exponent_less_1023));
>         }
>     }
>   else

  When exponent_less_1023 is exactly 52, the subexpression "((long int) msw <<
(exponent_less_1023 - 20))" becomes a left shift by 32.  On systems where a
long int is 32 bits, that's undefined behaviour; the shift amount must always
be less than the size of the left-hand operand (n1256, 6.5.7#3), and in
practice, this means that a standard x86 shift opcode is used, which has the
property of implicitly AND-ing the shift amount with 31; the result of a shift
by 32 bits thus isn't zero as the code expects, but the original unshifted
value of msw, which gets OR-d into completely the wrong place and splots bits
all over the correctly-shifted lsw part of the result.

  In this patch I define a couple of safe shift macros that return zero and
avoid the undefined behaviour when the shift amount is too great.  I've
actually been slightly over-cautious in applying it a couple of places, I
think, but this isn't the fast math code and anyway value-range analysis is a
job for a compiler, which it can do bearing in mind the actual type sizes in
effect during compilation.

newlib/ChangeLog:

	* libm/common/fdlibm.h (SAFE_LEFT_SHIFT): New macro definition.
	(SAFE_RIGHT_SHIFT): Likewise.
	* libm/common/s_llround.c (llround): Use them to avoid undefined
	behaviour.
	* libm/common/s_lround.c (lround): Likewise.

  OK?

    cheers,
      DaveK
-- 
(*) - http://cygwin.com/ml/cygwin/2010-06/threads.html#00315

Index: newlib/libm/common/fdlibm.h
===================================================================
RCS file: /cvs/src/src/newlib/libm/common/fdlibm.h,v
retrieving revision 1.7
diff -p -u -r1.7 fdlibm.h
--- newlib/libm/common/fdlibm.h	17 Nov 2009 22:35:46 -0000	1.7
+++ newlib/libm/common/fdlibm.h	16 Jul 2010 15:31:46 -0000
@@ -361,3 +361,13 @@ do {								\
   sf_u.word = (i);						\
   (d) = sf_u.value;						\
 } while (0)
+
+/* Macros to avoid undefined behaviour that can arise if the amount
+   of a shift is exactly equal to the size of the shifted operand.  */
+
+#define SAFE_LEFT_SHIFT(op,amt)					\
+  (((amt) < 8 * sizeof(op)) ? (op << amt) : 0)
+
+#define SAFE_RIGHT_SHIFT(op,amt)				\
+  (((amt) < 8 * sizeof(op)) ? (op >> amt) : 0)
+
Index: newlib/libm/common/s_llround.c
===================================================================
RCS file: /cvs/src/src/newlib/libm/common/s_llround.c,v
retrieving revision 1.1
diff -p -u -r1.1 s_llround.c
--- newlib/libm/common/s_llround.c	25 Mar 2009 19:13:01 -0000	1.1
+++ newlib/libm/common/s_llround.c	16 Jul 2010 15:31:46 -0000
@@ -49,13 +49,18 @@ llround(double x)
   else if (exponent_less_1023 < (8 * sizeof (long long int)) - 1)
     {
       if (exponent_less_1023 >= 52)
-        result = ((long long int) msw << (exponent_less_1023 - 20)) | (lsw << (exponent_less_1023 - 52));
+        result = SAFE_LEFT_SHIFT ((long long int) msw,
+                                        (exponent_less_1023 - 20))
+                    | SAFE_LEFT_SHIFT (lsw, (exponent_less_1023 - 52));
       else
         {
-          unsigned int tmp = lsw + (0x80000000 >> (exponent_less_1023 - 20));
+          unsigned int tmp = lsw
+                    + SAFE_RIGHT_SHIFT (0x80000000, (exponent_less_1023 - 20));
           if (tmp < lsw)
             ++msw;
-          result = ((long long int) msw << (exponent_less_1023 - 20)) | (tmp >> (52 - exponent_less_1023));
+          result = SAFE_LEFT_SHIFT ((long long int) msw,
+                                        (exponent_less_1023 - 20))
+                    | SAFE_RIGHT_SHIFT (tmp, (52 - exponent_less_1023));
         }
     }
   else
Index: newlib/libm/common/s_lround.c
===================================================================
RCS file: /cvs/src/src/newlib/libm/common/s_lround.c,v
retrieving revision 1.2
diff -p -u -r1.2 s_lround.c
--- newlib/libm/common/s_lround.c	25 Mar 2009 19:13:01 -0000	1.2
+++ newlib/libm/common/s_lround.c	16 Jul 2010 15:31:46 -0000
@@ -90,13 +90,16 @@ ANSI C, POSIX
   else if (exponent_less_1023 < (8 * sizeof (long int)) - 1)
     {
       if (exponent_less_1023 >= 52)
-        result = ((long int) msw << (exponent_less_1023 - 20)) | (lsw << (exponent_less_1023 - 52));
+        result = SAFE_LEFT_SHIFT ((long int) msw, (exponent_less_1023 - 20))
+                | SAFE_LEFT_SHIFT (lsw, (exponent_less_1023 - 52));
       else
         {
-          unsigned int tmp = lsw + (0x80000000 >> (exponent_less_1023 - 20));
+          unsigned int tmp = lsw
+                    + SAFE_RIGHT_SHIFT (0x80000000, (exponent_less_1023 - 20));
           if (tmp < lsw)
             ++msw;
-          result = ((long int) msw << (exponent_less_1023 - 20)) | (tmp >> (52 - exponent_less_1023));
+          result = SAFE_LEFT_SHIFT ((long int) msw, (exponent_less_1023 - 20))
+                    | SAFE_RIGHT_SHIFT (tmp, (52 - exponent_less_1023));
         }
     }
   else