From mboxrd@z Thu Jan 1 00:00:00 1970 X-Spam-Checker-Version: SpamAssassin 3.4.4 (2020-01-24) on inbox.vuxu.org X-Spam-Level: X-Spam-Status: No, score=-3.3 required=5.0 tests=MAILING_LIST_MULTI, RCVD_IN_DNSWL_MED,RCVD_IN_MSPIKE_H3,RCVD_IN_MSPIKE_WL, T_SCC_BODY_TEXT_LINE autolearn=ham autolearn_force=no version=3.4.4 Received: (qmail 31443 invoked from network); 11 Nov 2023 10:27:16 -0000 Received: from second.openwall.net (193.110.157.125) by inbox.vuxu.org with ESMTPUTF8; 11 Nov 2023 10:27:16 -0000 Received: (qmail 7177 invoked by uid 550); 11 Nov 2023 10:27:12 -0000 Mailing-List: contact musl-help@lists.openwall.com; run by ezmlm Precedence: bulk List-Post: List-Help: List-Unsubscribe: List-Subscribe: List-ID: Reply-To: musl@lists.openwall.com Received: (qmail 6115 invoked from network); 11 Nov 2023 10:27:11 -0000 Date: Sat, 11 Nov 2023 11:26:58 +0100 From: Szabolcs Nagy To: Damian McGuckin Cc: musl@lists.openwall.com Message-ID: <20231111102658.GH1427497@port70.net> Mail-Followup-To: Damian McGuckin , musl@lists.openwall.com References: <239849df-fbf3-95e5-f1bd-f1863d43fe8@esi.com.au> MIME-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline In-Reply-To: <239849df-fbf3-95e5-f1bd-f1863d43fe8@esi.com.au> Subject: Re: [musl] Double Rounding in .... * Damian McGuckin [2023-11-11 19:11:41 +1100]: > > hypot, hypotf, hypotl > and > scalbn, scalbnf, scalbnl > > Does anybody have a good/clear explanation for how double rounding occurs in > the case of squaring-of/multiplication-by small numbers and how the action > taken in the code addresses that? Not sure if it needs to mention n * 2^n is normally exact, but in the subnormal range there are less precision bits so rounding can occur. if you do the scaling in two steps (e.g. because the scale factor has such a big exponent that it is not representible in a single float) then you can get double rounding. e.g. scalbn(0x1.000000000000bp-1, -1024) == 0x1.0000000000008p-1025 normally the scaling is done in two steps: x*0x1p-1022*0x1p-2 because up to 0x1p-1022 the scaling is easy to construct from an int, while 0x1p-1024 is subnormal. with naive implementation the bit pattern at the end is ..001011 x ..00110 x*0x1p-1022 (rounded) ..010 x*0x1p-2 (rounded) with the fix ..001011 x*0x1p-969 (969 = 1022-53) ..001 x*0x1p-55 (rounded) so simply split the scale factor such that the first scaling can only do rounding if the second scaling (<0x1p-53) completely rounds the entire result away. > > FLT_EVAL_METHOD > > because certainly the float and double versions of those routines are > written in terms of float_t and double_t respectively. > > Thanks - Damian > > Pacific Engineering Systems International ..... 20D Grose St, Glebe NSW 2037 > Ph:+61-2-8571-0847 .. Fx:+61-2-9692-9623 | unsolicited email not wanted here > Views & opinions here are mine and not those of any past or present employer