path: root/gcc/ada/libgnat/a-nbnbre.adb

------------------------------------------------------------------------------
--                                                                          --
--                         GNAT RUN-TIME COMPONENTS                         --
--                                                                          --
--                   ADA.NUMERICS.BIG_NUMBERS.BIG_REALS                     --
--                                                                          --
--                                 B o d y                                  --
--                                                                          --
--             Copyright (C) 2019-2023, 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- --
-- ware  Foundation;  either version 3,  or (at your option) any later ver- --
-- sion.  GNAT is distributed in the hope that it will be useful, but WITH- --
-- OUT ANY WARRANTY;  without even the  implied warranty of MERCHANTABILITY --
-- or FITNESS FOR A PARTICULAR PURPOSE.                                     --
--                                                                          --
-- As a special exception under Section 7 of GPL version 3, you are granted --
-- additional permissions described in the GCC Runtime Library Exception,   --
-- version 3.1, as published by the Free Software Foundation.               --
--                                                                          --
-- You should have received a copy of the GNU General Public License and    --
-- a copy of the GCC Runtime Library Exception along with this program;     --
-- see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see    --
-- <http://www.gnu.org/licenses/>.                                          --
--                                                                          --
-- GNAT was originally developed  by the GNAT team at  New York University. --
-- Extensive contributions were provided by Ada Core Technologies Inc.      --
--                                                                          --
------------------------------------------------------------------------------

with System.Unsigned_Types; use System.Unsigned_Types;

package body Ada.Numerics.Big_Numbers.Big_Reals is

   use Big_Integers;

   procedure Normalize (Arg : in out Big_Real);
   --  Normalize Arg by ensuring that Arg.Den is always positive and that
   --  Arg.Num and Arg.Den always have a GCD of 1.

   --------------
   -- Is_Valid --
   --------------

   function Is_Valid (Arg : Big_Real) return Boolean is
     (Is_Valid (Arg.Num) and Is_Valid (Arg.Den));

   ---------
   -- "/" --
   ---------

   function "/" (Num, Den : Valid_Big_Integer) return Valid_Big_Real is
      Result : Big_Real;
   begin
      if Den = To_Big_Integer (0) then
         raise Constraint_Error with "divide by zero";
      end if;

      Result.Num := Num;
      Result.Den := Den;
      Normalize (Result);
      return Result;
   end "/";

   ---------------
   -- Numerator --
   ---------------

   function Numerator (Arg : Valid_Big_Real) return Valid_Big_Integer is
     (Arg.Num);

   -----------------
   -- Denominator --
   -----------------

   function Denominator (Arg : Valid_Big_Real) return Big_Positive is
     (Arg.Den);

   ---------
   -- "=" --
   ---------

   function "=" (L, R : Valid_Big_Real) return Boolean is
     (L.Num = R.Num and then L.Den = R.Den);

   ---------
   -- "<" --
   ---------

   function "<" (L, R : Valid_Big_Real) return Boolean is
     (L.Num * R.Den < R.Num * L.Den);
   --  The denominator is guaranteed to be positive since Normalized is
   --  always called when constructing a Valid_Big_Real

   ----------
   -- "<=" --
   ----------

   function "<=" (L, R : Valid_Big_Real) return Boolean is (not (R < L));

   ---------
   -- ">" --
   ---------

   function ">" (L, R : Valid_Big_Real) return Boolean is (R < L);

   ----------
   -- ">=" --
   ----------

   function ">=" (L, R : Valid_Big_Real) return Boolean is (not (L < R));

   -----------------------
   -- Float_Conversions --
   -----------------------

   package body Float_Conversions is

      package Conv is new
        Big_Integers.Unsigned_Conversions (Long_Long_Unsigned);

      -----------------
      -- 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

         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
         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;

         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

         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
         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;

   -----------------------
   -- Fixed_Conversions --
   -----------------------

   package body Fixed_Conversions is

      package Float_Aux is new Float_Conversions (Long_Float);

      subtype LLLI is Long_Long_Long_Integer;
      subtype LLLU is Long_Long_Long_Unsigned;

      Too_Large : constant Boolean :=
                    Num'Small_Numerator > LLLU'Last
                      or else Num'Small_Denominator > LLLU'Last;
      --  True if the Small is too large for Long_Long_Long_Unsigned, in which
      --  case we convert to/from Long_Float as an intermediate step.

      package Conv_I is new Big_Integers.Signed_Conversions (LLLI);
      package Conv_U is new Big_Integers.Unsigned_Conversions (LLLU);

      -----------------
      -- To_Big_Real --
      -----------------

      --  We just compute V * N / D where V is the mantissa value of the fixed
      --  point number, and N resp. D is the numerator resp. the denominator of
      --  the Small of the fixed-point type.

      function To_Big_Real (Arg : Num) return Valid_Big_Real is
         N, D, V : Big_Integer;

      begin
         if Too_Large then
            return Float_Aux.To_Big_Real (Long_Float (Arg));
         end if;

         N := Conv_U.To_Big_Integer (Num'Small_Numerator);
         D := Conv_U.To_Big_Integer (Num'Small_Denominator);
         V := Conv_I.To_Big_Integer (LLLI'Integer_Value (Arg));

         return V * N / D;
      end To_Big_Real;

      -------------------
      -- From_Big_Real --
      -------------------

      --  We first compute A / B = Arg * D / N where N resp. D is the numerator
      --  resp. the denominator of the Small of the fixed-point type. Then we
      --  divide A by B and convert the result to the mantissa value.

      function From_Big_Real (Arg : Big_Real) return Num is
         N, D, A, B, Q, X : Big_Integer;

      begin
         if Too_Large then
            return Num (Float_Aux.From_Big_Real (Arg));
         end if;

         N := Conv_U.To_Big_Integer (Num'Small_Numerator);
         D := Conv_U.To_Big_Integer (Num'Small_Denominator);
         A := Numerator (Arg) * D;
         B := Denominator (Arg) * N;

         Q := A / B;

         --  Round to nearest, ties to away, by comparing twice the remainder

         X := (A - Q * B) * To_Big_Integer (2);

         if X >= B then
            Q := Q + To_Big_Integer (1);

         elsif X <= -B then
            Q := Q - To_Big_Integer (1);
         end if;

         return Num'Fixed_Value (Conv_I.From_Big_Integer (Q));
      end From_Big_Real;

   end Fixed_Conversions;

   ---------------
   -- To_String --
   ---------------

   function To_String
     (Arg  : Valid_Big_Real;
      Fore : Field := 2;
      Aft  : Field := 3;
      Exp  : Field := 0) return String
   is
      Zero : constant Big_Integer := To_Big_Integer (0);
      Ten  : constant Big_Integer := To_Big_Integer (10);

      function Leading_Padding
        (Str        : String;
         Min_Length : Field;
         Char       : Character := ' ') return String;
      --  Return padding of Char concatenated with Str so that the resulting
      --  string is at least Min_Length long.

      function Trailing_Padding
        (Str    : String;
         Length : Field;
         Char   : Character := '0') return String;
      --  Return Str with trailing Char removed, and if needed either
      --  truncated or concatenated with padding of Char so that the resulting
      --  string is Length long.

      function Image (N : Natural) return String;
      --  Return image of N, with no leading space.

      function Numerator_Image
        (Num   : Big_Integer;
         After : Natural) return String;
      --  Return image of Num as a float value with After digits after the "."
      --  and taking Fore, Aft, Exp into account.

      -----------
      -- Image --
      -----------

      function Image (N : Natural) return String is
         S : constant String := Natural'Image (N);
      begin
         return S (2 .. S'Last);
      end Image;

      ---------------------
      -- Leading_Padding --
      ---------------------

      function Leading_Padding
        (Str        : String;
         Min_Length : Field;
         Char       : Character := ' ') return String is
      begin
         if Str = "" then
            return Leading_Padding ("0", Min_Length, Char);
         else
            return [1 .. Integer'Max (Integer (Min_Length) - Str'Length, 0)
                           => Char] & Str;
         end if;
      end Leading_Padding;

      ----------------------
      -- Trailing_Padding --
      ----------------------

      function Trailing_Padding
        (Str    : String;
         Length : Field;
         Char   : Character := '0') return String is
      begin
         if Str'Length > 0 and then Str (Str'Last) = Char then
            for J in reverse Str'Range loop
               if Str (J) /= '0' then
                  return Trailing_Padding
                    (Str (Str'First .. J), Length, Char);
               end if;
            end loop;
         end if;

         if Str'Length >= Length then
            return Str (Str'First .. Str'First + Length - 1);
         else
            return Str &
              [1 .. Integer'Max (Integer (Length) - Str'Length, 0)
                      => Char];
         end if;
      end Trailing_Padding;

      ---------------------
      -- Numerator_Image --
      ---------------------

      function Numerator_Image
        (Num   : Big_Integer;
         After : Natural) return String
      is
         Tmp   : constant String := To_String (Num);
         Str   : constant String (1 .. Tmp'Last - 1) := Tmp (2 .. Tmp'Last);
         Index : Integer;

      begin
         if After = 0 then
            return Leading_Padding (Str, Fore) & "."
                   & Trailing_Padding ("0", Aft);
         else
            Index := Str'Last - After;

            if Index < 0 then
               return Leading_Padding ("0", Fore)
                 & "."
                 & Trailing_Padding ([1 .. -Index => '0'] & Str, Aft)
                 & (if Exp = 0 then "" else "E+" & Image (Natural (Exp)));
            else
               return Leading_Padding (Str (Str'First .. Index), Fore)
                 & "."
                 & Trailing_Padding (Str (Index + 1 .. Str'Last), Aft)
                 & (if Exp = 0 then "" else "E+" & Image (Natural (Exp)));
            end if;
         end if;
      end Numerator_Image;

   begin
      if Arg.Num < Zero then
         declare
            Str : String := To_String (-Arg, Fore, Aft, Exp);
         begin
            if Str (1) = ' ' then
               for J in 1 .. Str'Last - 1 loop
                  if Str (J + 1) /= ' ' then
                     Str (J) := '-';
                     exit;
                  end if;
               end loop;

               return Str;
            else
               return '-' & Str;
            end if;
         end;
      else
         --  Compute Num * 10^Aft so that we get Aft significant digits
         --  in the integer part (rounded) to display.

         return Numerator_Image
           ((Arg.Num * Ten ** Aft) / Arg.Den, After => Exp + Aft);
      end if;
   end To_String;

   -----------------
   -- From_String --
   -----------------

   function From_String (Arg : String) return Valid_Big_Real is
      Ten   : constant Big_Integer := To_Big_Integer (10);
      Frac  : Big_Integer;
      Exp   : Integer := 0;
      Pow   : Natural := 0;
      Index : Natural := 0;
      Last  : Natural := Arg'Last;

   begin
      for J in reverse Arg'Range loop
         if Arg (J) in 'e' | 'E' then
            if Last /= Arg'Last then
               raise Constraint_Error with "multiple exponents specified";
            end if;

            Last := J - 1;
            Exp := Integer'Value (Arg (J + 1 .. Arg'Last));
            Pow := 0;

         elsif Arg (J) = '.' then
            Index := J - 1;
            exit;
         elsif Arg (J) /= '_' then
            Pow := Pow + 1;
         end if;
      end loop;

      if Index = 0 then
         raise Constraint_Error with "invalid real value";
      end if;

      declare
         Result : Big_Real;
      begin
         Result.Den := Ten ** Pow;
         Result.Num := From_String (Arg (Arg'First .. Index)) * Result.Den;
         Frac := From_String (Arg (Index + 2 .. Last));

         if Result.Num < To_Big_Integer (0) then
            Result.Num := Result.Num - Frac;
         else
            Result.Num := Result.Num + Frac;
         end if;

         if Exp > 0 then
            Result.Num := Result.Num * Ten ** Exp;
         elsif Exp < 0 then
            Result.Den := Result.Den * Ten ** (-Exp);
         end if;

         Normalize (Result);
         return Result;
      end;
   end From_String;

   --------------------------
   -- From_Quotient_String --
   --------------------------

   function From_Quotient_String (Arg : String) return Valid_Big_Real is
      Index : Natural := 0;
   begin
      for J in Arg'First + 1 .. Arg'Last - 1 loop
         if Arg (J) = '/' then
            Index := J;
            exit;
         end if;
      end loop;

      if Index = 0 then
         raise Constraint_Error with "no quotient found";
      end if;

      return Big_Integers.From_String (Arg (Arg'First .. Index - 1)) /
        Big_Integers.From_String (Arg (Index + 1 .. Arg'Last));
   end From_Quotient_String;

   ---------------
   -- Put_Image --
   ---------------

   procedure Put_Image (S : in out Root_Buffer_Type'Class; V : Big_Real) is
      --  This is implemented in terms of To_String. It might be more elegant
      --  and more efficient to do it the other way around, but this is the
      --  most expedient implementation for now.
   begin
      Strings.Text_Buffers.Put_UTF_8 (S, To_String (V));
   end Put_Image;

   ---------
   -- "+" --
   ---------

   function "+" (L : Valid_Big_Real) return Valid_Big_Real is
      Result : Big_Real;
   begin
      Result.Num := L.Num;
      Result.Den := L.Den;
      return Result;
   end "+";

   ---------
   -- "-" --
   ---------

   function "-" (L : Valid_Big_Real) return Valid_Big_Real is
     (Num => -L.Num, Den => L.Den);

   -----------
   -- "abs" --
   -----------

   function "abs" (L : Valid_Big_Real) return Valid_Big_Real is
     (Num => abs L.Num, Den => L.Den);

   ---------
   -- "+" --
   ---------

   function "+" (L, R : Valid_Big_Real) return Valid_Big_Real is
      Result : Big_Real;
   begin
      Result.Num := L.Num * R.Den + R.Num * L.Den;
      Result.Den := L.Den * R.Den;
      Normalize (Result);
      return Result;
   end "+";

   ---------
   -- "-" --
   ---------

   function "-" (L, R : Valid_Big_Real) return Valid_Big_Real is
      Result : Big_Real;
   begin
      Result.Num := L.Num * R.Den - R.Num * L.Den;
      Result.Den := L.Den * R.Den;
      Normalize (Result);
      return Result;
   end "-";

   ---------
   -- "*" --
   ---------

   function "*" (L, R : Valid_Big_Real) return Valid_Big_Real is
      Result : Big_Real;
   begin
      Result.Num := L.Num * R.Num;
      Result.Den := L.Den * R.Den;
      Normalize (Result);
      return Result;
   end "*";

   ---------
   -- "/" --
   ---------

   function "/" (L, R : Valid_Big_Real) return Valid_Big_Real is
      Result : Big_Real;
   begin
      Result.Num := L.Num * R.Den;
      Result.Den := L.Den * R.Num;
      Normalize (Result);
      return Result;
   end "/";

   ----------
   -- "**" --
   ----------

   function "**" (L : Valid_Big_Real; R : Integer) return Valid_Big_Real is
      Result : Big_Real;
   begin
      if R = 0 then
         Result.Num := To_Big_Integer (1);
         Result.Den := To_Big_Integer (1);
      else
         if R < 0 then
            Result.Num := L.Den ** (-R);
            Result.Den := L.Num ** (-R);
         else
            Result.Num := L.Num ** R;
            Result.Den := L.Den ** R;
         end if;

         Normalize (Result);
      end if;

      return Result;
   end "**";

   ---------
   -- Min --
   ---------

   function Min (L, R : Valid_Big_Real) return Valid_Big_Real is
     (if L < R then L else R);

   ---------
   -- Max --
   ---------

   function Max (L, R : Valid_Big_Real) return Valid_Big_Real is
     (if L > R then L else R);

   ---------------
   -- Normalize --
   ---------------

   procedure Normalize (Arg : in out Big_Real) is
      Zero : constant Big_Integer := To_Big_Integer (0);
   begin
      if Arg.Den < Zero then
         Arg.Num := -Arg.Num;
         Arg.Den := -Arg.Den;
      end if;

      if Arg.Num = Zero then
         Arg.Den := To_Big_Integer (1);
      else
         declare
            GCD : constant Big_Integer :=
              Greatest_Common_Divisor (Arg.Num, Arg.Den);
         begin
            Arg.Num := Arg.Num / GCD;
            Arg.Den := Arg.Den / GCD;
         end;
      end if;
   end Normalize;

end Ada.Numerics.Big_Numbers.Big_Reals;