libgcc2.c floatdisf is broken

[Alan Modra <amodra at bigpond dot net dot au>]

floatdisf, 64 bit int to float (32 bit IEEE, ie. 23 bit fraction)

If implemented as

  DI -> DF
  DF -> SF

then there is a possible double rounding problem, when the first
conversion is not exact, and the rounding mode is "round to nearest".
Rounding modes of "toward zero", "toward +inf" or "toward -inf" do not
suffer this problem.  The following was written with rs6000.md
floatdisf2 pattern in mind, but the same analysis applies to libgcc2.c
__floatdisf function where the DI -> DF rounding occurs when the low
part of the DImode value is added to the high part.


Examples:

2**60 is representable as a SF number, the next higher SF value being
2**60 + 2**37.  Considering integers in this range, according to the
"round to nearest" rules, values in the range [2**60, 2**60 + 2**36]
should be converted to a result of 2**60, while integers in the range
[2**60 + 2**36 + 1, 2**60 + 2**37] should convert to 2**60 + 2**37.
Note that since powerpc hardware uses "round to nearest even", the
midpoint, 2**60 + 2**36, is converted to the lower value in this case
as the lower value has a zero LSB.  If the lower value had the LSB
set, we'd convert the midpoint up.

Let's see what happens with 2**60 + 2**36 + 1.

The nearest DF value (52 bit fraction) is 2**60 + 2**36 (the next
highest value is 2**60 + 2**36 + 2**8), so the result of the first
conversion is a number that is further rounded down in the DF->SF
conversion, and we end up with an incorrect result of 2**60.  It's
fairly easy to see that integers in the range [2**60 + 2**36 + 1,
2**60 + 2**36 + 2**7] will similarly produce the wrong result.

The next higher SF range is [2**60 + 2**37, 2**60 + 2**38].  Here, the
midpoint value of 2**60 + 3*2**36 converts up to 2**60 + 2**38.

What happens with 2**60 + 3*2**36 - 1?  Well, the DI->DF conversion
results in 2**60 + 3*2**36, which goes to 2**60 + 2**38 on the DF->SF
conversion.  Wrong again.  This will happen for integers in the range
[2**60 + 3*2**36 - 2**7, 2**60 + 3*2**36 - 1].


Details:

As the examples indicate, the problem occurs when the first conversion
results in a number that falls on the midpoint between two SF values,
and the midpoint is rounded away from the original value.
ie. The original integer is a positive number of the form

(A) 0..01XX..X010000000000000000000000000000Y..Y
where XX..X is a string of 22 binary digits, and Y..Y is a
string of 1 or more digits with value 0..01 to 10..0

or

(B) 0..01XX..X101111111111111111111111111111Y..Y
where XX..X is a string of 22 binary digits, and Y..Y is a
string of 1 or more digits with value 10..0 to 11..1

or a negative number with a bit pattern matching the complement of the
above patterns.  Values matching (A) will be incorrectly rounded down,
while values matching (B) will be incorrectly rounded up.


Conclusion:

The attempted fix in libgcc2.c:__floatdisf for the double rounding
problem won't work.  Case (B) cannot be fixed by setting a bit.
I think the following patch will cure case (B) as well as (A), but
haven't properly tested for compliance with IEEE rounding requirements
yet.  ie. I'm not at this point requesting permission to apply the
patch.  It's more a case of inviting comment.  Also, can someone
suggest a test program?

	* libgcc2.c (__floatdisf): Mask off bits below REP_BIT.
	* config/rs6000/rs6000.md (floatdisf2): Rename to floatdisf2_int1.
	(floatdisf2): New define_expand.
	(floatdisf2_int2): Likewise.

Index: gcc/libgcc2.c
===================================================================
RCS file: /cvs/gcc/gcc/gcc/libgcc2.c,v
retrieving revision 1.149
diff -u -p -r1.149 libgcc2.c
--- gcc/libgcc2.c	8 Sep 2002 12:47:26 -0000	1.149
+++ gcc/libgcc2.c	17 Sep 2002 13:42:09 -0000
@@ -1091,7 +1091,10 @@ __floatdisf (DWtype u)
 	     && u < ((DWtype) 1 << DF_SIZE)))
 	{
 	  if ((UDWtype) u & (REP_BIT - 1))
-	    u |= REP_BIT;
+	    {
+	      u &= ~ (REP_BIT - 1);
+	      u |= REP_BIT;
+	    }
 	}
     }
   f = (Wtype) (u >> WORD_SIZE);
Index: gcc/config/rs6000/rs6000.md
===================================================================
RCS file: /cvs/gcc/gcc/gcc/config/rs6000/rs6000.md,v
retrieving revision 1.206
diff -u -p -r1.206 rs6000.md
--- gcc/config/rs6000/rs6000.md	13 Sep 2002 01:40:44 -0000	1.206
+++ gcc/config/rs6000/rs6000.md	17 Sep 2002 14:55:58 -0000
@@ -5890,13 +5890,30 @@
   "fctidz %0,%1"
   [(set_attr "type" "fp")])
 
-;; This only is safe if rounding mode set appropriately.
-(define_insn_and_split "floatdisf2"
+(define_expand "floatdisf2"
+  [(set (match_operand:SF 0 "gpc_reg_operand" "")
+        (float:SF (match_operand:DI 1 "gpc_reg_operand" "")))]
+  "TARGET_POWERPC64 && TARGET_HARD_FLOAT && TARGET_FPRS"
+  "
+{
+  if (!flag_unsafe_math_optimizations)
+    {
+      rtx label = gen_label_rtx ();
+      emit_insn (gen_floatdisf2_int2 (operands[1], label));
+      emit_label (label);
+    }
+  emit_insn (gen_floatdisf2_int1 (operands[0], operands[1]));
+  DONE;
+}")
+
+;; This is not IEEE compliant if rounding mode is "round to nearest".
+;; If the DI->DF conversion is inexact, then it's possible to suffer
+;; from double rounding.
+(define_insn_and_split "floatdisf2_int1"
   [(set (match_operand:SF 0 "gpc_reg_operand" "=f")
         (float:SF (match_operand:DI 1 "gpc_reg_operand" "*f")))
    (clobber (match_scratch:DF 2 "=f"))]
-  "TARGET_POWERPC64 && TARGET_HARD_FLOAT && TARGET_FPRS
-   && flag_unsafe_math_optimizations"
+  "TARGET_POWERPC64 && TARGET_HARD_FLOAT && TARGET_FPRS"
   "#"
   "&& reload_completed"
   [(set (match_dup 2)
@@ -5904,6 +5921,31 @@
    (set (match_dup 0)
         (float_truncate:SF (match_dup 2)))]
   "")
+
+;; Twiddles bits to avoid double rounding.
+;; See libgcc2.c __floatdisf, except that this is somewhat optimized.
+(define_expand "floatdisf2_int2"
+  [(set (match_dup 2) (and:DI (match_operand:DI 0 "" "") (const_int 2047)))
+   (set (match_dup 4) (compare:CC (match_dup 2) (const_int 0)))
+   (set (match_dup 3) (ashiftrt:DI (match_dup 0) (const_int 53)))
+   (set (match_dup 3) (plus:DI (match_dup 3) (const_int 1)))
+   (set (pc) (if_then_else (eq (match_dup 4) (const_int 0))
+			   (label_ref (match_operand:DI 1 "" ""))
+			   (pc)))
+   (set (match_dup 5) (compare:CCUNS (match_dup 3) (const_int 2)))
+   (set (pc) (if_then_else (ltu (match_dup 5) (const_int 0))
+			   (label_ref (match_dup 1))
+			   (pc)))
+   (set (match_dup 0) (xor:DI (match_dup 0) (match_dup 2)))
+   (set (match_dup 0) (ior:DI (match_dup 0) (const_int 2048)))]
+  "TARGET_POWERPC64 && TARGET_HARD_FLOAT && TARGET_FPRS"
+  "
+{
+  operands[2] = gen_reg_rtx (DImode);
+  operands[3] = gen_reg_rtx (DImode);
+  operands[4] = gen_reg_rtx (CCmode);
+  operands[5] = gen_reg_rtx (CCUNSmode);
+}")
 
 ;; Define the DImode operations that can be done in a small number
 ;; of instructions.  The & constraints are to prevent the register

-- 
Alan Modra
IBM OzLabs - Linux Technology Centre