diff options
Diffstat (limited to 'gcc/ada/libgnat/s-valrea.adb')
| -rw-r--r-- | gcc/ada/libgnat/s-valrea.adb | 198 |
1 files changed, 159 insertions, 39 deletions
diff --git a/gcc/ada/libgnat/s-valrea.adb b/gcc/ada/libgnat/s-valrea.adb index 0ac3846..bc5465c 100644 --- a/gcc/ada/libgnat/s-valrea.adb +++ b/gcc/ada/libgnat/s-valrea.adb @@ -6,7 +6,7 @@ -- -- -- B o d y -- -- -- --- Copyright (C) 1992-2020, Free Software Foundation, Inc. -- +-- Copyright (C) 1992-2021, Free Software Foundation, Inc. -- -- -- -- GNAT is free software; you can redistribute it and/or modify it under -- -- terms of the GNU General Public License as published by the Free Soft- -- @@ -29,6 +29,7 @@ -- -- ------------------------------------------------------------------------------ +with System.Double_Real; with System.Float_Control; with System.Unsigned_Types; use System.Unsigned_Types; with System.Val_Util; use System.Val_Util; @@ -46,7 +47,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 +77,12 @@ 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); + package Double_Real is new System.Double_Real (Num); + use type Double_Real.Double_T; + + subtype Double_T is Double_Real.Double_T; + -- The double floating-point type + function Integer_to_Real (Str : String; Val : Uns; @@ -85,6 +92,9 @@ package body System.Val_Real is Minus : Boolean) return Num; -- Convert the real value from integer to real representation + function Large_Powten (Exp : Natural) return Double_T; + -- Return 10.0**Exp as a double number, where Exp > Maxpow + --------------------- -- Integer_to_Real -- --------------------- @@ -110,9 +120,8 @@ 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; + D_Val : Double_T; S : Integer := Scale; begin @@ -125,75 +134,151 @@ package body System.Val_Real is System.Float_Control.Reset; end if; - -- Do the conversion + -- Take into account the extra digit, i.e. do the two computations - R_Val := Num (Val); + -- (1) R_Val := R_Val * Num (B) + Num (Extra) + -- (2) S := S - 1 - -- Take into account the extra digit, if need be. In this case, the - -- three operands are exact, so using an FMA would be ideal. + -- 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 := Num (Val); + DS : Double_T; + + 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 + if Acc = 0.0 then + Acc := Fac; + else + DS := Double_Real.Quick_Two_Sum (Fac, Acc); + Acc := DS.Hi; + Err := Err + DS.Lo; + end if; + 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 + + D_Val := Double_Real.Quick_Two_Sum (Acc, Err + Num (Extra)); + + S := S - 1; + end; + + -- Or else, if the Extra digit is zero, do the exact conversion + + elsif Need_Extra then + D_Val := Double_Real.To_Double (Num (Val)); + + -- Otherwise, the value contains more bits than the mantissa so do the + -- conversion in two steps. + + else + declare + Mask : constant Uns := 2**(Uns'Size - Num'Machine_Mantissa) - 1; + Hi : constant Uns := Val and not Mask; + Lo : constant Uns := Val and Mask; + + begin + if Hi = 0 then + D_Val := Double_Real.To_Double (Num (Lo)); + else + D_Val := Double_Real.Quick_Two_Sum (Num (Hi), Num (Lo)); + end if; + end; end if; - -- Compute the final value + -- Compute the final value by applying the scaling, if any + + if Val = 0 or else S = 0 then + R_Val := Double_Real.To_Single (D_Val); - if R_Val /= 0.0 and then S /= 0 then + else case Base is -- If the base is a power of two, we use the efficient Scaling -- attribute with an overflow check, if it is not 2, to catch -- ludicrous exponents that would result in an infinity or zero. when 2 => - R_Val := Num'Scaling (R_Val, S); + R_Val := Num'Scaling (Double_Real.To_Single (D_Val), S); when 4 => if Integer'First / 2 <= S and then S <= Integer'Last / 2 then S := S * 2; end if; - R_Val := Num'Scaling (R_Val, S); + R_Val := Num'Scaling (Double_Real.To_Single (D_Val), S); when 8 => if Integer'First / 3 <= S and then S <= Integer'Last / 3 then S := S * 3; end if; - R_Val := Num'Scaling (R_Val, S); + R_Val := Num'Scaling (Double_Real.To_Single (D_Val), S); when 16 => if Integer'First / 4 <= S and then S <= Integer'Last / 4 then S := S * 4; end if; - R_Val := Num'Scaling (R_Val, S); + R_Val := Num'Scaling (Double_Real.To_Single (D_Val), S); - -- If the base is 10, we use a table of powers for accuracy's sake + -- If the base is 10, use a double implementation for the sake + -- of accuracy, to be removed when exponentiation is improved. + + -- When the exponent is positive, we can do the computation + -- directly because, if the exponentiation overflows, then + -- the final value overflows as well. But when the exponent + -- is negative, we may need to do it in two steps to avoid + -- an artificial underflow. when 10 => declare - Powten : constant array (0 .. Maxpow) of Num; + Powten : constant array (0 .. Maxpow) of Double_T; pragma Import (Ada, Powten); for Powten'Address use Powten_Address; begin if S > 0 then - while S > Maxpow loop - R_Val := R_Val * Powten (Maxpow); - S := S - Maxpow; - end loop; - - R_Val := R_Val * Powten (S); + if S <= Maxpow then + D_Val := D_Val * Powten (S); + else + D_Val := D_Val * Large_Powten (S); + end if; else - while S < -Maxpow loop - R_Val := R_Val / Powten (Maxpow); - S := S + Maxpow; - end loop; - - R_Val := R_Val / Powten (-S); + if S < -Maxexp then + D_Val := D_Val / Large_Powten (Maxexp); + S := S + Maxexp; + end if; + + if S >= -Maxpow then + D_Val := D_Val / Powten (-S); + else + D_Val := D_Val / Large_Powten (-S); + end if; end if; + + R_Val := Double_Real.To_Single (D_Val); end; -- Implementation for other bases with exponentiation @@ -205,17 +290,24 @@ 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 := Double_Real.To_Single (D_Val); - 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; @@ -228,6 +320,34 @@ package body System.Val_Real is when Constraint_Error => Bad_Value (Str); end Integer_to_Real; + ------------------ + -- Large_Powten -- + ------------------ + + function Large_Powten (Exp : Natural) return Double_T is + Powten : constant array (0 .. Maxpow) of Double_T; + pragma Import (Ada, Powten); + for Powten'Address use Powten_Address; + + R : Double_T; + E : Natural; + + begin + pragma Assert (Exp > Maxpow); + + R := Powten (Maxpow); + E := Exp - Maxpow; + + while E > Maxpow loop + R := R * Powten (Maxpow); + E := E - Maxpow; + end loop; + + R := R * Powten (E); + + return R; + end Large_Powten; + --------------- -- Scan_Real -- --------------- |