From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Received: from rock.gnat.com (rock.gnat.com [205.232.38.15]) by sourceware.org (Postfix) with ESMTP id 455D3399A419 for ; Wed, 28 Apr 2021 09:41:42 +0000 (GMT) DMARC-Filter: OpenDMARC Filter v1.3.2 sourceware.org 455D3399A419 Authentication-Results: sourceware.org; dmarc=none (p=none dis=none) header.from=adacore.com Authentication-Results: sourceware.org; spf=pass smtp.mailfrom=derodat@adacore.com Received: from localhost (localhost.localdomain [127.0.0.1]) by filtered-rock.gnat.com (Postfix) with ESMTP id 0D73F56149; Wed, 28 Apr 2021 05:41:40 -0400 (EDT) X-Virus-Scanned: Debian amavisd-new at gnat.com Received: from rock.gnat.com ([127.0.0.1]) by localhost (rock.gnat.com [127.0.0.1]) (amavisd-new, port 10024) with LMTP id Rz91Qa1NvHtv; Wed, 28 Apr 2021 05:41:39 -0400 (EDT) Received: from tron.gnat.com (tron.gnat.com [205.232.38.10]) by rock.gnat.com (Postfix) with ESMTP id D1B2F56140; Wed, 28 Apr 2021 05:41:39 -0400 (EDT) Received: by tron.gnat.com (Postfix, from userid 4862) id D0DC312D; Wed, 28 Apr 2021 05:41:39 -0400 (EDT) Date: Wed, 28 Apr 2021 05:41:39 -0400 From: Pierre-Marie de Rodat To: gcc-patches@gcc.gnu.org Cc: Eric Botcazou Subject: [Ada] Eliminate early roundoff error for Long_Long_Float on x86 Message-ID: <20210428094139.GA139939@adacore.com> MIME-Version: 1.0 Content-Type: multipart/mixed; boundary="ZGiS0Q5IWpPtfppv" Content-Disposition: inline User-Agent: Mutt/1.5.23 (2014-03-12) X-Spam-Status: No, score=-12.1 required=5.0 tests=BAYES_00, GIT_PATCH_0, KAM_DMARC_STATUS, KAM_NUMSUBJECT, SPF_HELO_NONE, SPF_PASS, TXREP autolearn=ham autolearn_force=no version=3.4.2 X-Spam-Checker-Version: SpamAssassin 3.4.2 (2018-09-13) on server2.sourceware.org X-BeenThere: gcc-patches@gcc.gnu.org X-Mailman-Version: 2.1.29 Precedence: list List-Id: Gcc-patches mailing list List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Wed, 28 Apr 2021 09:41:45 -0000 --ZGiS0Q5IWpPtfppv Content-Type: text/plain; charset=us-ascii Content-Disposition: inline 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. --ZGiS0Q5IWpPtfppv Content-Type: text/x-diff; charset=us-ascii Content-Disposition: attachment; filename="patch.diff" 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; --ZGiS0Q5IWpPtfppv--