[Ada] Reimplement Ada.Numerics.Big_Numbers.Big_Reals.Float_Conversions
authorEric Botcazou <ebotcazou@adacore.com>
Tue, 3 Nov 2020 10:50:40 +0000 (11:50 +0100)
committerPierre-Marie de Rodat <derodat@adacore.com>
Fri, 27 Nov 2020 09:15:48 +0000 (04:15 -0500)
gcc/ada/

* libgnat/a-nbnbre.adb: Remove clauses for System.Img_Real and
add them for System.Unsigned_Types.
(Float_Conversions.To_Big_Real): Reimplement.
(Float_Conversions.From_Big_Real): Likewise.

gcc/ada/libgnat/a-nbnbre.adb

index e5d810e4db5c601b09ed11821eaf7e9282d4b96d..4254b302faff366efcca59c380728f5272d06461 100644 (file)
@@ -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;