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