This overcomes the lack of fused multiply-add instruction on the x87
FPU by doing an iterated addition with exact error handling for the
last digit taken into account for the mantissa.
Tested on x86_64-pc-linux-gnu, committed on trunk
gcc/ada/
* libgnat/s-valrea.adb (Fast2Sum): New function.
(Integer_to_Real): Use it in an iterated addition with exact
error handling for the case where an extra digit is needed.
Move local variable now only used in the exponentiation case.diff --git a/gcc/ada/libgnat/s-valrea.adb b/gcc/ada/libgnat/s-valrea.adb
--- a/gcc/ada/libgnat/s-valrea.adb
+++ b/gcc/ada/libgnat/s-valrea.adb
@@ -46,7 +46,7 @@ package body System.Val_Real is
-- If the mantissa of the floating-point type is almost as large as the
-- unsigned type, we do not have enough space for an extra digit in the
-- unsigned type so we handle the extra digit separately, at the cost of
- -- a potential roundoff error.
+ -- a bit more work in Integer_to_Real.
Precision_Limit : constant Uns :=
(if Need_Extra then 2**Num'Machine_Mantissa - 1 else 2**Uns'Size - 1);
@@ -76,6 +76,10 @@ package body System.Val_Real is
7 => 5836, 8 => 5461, 9 => 5168, 10 => 4932, 11 => 4736,
12 => 4570, 13 => 4427, 14 => 4303, 15 => 4193, 16 => 4095);
+ function Fast2Sum (A, B : Num; Err : in out Num) return Num;
+ -- This is the classical Fast2Sum function assuming round to nearest,
+ -- with the error accumulated into Err.
+
function Integer_to_Real
(Str : String;
Val : Uns;
@@ -85,6 +89,25 @@ package body System.Val_Real is
Minus : Boolean) return Num;
-- Convert the real value from integer to real representation
+ --------------
+ -- Fast2Sum --
+ --------------
+
+ function Fast2Sum (A, B : Num; Err : in out Num) return Num is
+ S, Z : Num;
+
+ begin
+ pragma Assert (abs (A) >= abs (B));
+
+ S := A + B;
+ Z := S - A;
+ Z := B - Z;
+
+ Err := Err + Z;
+
+ return S;
+ end Fast2Sum;
+
---------------------
-- Integer_to_Real --
---------------------
@@ -110,8 +133,6 @@ package body System.Val_Real is
else raise Program_Error);
-- Maximum exponent of the base that can fit in Num
- B : constant Num := Num (Base);
-
R_Val : Num;
S : Integer := Scale;
@@ -129,12 +150,53 @@ package body System.Val_Real is
R_Val := Num (Val);
- -- Take into account the extra digit, if need be. In this case, the
- -- three operands are exact, so using an FMA would be ideal.
+ -- Take into account the extra digit, i.e. do the two computations
+
+ -- (1) R_Val := R_Val * Num (B) + Num (Extra)
+ -- (2) S := S - 1
+
+ -- In the first, the three operands are exact, so using an FMA would
+ -- be ideal, but we are most likely running on the x87 FPU, hence we
+ -- may not have one. That is why we turn the multiplication into an
+ -- iterated addition with exact error handling, so that we can do a
+ -- single rounding at the end.
if Need_Extra and then Extra > 0 then
- R_Val := R_Val * B + Num (Extra);
- S := S - 1;
+ declare
+ B : Unsigned := Base;
+
+ Acc : Num := 0.0;
+ Err : Num := 0.0;
+ Fac : Num := R_Val;
+
+ begin
+ loop
+ -- If B is odd, add one factor. Note that the accumulator is
+ -- never larger than the factor at this point (it is in fact
+ -- never larger than the factor minus the initial value).
+
+ if B rem 2 /= 0 then
+ Acc := (if Acc = 0.0 then Fac else Fast2Sum (Fac, Acc, Err));
+ exit when B = 1;
+ end if;
+
+ -- Now B is (morally) even, halve it and double the factor,
+ -- which is always an exact operation.
+
+ B := B / 2;
+ Fac := Fac * 2.0;
+ end loop;
+
+ -- Add Extra to the error, which are both small integers
+
+ Err := Err + Num (Extra);
+
+ -- Acc + Err is the exact result before rounding
+
+ R_Val := Acc + Err;
+
+ S := S - 1;
+ end;
end if;
-- Compute the final value
@@ -207,17 +269,23 @@ package body System.Val_Real is
-- an artificial underflow.
when others =>
- if S > 0 then
- R_Val := R_Val * B ** S;
+ declare
+ B : constant Num := Num (Base);
- else
- if S < -Maxexp then
- R_Val := R_Val / B ** Maxexp;
- S := S + Maxexp;
- end if;
+ begin
- R_Val := R_Val / B ** (-S);
- end if;
+ if S > 0 then
+ R_Val := R_Val * B ** S;
+
+ else
+ if S < -Maxexp then
+ R_Val := R_Val / B ** Maxexp;
+ S := S + Maxexp;
+ end if;
+
+ R_Val := R_Val / B ** (-S);
+ end if;
+ end;
end case;
end if;
