------------------------------------------------------------------------------
--                                                                          --
--                         GNAT COMPILER COMPONENTS                         --
--                                                                          --
--                    S Y S T E M . D O U B L E _ R E A L                   --
--                                                                          --
--                                 B o d y                                  --
--                                                                          --
--            Copyright (C) 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- --
-- 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.      --
--                                                                          --
------------------------------------------------------------------------------

package body System.Double_Real is

   function Is_NaN (N : Num) return Boolean is (N /= N);
   --  Return True if N is a NaN

   function Is_Infinity (N : Num) return Boolean is (Is_NaN (N - N));
   --  Return True if N is an infinity. Used to avoid propagating meaningless
   --  errors when the result of a product is an infinity.

   function Is_Zero (N : Num) return Boolean is (N = -N);
   --  Return True if N is a Zero. Used to preserve the sign when the result of
   --  a product is a zero.

   package Product is
      function Two_Prod (A, B : Num) return Double_T;
      function Two_Sqr (A : Num) return Double_T;
   end Product;
   --  The low-level implementation of multiplicative operations

   package body Product is separate;
   --  This is a separate body because the implementation depends on whether a
   --  Fused Multiply-Add instruction is available on the target.

   -------------------
   -- Quick_Two_Sum --
   -------------------

   function Quick_Two_Sum (A, B : Num) return Double_T is
      S : constant Num := A + B;
      V : constant Num := S - A;
      E : constant Num := B - V;

   begin
      return (S, E);
   end Quick_Two_Sum;

   -------------
   -- Two_Sum --
   -------------

   function Two_Sum (A, B : Num) return Double_T is
      S : constant Num := A + B;
      V : constant Num := S - A;
      E : constant Num := (A - (S - V)) + (B - V);

   begin
      return (S, E);
   end Two_Sum;

   --------------
   -- Two_Diff --
   --------------

   function Two_Diff (A, B : Num) return Double_T is
      S : constant Num := A - B;
      V : constant Num := S - A;
      E : constant Num := (A - (S - V)) - (B + V);

   begin
      return (S, E);
   end Two_Diff;

   --------------
   -- Two_Prod --
   --------------

   function Two_Prod (A, B : Num) return Double_T renames Product.Two_Prod;

   -------------
   -- Two_Sqr --
   -------------

   function Two_Sqr (A : Num) return Double_T renames Product.Two_Sqr;

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

   function "+" (A : Double_T; B : Num) return Double_T is
      S : constant Double_T := Two_Sum (A.Hi, B);

   begin
      return Quick_Two_Sum (S.Hi, S.Lo + A.Lo);
   end "+";

   function "+" (A, B : Double_T) return Double_T is
      S1 : constant Double_T := Two_Sum (A.Hi, B.Hi);
      S2 : constant Double_T := Two_Sum (A.Lo, B.Lo);
      S3 : constant Double_T := Quick_Two_Sum (S1.Hi, S1.Lo + S2.Hi);

   begin
      return Quick_Two_Sum (S3.Hi, S3.Lo + S2.Lo);
   end "+";

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

   function "-" (A : Double_T; B : Num) return Double_T is
      D : constant Double_T := Two_Diff (A.Hi, B);

   begin
      return Quick_Two_Sum (D.Hi, D.Lo + A.Lo);
   end "-";

   function "-" (A, B : Double_T) return Double_T is
      D1 : constant Double_T := Two_Diff (A.Hi, B.Hi);
      D2 : constant Double_T := Two_Diff (A.Lo, B.Lo);
      D3 : constant Double_T := Quick_Two_Sum (D1.Hi, D1.Lo + D2.Hi);

   begin
      return Quick_Two_Sum (D3.Hi, D3.Lo + D2.Lo);
   end "-";

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

   function "*" (A : Double_T; B : Num) return Double_T is
      P : constant Double_T := Two_Prod (A.Hi, B);

   begin
      if Is_Infinity (P.Hi) or else Is_Zero (P.Hi) then
         return (P.Hi, 0.0);
      else
         return Quick_Two_Sum (P.Hi, P.Lo + A.Lo * B);
      end if;
   end "*";

   function "*" (A, B : Double_T) return Double_T is
      P : constant Double_T := Two_Prod (A.Hi, B.Hi);

   begin
      if Is_Infinity (P.Hi) or else Is_Zero (P.Hi) then
         return (P.Hi, 0.0);
      else
         return Quick_Two_Sum (P.Hi, P.Lo + A.Hi * B.Lo + A.Lo * B.Hi);
      end if;
   end "*";

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

   function "/" (A : Double_T; B : Num) return Double_T is
      Q1, Q2 : Num;
      P, R   : Double_T;

   begin
      Q1 := A.Hi / B;

      --  Compute R = A - B * Q1

      P := Two_Prod (B, Q1);
      R := Two_Diff (A.Hi, P.Hi);
      R.Lo := (R.Lo + A.Lo) - P.Lo;

      Q2 := (R.Hi + R.Lo) / B;

      return Quick_Two_Sum (Q1, Q2);
   end "/";

   function "/" (A, B : Double_T) return Double_T is
      Q1, Q2, Q3 : Num;
      R, S       : Double_T;

   begin
      Q1 := A.Hi / B.Hi;
      R := A - B * Q1;

      Q2 := R.Hi / B.Hi;
      R := R - B * Q2;

      Q3 := R.Hi / B.Hi;

      S := Quick_Two_Sum (Q1, Q2);
      return Quick_Two_Sum (S.Hi, S.Lo + Q3);
   end "/";

   ---------
   -- Sqr --
   ---------

   function Sqr (A : Double_T) return Double_T is
      Q : constant Double_T := Two_Sqr (A.Hi);

   begin
      if Is_Infinity (Q.Hi) or else Is_Zero (Q.Hi) then
         return (Q.Hi, 0.0);
      else
         return Quick_Two_Sum (Q.Hi, Q.Lo + 2.0 * A.Hi * A.Lo + A.Lo * A.Lo);
      end if;
   end Sqr;

   -------------------
   -- From_Unsigned --
   -------------------

   function From_Unsigned (U : Uns) return Double_T is
   begin
      return To_Double (Num (U));
   end From_Unsigned;

   -----------------
   -- To_Unsigned --
   -----------------

   function To_Unsigned (D : Double_T) return Uns is
      Hi : constant Num := Num'Truncation (D.Hi);

   begin
      --  If the high part is already an integer, add Floor of the low part,
      --  which means subtract Ceiling of its opposite if it is negative.

      if Hi = D.Hi then
         if D.Lo < 0.0 then
            return Uns (Hi) - Uns (Num'Ceiling (-D.Lo));
         else
            return Uns (Hi) + Uns (Num'Floor (D.Lo));
         end if;

      else
         return Uns (Hi);
      end if;
   end To_Unsigned;

end System.Double_Real;