diff --git a/gcc/ada/libgnat/a-nbnbre.adb b/gcc/ada/libgnat/a-nbnbre.adb --- a/gcc/ada/libgnat/a-nbnbre.adb +++ b/gcc/ada/libgnat/a-nbnbre.adb @@ -30,7 +30,7 @@ ------------------------------------------------------------------------------ with Ada.Strings.Text_Output.Utils; -with System.Img_Real; use System.Img_Real; +with System.Unsigned_Types; use System.Unsigned_Types; package body Ada.Numerics.Big_Numbers.Big_Reals is @@ -122,25 +122,184 @@ package body Ada.Numerics.Big_Numbers.Big_Reals is -- To_Big_Real -- ----------------- + -- We get the fractional representation of the floating-point number by + -- multiplying Num'Fraction by 2.0**M, with M the size of the mantissa, + -- which gives zero or a number in the range [2.0**(M-1)..2.0**M), which + -- means that it is an integer N of M bits. The floating-point number is + -- thus equal to N / 2**(M-E) where E is its Num'Exponent. + function To_Big_Real (Arg : Num) return Valid_Big_Real is - S : String (1 .. Max_Real_Image_Length); - P : Natural := 0; + + package Conv is new + Big_Integers.Unsigned_Conversions (Long_Long_Unsigned); + + A : constant Num'Base := abs (Arg); + E : constant Integer := Num'Exponent (A); + F : constant Num'Base := Num'Fraction (A); + M : constant Natural := Num'Machine_Mantissa; + + N, D : Big_Integer; + begin - -- Use Long_Long_Unsigned'Width - 1 digits = 20 which is sufficient - -- for the largest floating point format. + pragma Assert (Num'Machine_Radix = 2); + -- This implementation does not handle radix 16 + + pragma Assert (M <= 64); + -- This implementation handles only 80-bit IEEE Extended or smaller + + N := Conv.To_Big_Integer (Long_Long_Unsigned (F * 2.0**M)); + + -- If E is smaller than M, the denominator is 2**(M-E) + + if E < M then + D := To_Big_Integer (2) ** (M - E); + + -- Or else, if E is larger than M, multiply the numerator by 2**(E-M) + + elsif E > M then + N := N * To_Big_Integer (2) ** (E - M); + D := To_Big_Integer (1); + + -- Otherwise E is equal to M and the result is just N + + else + D := To_Big_Integer (1); + end if; - Set_Image_Real - (Long_Long_Float (Arg), S, P, Fore => 1, Aft => 20, Exp => 5); - return From_String (S (1 .. P)); + return (if Arg >= 0.0 then N / D else -N / D); end To_Big_Real; ------------------- -- From_Big_Real -- ------------------- + -- We get the (Frac, Exp) representation of the real number by finding + -- the exponent E such that it lies in the range [2.0**(E-1)..2.0**E), + -- multiplying the number by 2.0**(M-E) with M the size of the mantissa, + -- and converting the result to integer N in the range [2**(M-1)..2**M) + -- with rounding to nearest, ties to even, and finally call Num'Compose. + -- This does not apply to the zero, for which we return 0.0 early. + function From_Big_Real (Arg : Big_Real) return Num is + + package Conv is new + Big_Integers.Unsigned_Conversions (Long_Long_Unsigned); + + M : constant Natural := Num'Machine_Mantissa; + One : constant Big_Real := To_Real (1); + Two : constant Big_Real := To_Real (2); + Half : constant Big_Real := One / Two; + TwoI : constant Big_Integer := To_Big_Integer (2); + + function Log2_Estimate (V : Big_Real) return Natural; + -- Return an integer not larger than Log2 (V) for V >= 1.0 + + function Minus_Log2_Estimate (V : Big_Real) return Natural; + -- Return an integer not larger than -Log2 (V) for V < 1.0 + + ------------------- + -- Log2_Estimate -- + ------------------- + + function Log2_Estimate (V : Big_Real) return Natural is + Log : Natural := 1; + Pow : Big_Real := Two; + + begin + while V >= Pow loop + Pow := Pow * Pow; + Log := Log + Log; + end loop; + + return Log / 2; + end Log2_Estimate; + + ------------------------- + -- Minus_Log2_Estimate -- + ------------------------- + + function Minus_Log2_Estimate (V : Big_Real) return Natural is + Log : Natural := 1; + Pow : Big_Real := Half; + + begin + while V <= Pow loop + Pow := Pow * Pow; + Log := Log + Log; + end loop; + + return Log / 2; + end Minus_Log2_Estimate; + + -- Local variables + + V : Big_Real := abs (Arg); + E : Integer := 0; + L : Integer; + + A, B, Q, X : Big_Integer; + N : Long_Long_Unsigned; + R : Num'Base; + begin - return Num'Value (To_String (Arg)); + pragma Assert (Num'Machine_Radix = 2); + -- This implementation does not handle radix 16 + + pragma Assert (M <= 64); + -- This implementation handles only 80-bit IEEE Extended or smaller + + -- Protect from degenerate case + + if Numerator (V) = To_Big_Integer (0) then + return 0.0; + end if; + + -- Use a binary search to compute exponent E + + while V < Half loop + L := Minus_Log2_Estimate (V); + V := V * (Two ** L); + E := E - L; + end loop; + + -- The dissymetry with above is expected since we go below 2 + + while V >= One loop + L := Log2_Estimate (V) + 1; + V := V / (Two ** L); + E := E + L; + end loop; + + -- The multiplication by 2.0**(-E) has already been done in the loops + + V := V * To_Big_Real (TwoI ** M); + + -- Now go into the integer domain and divide + + A := Numerator (V); + B := Denominator (V); + + Q := A / B; + N := Conv.From_Big_Integer (Q); + + -- Round to nearest, ties to even, by comparing twice the remainder + + X := (A - Q * B) * TwoI; + + if X > B or else (X = B and then (N mod 2) = 1) then + N := N + 1; + + -- If the adjusted quotient overflows the mantissa, scale up + + if N = 2**M then + N := 1; + E := E + 1; + end if; + end if; + + R := Num'Compose (Num'Base (N), E); + + return (if Numerator (Arg) >= To_Big_Integer (0) then R else -R); end From_Big_Real; end Float_Conversions;