diff options
Diffstat (limited to 'gcc/ada/libgnat/s-dourea.adb')
| -rw-r--r-- | gcc/ada/libgnat/s-dourea.adb | 258 |
1 files changed, 258 insertions, 0 deletions
diff --git a/gcc/ada/libgnat/s-dourea.adb b/gcc/ada/libgnat/s-dourea.adb new file mode 100644 index 0000000..53bed1d --- /dev/null +++ b/gcc/ada/libgnat/s-dourea.adb @@ -0,0 +1,258 @@ +------------------------------------------------------------------------------ +-- -- +-- 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; |