From: Pierre-Marie de Rodat <derodat@adacore.com>
To: gcc-patches@gcc.gnu.org
Cc: Eric Botcazou <ebotcazou@adacore.com>
Subject: [Ada] Eliminate early roundoff error for Long_Long_Float on x86
Date: Wed, 28 Apr 2021 05:41:39 -0400 [thread overview]
Message-ID: <20210428094139.GA139939@adacore.com> (raw)
[-- Attachment #1: Type: text/plain, Size: 482 bytes --]
This overcomes the lack of fused multiply-add instruction on the x87
FPU by doing an iterated addition with exact error handling for the
last digit taken into account for the mantissa.
Tested on x86_64-pc-linux-gnu, committed on trunk
gcc/ada/
* libgnat/s-valrea.adb (Fast2Sum): New function.
(Integer_to_Real): Use it in an iterated addition with exact
error handling for the case where an extra digit is needed.
Move local variable now only used in the exponentiation case.
[-- Attachment #2: patch.diff --]
[-- Type: text/x-diff, Size: 4827 bytes --]
diff --git a/gcc/ada/libgnat/s-valrea.adb b/gcc/ada/libgnat/s-valrea.adb
--- a/gcc/ada/libgnat/s-valrea.adb
+++ b/gcc/ada/libgnat/s-valrea.adb
@@ -46,7 +46,7 @@ package body System.Val_Real is
-- If the mantissa of the floating-point type is almost as large as the
-- unsigned type, we do not have enough space for an extra digit in the
-- unsigned type so we handle the extra digit separately, at the cost of
- -- a potential roundoff error.
+ -- a bit more work in Integer_to_Real.
Precision_Limit : constant Uns :=
(if Need_Extra then 2**Num'Machine_Mantissa - 1 else 2**Uns'Size - 1);
@@ -76,6 +76,10 @@ package body System.Val_Real is
7 => 5836, 8 => 5461, 9 => 5168, 10 => 4932, 11 => 4736,
12 => 4570, 13 => 4427, 14 => 4303, 15 => 4193, 16 => 4095);
+ function Fast2Sum (A, B : Num; Err : in out Num) return Num;
+ -- This is the classical Fast2Sum function assuming round to nearest,
+ -- with the error accumulated into Err.
+
function Integer_to_Real
(Str : String;
Val : Uns;
@@ -85,6 +89,25 @@ package body System.Val_Real is
Minus : Boolean) return Num;
-- Convert the real value from integer to real representation
+ --------------
+ -- Fast2Sum --
+ --------------
+
+ function Fast2Sum (A, B : Num; Err : in out Num) return Num is
+ S, Z : Num;
+
+ begin
+ pragma Assert (abs (A) >= abs (B));
+
+ S := A + B;
+ Z := S - A;
+ Z := B - Z;
+
+ Err := Err + Z;
+
+ return S;
+ end Fast2Sum;
+
---------------------
-- Integer_to_Real --
---------------------
@@ -110,8 +133,6 @@ package body System.Val_Real is
else raise Program_Error);
-- Maximum exponent of the base that can fit in Num
- B : constant Num := Num (Base);
-
R_Val : Num;
S : Integer := Scale;
@@ -129,12 +150,53 @@ package body System.Val_Real is
R_Val := Num (Val);
- -- Take into account the extra digit, if need be. In this case, the
- -- three operands are exact, so using an FMA would be ideal.
+ -- Take into account the extra digit, i.e. do the two computations
+
+ -- (1) R_Val := R_Val * Num (B) + Num (Extra)
+ -- (2) S := S - 1
+
+ -- In the first, the three operands are exact, so using an FMA would
+ -- be ideal, but we are most likely running on the x87 FPU, hence we
+ -- may not have one. That is why we turn the multiplication into an
+ -- iterated addition with exact error handling, so that we can do a
+ -- single rounding at the end.
if Need_Extra and then Extra > 0 then
- R_Val := R_Val * B + Num (Extra);
- S := S - 1;
+ declare
+ B : Unsigned := Base;
+
+ Acc : Num := 0.0;
+ Err : Num := 0.0;
+ Fac : Num := R_Val;
+
+ begin
+ loop
+ -- If B is odd, add one factor. Note that the accumulator is
+ -- never larger than the factor at this point (it is in fact
+ -- never larger than the factor minus the initial value).
+
+ if B rem 2 /= 0 then
+ Acc := (if Acc = 0.0 then Fac else Fast2Sum (Fac, Acc, Err));
+ exit when B = 1;
+ end if;
+
+ -- Now B is (morally) even, halve it and double the factor,
+ -- which is always an exact operation.
+
+ B := B / 2;
+ Fac := Fac * 2.0;
+ end loop;
+
+ -- Add Extra to the error, which are both small integers
+
+ Err := Err + Num (Extra);
+
+ -- Acc + Err is the exact result before rounding
+
+ R_Val := Acc + Err;
+
+ S := S - 1;
+ end;
end if;
-- Compute the final value
@@ -207,17 +269,23 @@ package body System.Val_Real is
-- an artificial underflow.
when others =>
- if S > 0 then
- R_Val := R_Val * B ** S;
+ declare
+ B : constant Num := Num (Base);
- else
- if S < -Maxexp then
- R_Val := R_Val / B ** Maxexp;
- S := S + Maxexp;
- end if;
+ begin
- R_Val := R_Val / B ** (-S);
- end if;
+ if S > 0 then
+ R_Val := R_Val * B ** S;
+
+ else
+ if S < -Maxexp then
+ R_Val := R_Val / B ** Maxexp;
+ S := S + Maxexp;
+ end if;
+
+ R_Val := R_Val / B ** (-S);
+ end if;
+ end;
end case;
end if;
reply other threads:[~2021-04-28 9:41 UTC|newest]
Thread overview: [no followups] expand[flat|nested] mbox.gz Atom feed
Reply instructions:
You may reply publicly to this message via plain-text email
using any one of the following methods:
* Save the following mbox file, import it into your mail client,
and reply-to-all from there: mbox
Avoid top-posting and favor interleaved quoting:
https://en.wikipedia.org/wiki/Posting_style#Interleaved_style
* Reply using the --to, --cc, and --in-reply-to
switches of git-send-email(1):
git send-email \
--in-reply-to=20210428094139.GA139939@adacore.com \
--to=derodat@adacore.com \
--cc=ebotcazou@adacore.com \
--cc=gcc-patches@gcc.gnu.org \
/path/to/YOUR_REPLY
https://kernel.org/pub/software/scm/git/docs/git-send-email.html
* If your mail client supports setting the In-Reply-To header
via mailto: links, try the mailto: link
Be sure your reply has a Subject: header at the top and a blank line
before the message body.
This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox;
as well as URLs for read-only IMAP folder(s) and NNTP newsgroup(s).