From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Received: from EUR04-DB3-obe.outbound.protection.outlook.com (mail-db3eur04on2055.outbound.protection.outlook.com [40.107.6.55]) by sourceware.org (Postfix) with ESMTPS id C66813858D37 for ; Mon, 3 Apr 2023 13:33:28 +0000 (GMT) DMARC-Filter: OpenDMARC Filter v1.4.2 sourceware.org C66813858D37 Authentication-Results: sourceware.org; dmarc=pass (p=none dis=none) header.from=arm.com Authentication-Results: sourceware.org; spf=pass smtp.mailfrom=arm.com DKIM-Signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=armh.onmicrosoft.com; s=selector2-armh-onmicrosoft-com; h=From:Date:Subject:Message-ID:Content-Type:MIME-Version:X-MS-Exchange-SenderADCheck; bh=icy0J/wzGNPIw1fmJIeq36y//OmFsUqdLhHaGBjDOL8=; b=vl8g1wCo4fo7yM890cZ+05MOGDw9+OFMaS6GxtAvGXsf5bzTuwr7U8pQSrpiADLR5G6xxKrS2HpFoxQYBY0X61U9xcXP23IEcMrFXMhJXnAHwFNG1az9gtkezuThqnZ680nFSZ3CBvPLmF13p0WW6Vr8NiHMx5vLKkfOLf4wPIU= Received: from AM6PR01CA0072.eurprd01.prod.exchangelabs.com (2603:10a6:20b:e0::49) by DU0PR08MB9485.eurprd08.prod.outlook.com (2603:10a6:10:42e::20) with Microsoft SMTP Server (version=TLS1_2, cipher=TLS_ECDHE_RSA_WITH_AES_256_GCM_SHA384) id 15.20.6254.30; Mon, 3 Apr 2023 13:33:25 +0000 Received: from AM7EUR03FT043.eop-EUR03.prod.protection.outlook.com (2603:10a6:20b:e0:cafe::85) by AM6PR01CA0072.outlook.office365.com (2603:10a6:20b:e0::49) with Microsoft SMTP Server (version=TLS1_2, cipher=TLS_ECDHE_RSA_WITH_AES_256_GCM_SHA384) id 15.20.6254.22 via Frontend Transport; Mon, 3 Apr 2023 13:33:25 +0000 X-MS-Exchange-Authentication-Results: spf=pass (sender IP is 63.35.35.123) smtp.mailfrom=arm.com; dkim=pass (signature was verified) header.d=armh.onmicrosoft.com;dmarc=pass action=none header.from=arm.com; Received-SPF: Pass (protection.outlook.com: domain of arm.com designates 63.35.35.123 as permitted sender) receiver=protection.outlook.com; client-ip=63.35.35.123; helo=64aa7808-outbound-1.mta.getcheckrecipient.com; pr=C Received: from 64aa7808-outbound-1.mta.getcheckrecipient.com (63.35.35.123) by AM7EUR03FT043.mail.protection.outlook.com (100.127.140.160) with Microsoft SMTP Server (version=TLS1_2, cipher=TLS_ECDHE_RSA_WITH_AES_256_GCM_SHA384) id 15.20.6277.26 via Frontend Transport; Mon, 3 Apr 2023 13:33:25 +0000 Received: ("Tessian outbound 945aec65ec65:v136"); Mon, 03 Apr 2023 13:33:24 +0000 X-CheckRecipientChecked: true X-CR-MTA-CID: ceaae2b092588c56 X-CR-MTA-TID: 64aa7808 Received: from 551cc47995aa.1 by 64aa7808-outbound-1.mta.getcheckrecipient.com id 7989B3DA-9710-41A5-BEBA-84969AEC9629.1; Mon, 03 Apr 2023 13:33:19 +0000 Received: from EUR04-DB3-obe.outbound.protection.outlook.com by 64aa7808-outbound-1.mta.getcheckrecipient.com with ESMTPS id 551cc47995aa.1 (version=TLSv1.2 cipher=ECDHE-RSA-AES256-GCM-SHA384); Mon, 03 Apr 2023 13:33:19 +0000 ARC-Seal: i=1; a=rsa-sha256; s=arcselector9901; d=microsoft.com; cv=none; b=NlpdOv8D3FIZj1PIGN2m1XTmzZpJwqOBmOdtskn9qrPH6nSq2KquTiGjSsern9gUWssXF/l54vBZG+Py+90FfWhWeWL5qTMVGS54NnCUrTDRCrz38CNOoZTjNJU43dgEUwn6p7lPfqNvG7wWqBKyZz6xstbwQwRHxH7ouSiBlMu9cqWjeFRR+pMd8VW+CLheNi5xxE3SvYb2m9xIeODpQZUlg331b9Ta5Ou4a3JyDJxOsQAltiVx2rrrYai/6Z9eaiMGH2y11jmfdPoyhgHg4bjDaZAWNzhJRveO4oOM+FtW7ZsGhtUzc7+xwagtQqtBT3Ok2HyJ77sWkpTP8PV+YA== ARC-Message-Signature: i=1; a=rsa-sha256; c=relaxed/relaxed; d=microsoft.com; s=arcselector9901; h=From:Date:Subject:Message-ID:Content-Type:MIME-Version:X-MS-Exchange-AntiSpam-MessageData-ChunkCount:X-MS-Exchange-AntiSpam-MessageData-0:X-MS-Exchange-AntiSpam-MessageData-1; bh=icy0J/wzGNPIw1fmJIeq36y//OmFsUqdLhHaGBjDOL8=; b=EljaANdJ+gu9v3vwaak8WxEybNnwiSjZcUxodBcUfxvEB6g+idvC05hZTWGVbLPh/YRmxmZNN0LeNXWq+EISIbItwKLSEj9jWyp9IIjMwiOuVgYCBtMhUoMw8X9yqhJmeYa7o1hPj34aXkEre3ziOCGg17Lg189N/FTR9noJAhhbjp31l++qZls7CAsM01M8ryORSUTXyNnFeLfiOlyBCJO3co50ppM2ukY/oaFPXrE8/xczec1um4GG1HwmGvHwXhwMtpafjNTCLhvvxSmSWYbL24DAzqfO9qRqKNXaA6u39sBAS3s5TSfPLtIg3pSZBKQ+TupRYJnUVqs3hvqFjg== ARC-Authentication-Results: i=1; mx.microsoft.com 1; spf=pass smtp.mailfrom=arm.com; dmarc=pass action=none header.from=arm.com; dkim=pass header.d=arm.com; arc=none DKIM-Signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=armh.onmicrosoft.com; s=selector2-armh-onmicrosoft-com; h=From:Date:Subject:Message-ID:Content-Type:MIME-Version:X-MS-Exchange-SenderADCheck; bh=icy0J/wzGNPIw1fmJIeq36y//OmFsUqdLhHaGBjDOL8=; b=vl8g1wCo4fo7yM890cZ+05MOGDw9+OFMaS6GxtAvGXsf5bzTuwr7U8pQSrpiADLR5G6xxKrS2HpFoxQYBY0X61U9xcXP23IEcMrFXMhJXnAHwFNG1az9gtkezuThqnZ680nFSZ3CBvPLmF13p0WW6Vr8NiHMx5vLKkfOLf4wPIU= Received: from PAWPR08MB8982.eurprd08.prod.outlook.com (2603:10a6:102:33f::20) by DU0PR08MB7787.eurprd08.prod.outlook.com (2603:10a6:10:3b8::14) with Microsoft SMTP Server (version=TLS1_2, cipher=TLS_ECDHE_RSA_WITH_AES_256_GCM_SHA384) id 15.20.6254.33; Mon, 3 Apr 2023 13:33:16 +0000 Received: from PAWPR08MB8982.eurprd08.prod.outlook.com ([fe80::13be:967d:6e80:432f]) by PAWPR08MB8982.eurprd08.prod.outlook.com ([fe80::13be:967d:6e80:432f%9]) with mapi id 15.20.6254.030; Mon, 3 Apr 2023 13:33:16 +0000 From: Wilco Dijkstra To: Adhemerval Zanella , "libc-alpha@sourceware.org" , "H . J . Lu" CC: kirill Subject: Re: [PATCH v4 4/5] math: Improve fmodf Thread-Topic: [PATCH v4 4/5] math: Improve fmodf Thread-Index: AQHZW0VEZalzvXDsVEy4T2C03hdFPa8ZqqQi Date: Mon, 3 Apr 2023 13:33:16 +0000 Message-ID: References: <20230320160118.352206-1-adhemerval.zanella@linaro.org> <20230320160118.352206-5-adhemerval.zanella@linaro.org> In-Reply-To: <20230320160118.352206-5-adhemerval.zanella@linaro.org> Accept-Language: en-GB, en-US Content-Language: en-GB X-MS-Has-Attach: X-MS-TNEF-Correlator: msip_labels: Authentication-Results-Original: dkim=none (message not signed) header.d=none;dmarc=none action=none header.from=arm.com; x-ms-traffictypediagnostic: PAWPR08MB8982:EE_|DU0PR08MB7787:EE_|AM7EUR03FT043:EE_|DU0PR08MB9485:EE_ X-MS-Office365-Filtering-Correlation-Id: cd9bc2dd-84a6-447f-b685-08db34480064 x-checkrecipientrouted: true nodisclaimer: true X-MS-Exchange-SenderADCheck: 1 X-MS-Exchange-AntiSpam-Relay: 0 X-Microsoft-Antispam-Untrusted: BCL:0; X-Microsoft-Antispam-Message-Info-Original: qWqg+RkeOWOLnoD/UKvjtM36LrcUzg9WcDA3iKuUw8Wau/N3iVFeEcJ+yMvebSUTPPY/D6PdtAbUSPiiZ5eKW+xtPkCFkguoLxorFX+SRRBqBIGHfG4xn+6b99GgrX7JVeHBtmkxSj99pBlB8M71fzCWP5ykKASI+Os3Gh3TAYFBJiV0HRzVTMRVXGJrMX5RqEeKLTYHoFwwmIv5f7fOGqQP47m32n5DctuDiYgz8NEjWB+sFbOLSleQxnG6neipLWNkC6UnCwb7kVA32YqmkU57D7kHXMmuPtJYMwth6S5KKkCKTIrIhQMjN+BdjsxQXOeq1fHN0WtgWXeh9PmTqcJD+ojyNKdPufZYvJ7rJSSh/l3kN8C69HUjnoQN74y9VqmhoGO91p7UM2Y40lRA9AUrA0ppvDLTIjO5fqTvzE6pLut6Smul9jF4uC8ACKuKMNpmJV2vihQOWjZ8JAYZebsnwmFXZ3exTQgi0RT3aGgvYICMW9N1UKa3GcEOdJgr+vOclK/HQRFwtXzdWWWOMha83qL4B+tbhd4TMRTMwHC0+ETI0rv+t67Y0McUE5rlghvnQSjtrXQAb+Lbys5+VYp/epezw07SyS7WATjhne17ZVuF6Iok0VHof60S+wzR X-Forefront-Antispam-Report-Untrusted: CIP:255.255.255.255;CTRY:;LANG:en;SCL:1;SRV:;IPV:NLI;SFV:NSPM;H:PAWPR08MB8982.eurprd08.prod.outlook.com;PTR:;CAT:NONE;SFS:(13230028)(4636009)(136003)(346002)(39860400002)(396003)(376002)(366004)(451199021)(55016003)(186003)(26005)(2906002)(5660300002)(86362001)(38070700005)(9686003)(8936002)(66556008)(66946007)(478600001)(76116006)(66476007)(91956017)(316002)(110136005)(4326008)(52536014)(83380400001)(6506007)(41300700001)(38100700002)(33656002)(122000001)(71200400001)(7696005)(64756008)(8676002)(66446008)(2004002);DIR:OUT;SFP:1101; Content-Type: text/plain; charset="iso-8859-1" Content-Transfer-Encoding: quoted-printable MIME-Version: 1.0 X-MS-Exchange-Transport-CrossTenantHeadersStamped: DU0PR08MB7787 Original-Authentication-Results: dkim=none (message not signed) header.d=none;dmarc=none action=none header.from=arm.com; X-EOPAttributedMessage: 0 X-MS-Exchange-Transport-CrossTenantHeadersStripped: AM7EUR03FT043.eop-EUR03.prod.protection.outlook.com X-MS-PublicTrafficType: Email X-MS-Office365-Filtering-Correlation-Id-Prvs: 2b77e337-5b4a-4161-c3dd-08db3447fb06 X-Microsoft-Antispam: BCL:0; X-Microsoft-Antispam-Message-Info: UOfB5Qsy5wwSryjA73cLr30B83bSEokI5PzKToZ/eFabzH1ewI8YjRPNcjutJkwP/3cb9SkTDnsBbXH/lNlPVTAjEUD0SUaUV20OwRaYTY/CWljgjN1HxMwEezYTWZG6NFxmNljUfvvRt6VhPTTKJSuVEAibvKMde8nhcI0cVv7a6sokdcJFLlHuIdWPtR/w5L0eOdDrV4o42CynSaQE/xlHGp0g99RvZgwgZeDpNV/CR6QKCN2yZ9A3zH1jrbZpLSfku7spSbGzCwuAhhTNKBoDXJTD9aFIHMyaysQndPmAeh2E4eYL67vcwTfXxO7gY0nEAwHZZ87EZR+q+OE+chFjv1pixJ8GT0ceHsYT+YcMzbmG6xlnnqSV7QZ71AOkKlP2MuXuLVdh24UlXFwIk3dNRyz9SYk4yPYM6n8UiD60XdAU7Yjaox16pRyXWlnhf7+3PsOCaR6mzOR8Je9LGAEmjBxTV415ys7MdTDCnS4wpFLXnOnYLdOHXVatrR7gDZ5kzN+8s6cZkqt9b6boqaUUDMAWggQ4EZI3x7WK/EP3Ehys8UhgSHtT02hGVcj7qcuiUz0nWUMn0bdpBXrB+hk/637m4ruJ7kXnEKeyB53vIVQHZpYGdz2Z+iI4tBr3tj5uAZrx5QkIwabibU4KXDPYajjjr+0lw27FtRY1w5aQzNIBsaPWXXX5s1Meh55Abh581E55vf8YWKF77odP6w== X-Forefront-Antispam-Report: CIP:63.35.35.123;CTRY:IE;LANG:en;SCL:1;SRV:;IPV:CAL;SFV:NSPM;H:64aa7808-outbound-1.mta.getcheckrecipient.com;PTR:ec2-63-35-35-123.eu-west-1.compute.amazonaws.com;CAT:NONE;SFS:(13230028)(4636009)(376002)(39860400002)(346002)(396003)(136003)(451199021)(46966006)(36840700001)(40470700004)(7696005)(55016003)(40480700001)(2906002)(33656002)(41300700001)(4326008)(40460700003)(6506007)(107886003)(26005)(186003)(70586007)(70206006)(478600001)(5660300002)(52536014)(316002)(110136005)(8936002)(8676002)(82740400003)(36860700001)(356005)(83380400001)(47076005)(336012)(9686003)(81166007)(82310400005)(86362001)(2004002);DIR:OUT;SFP:1101; X-OriginatorOrg: arm.com X-MS-Exchange-CrossTenant-OriginalArrivalTime: 03 Apr 2023 13:33:25.1165 (UTC) X-MS-Exchange-CrossTenant-Network-Message-Id: cd9bc2dd-84a6-447f-b685-08db34480064 X-MS-Exchange-CrossTenant-Id: f34e5979-57d9-4aaa-ad4d-b122a662184d X-MS-Exchange-CrossTenant-OriginalAttributedTenantConnectingIp: TenantId=f34e5979-57d9-4aaa-ad4d-b122a662184d;Ip=[63.35.35.123];Helo=[64aa7808-outbound-1.mta.getcheckrecipient.com] X-MS-Exchange-CrossTenant-AuthSource: AM7EUR03FT043.eop-EUR03.prod.protection.outlook.com X-MS-Exchange-CrossTenant-AuthAs: Anonymous X-MS-Exchange-CrossTenant-FromEntityHeader: HybridOnPrem X-MS-Exchange-Transport-CrossTenantHeadersStamped: DU0PR08MB9485 X-Spam-Status: No, score=-9.1 required=5.0 tests=BAYES_00,BODY_8BITS,DKIM_SIGNED,DKIM_VALID,FORGED_SPF_HELO,GIT_PATCH_0,KAM_ASCII_DIVIDERS,KAM_DMARC_NONE,KAM_SHORT,RCVD_IN_DNSWL_NONE,RCVD_IN_MSPIKE_H2,SPF_HELO_PASS,SPF_NONE,TXREP,UNPARSEABLE_RELAY autolearn=ham autolearn_force=no version=3.4.6 X-Spam-Checker-Version: SpamAssassin 3.4.6 (2021-04-09) on server2.sourceware.org List-Id: Hi Adhemerval,=0A= =0A= LGTM=0A= =0A= Reviewed-by: Wilco Dijkstra =0A= =0A= Cheers,=0A= Wilco=0A= =0A= diff --git a/math/libm-test-fmod.inc b/math/libm-test-fmod.inc=0A= index 8841c1332c..43376c3df2 100644=0A= --- a/math/libm-test-fmod.inc=0A= +++ b/math/libm-test-fmod.inc=0A= @@ -213,6 +213,10 @@ static const struct test_ff_f_data fmod_test_data[] = =3D=0A= =A0=A0=A0=A0 TEST_ff_f (fmod, -0x1p127L, -0x3p-148L, -0x1p-147L, NO_INEXACT= _EXCEPTION|ERRNO_UNCHANGED),=0A= =A0=A0=A0=A0 TEST_ff_f (fmod, -0x1p127L, 0x3p-126L, -0x1p-125L, NO_INEXACT_= EXCEPTION|ERRNO_UNCHANGED),=0A= =A0=A0=A0=A0 TEST_ff_f (fmod, -0x1p127L, -0x3p-126L, -0x1p-125L, NO_INEXACT= _EXCEPTION|ERRNO_UNCHANGED),=0A= +=A0=A0=A0 TEST_ff_f (fmod, 0x1.3a3e6p-127, 0x1.8b8338p-128, 0x1.d1f31p-129= , NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),=0A= +=A0=A0=A0 TEST_ff_f (fmod, 0x1.3a3e6p-127, -0x1.8b8338p-128, 0x1.d1f31p-12= 9, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),=0A= +=A0=A0=A0 TEST_ff_f (fmod, -0x1.3a3e6p-127, 0x1.8b8338p-128, -0x1.d1f31p-1= 29, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),=0A= +=A0=A0=A0 TEST_ff_f (fmod, -0x1.3a3e6p-127, -0x1.8b8338p-128, -0x1.d1f31p-= 129, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),=0A= =A0#if !TEST_COND_binary32=0A= =A0=A0=A0=A0 TEST_ff_f (fmod, 0x1p1023L, 0x3p-1074L, 0x1p-1073L, NO_INEXACT= _EXCEPTION|ERRNO_UNCHANGED),=0A= =A0=A0=A0=A0 TEST_ff_f (fmod, 0x1p1023L, -0x3p-1074L, 0x1p-1073L, NO_INEXAC= T_EXCEPTION|ERRNO_UNCHANGED),=0A= =0A= OK=0A= =0A= diff --git a/sysdeps/ieee754/flt-32/e_fmodf.c b/sysdeps/ieee754/flt-32/e_fm= odf.c=0A= index b71c4f754f..53db32cbe2 100644=0A= --- a/sysdeps/ieee754/flt-32/e_fmodf.c=0A= +++ b/sysdeps/ieee754/flt-32/e_fmodf.c=0A= @@ -1,102 +1,144 @@=0A= -/* e_fmodf.c -- float version of e_fmod.c.=0A= - */=0A= -=0A= -/*=0A= - * =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=0A= - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.=0A= - *=0A= - * Developed at SunPro, a Sun Microsystems, Inc. business.=0A= - * Permission to use, copy, modify, and distribute this=0A= - * software is freely granted, provided that this notice=0A= - * is preserved.=0A= - * =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=0A= - */=0A= -=0A= -/*=0A= - * __ieee754_fmodf(x,y)=0A= - * Return x mod y in exact arithmetic=0A= - * Method: shift and subtract=0A= - */=0A= +/* Floating-point remainder function.=0A= +=A0=A0 Copyright (C) 2023 Free Software Foundation, Inc.=0A= +=A0=A0 This file is part of the GNU C Library.=0A= +=0A= +=A0=A0 The GNU C Library is free software; you can redistribute it and/or= =0A= +=A0=A0 modify it under the terms of the GNU Lesser General Public=0A= +=A0=A0 License as published by the Free Software Foundation; either=0A= +=A0=A0 version 2.1 of the License, or (at your option) any later version.= =0A= +=0A= +=A0=A0 The GNU C Library is distributed in the hope that it will be useful= ,=0A= +=A0=A0 but WITHOUT ANY WARRANTY; without even the implied warranty of=0A= +=A0=A0 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.=A0 See the GNU= =0A= +=A0=A0 Lesser General Public License for more details.=0A= +=0A= +=A0=A0 You should have received a copy of the GNU Lesser General Public=0A= +=A0=A0 License along with the GNU C Library; if not, see=0A= +=A0=A0 .=A0 */=0A= =A0=0A= -#include =0A= -#include =0A= =A0#include =0A= +#include =0A= +#include "math_config.h"=0A= +=0A= +/* With x =3D mx * 2^ex and y =3D my * 2^ey (mx, my, ex, ey being integers= ), the=0A= +=A0=A0 simplest implementation is:=0A= +=0A= +=A0=A0 mx * 2^ex =3D=3D 2 * mx * 2^(ex - 1)=0A= =A0=0A= -static const float one =3D 1.0, Zero[] =3D {0.0, -0.0,};=0A= +=A0=A0 while (ex > ey)=0A= +=A0=A0=A0=A0 {=0A= +=A0=A0=A0=A0=A0=A0 mx *=3D 2;=0A= +=A0=A0=A0=A0=A0=A0 --ex;=0A= +=A0=A0=A0=A0=A0=A0 mx %=3D my;=0A= +=A0=A0=A0=A0 }=0A= +=0A= +=A0=A0 With mx/my being mantissa of double floating pointer, on each step = the=0A= +=A0=A0 argument reduction can be improved 11 (which is sizeo of uint64_t m= inus=0A= +=A0=A0 MANTISSA_WIDTH plus the signal bit):=0A= +=0A= +=A0=A0 while (ex > ey)=0A= +=A0=A0=A0=A0 {=0A= +=A0=A0=A0=A0=A0=A0 mx << 11;=0A= +=A0=A0=A0=A0=A0=A0 ex -=3D 11;=0A= +=A0=A0=A0=A0=A0=A0 mx %=3D my;=0A= +=A0=A0=A0=A0 }=A0 */=0A= =A0=0A= =A0float=0A= =A0__ieee754_fmodf (float x, float y)=0A= =A0{=0A= -=A0=A0=A0=A0=A0=A0 int32_t n,hx,hy,hz,ix,iy,sx,i;=0A= -=0A= -=A0=A0=A0=A0=A0=A0 GET_FLOAT_WORD(hx,x);=0A= -=A0=A0=A0=A0=A0=A0 GET_FLOAT_WORD(hy,y);=0A= -=A0=A0=A0=A0=A0=A0 sx =3D hx&0x80000000;=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0= =A0 /* sign of x */=0A= -=A0=A0=A0=A0=A0=A0 hx ^=3Dsx;=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0= /* |x| */=0A= -=A0=A0=A0=A0=A0=A0 hy &=3D 0x7fffffff;=A0=A0=A0=A0=A0=A0 /* |y| */=0A= -=0A= -=A0=A0=A0 /* purge off exception values */=0A= -=A0=A0=A0=A0=A0=A0 if(hy=3D=3D0||(hx>=3D0x7f800000)||=A0=A0=A0=A0=A0=A0=A0= =A0=A0=A0=A0 /* y=3D0,or x not finite */=0A= -=A0=A0=A0=A0=A0=A0=A0=A0=A0 (hy>0x7f800000))=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0= =A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0 /* or y is NaN */=0A= -=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0 return (x*y)/(x*y);=0A= -=A0=A0=A0=A0=A0=A0 if(hx>31];=A0=A0=A0=A0= =A0=A0 /* |x|=3D|y| return x*0*/=0A= -=0A= -=A0=A0=A0 /* determine ix =3D ilogb(x) */=0A= -=A0=A0=A0=A0=A0=A0 if(hx<0x00800000) {=A0=A0=A0=A0 /* subnormal x */=0A= -=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0 for (ix =3D -126,i=3D(hx<<8); i>0; i<<=3D1)= ix -=3D1;=0A= -=A0=A0=A0=A0=A0=A0 } else ix =3D (hx>>23)-127;=0A= -=0A= -=A0=A0=A0 /* determine iy =3D ilogb(y) */=0A= -=A0=A0=A0=A0=A0=A0 if(hy<0x00800000) {=A0=A0=A0=A0 /* subnormal y */=0A= -=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0 for (iy =3D -126,i=3D(hy<<8); i>=3D0; i<<= =3D1) iy -=3D1;=0A= -=A0=A0=A0=A0=A0=A0 } else iy =3D (hy>>23)-127;=0A= -=0A= -=A0=A0=A0 /* set up {hx,lx}, {hy,ly} and align y to x */=0A= -=A0=A0=A0=A0=A0=A0 if(ix >=3D -126)=0A= -=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0 hx =3D 0x00800000|(0x007fffff&hx);=0A= -=A0=A0=A0=A0=A0=A0 else {=A0=A0=A0=A0=A0=A0=A0=A0=A0 /* subnormal x, shift= x to normal */=0A= -=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0 n =3D -126-ix;=0A= -=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0 hx =3D hx<=3D -126)=0A= -=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0 hy =3D 0x00800000|(0x007fffff&hy);=0A= -=A0=A0=A0=A0=A0=A0 else {=A0=A0=A0=A0=A0=A0=A0=A0=A0 /* subnormal y, shift= y to normal */=0A= -=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0 n =3D -126-iy;=0A= -=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0 hy =3D hy<>31];=0A= -=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0 hx =3D hz+hz;=0A= -=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0 }=0A= -=A0=A0=A0=A0=A0=A0 }=0A= -=A0=A0=A0=A0=A0=A0 hz=3Dhx-hy;=0A= -=A0=A0=A0=A0=A0=A0 if(hz>=3D0) {hx=3Dhz;}=0A= -=0A= -=A0=A0=A0 /* convert back to floating value and restore the sign */=0A= -=A0=A0=A0=A0=A0=A0 if(hx=3D=3D0)=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0= =A0=A0=A0=A0=A0=A0=A0=A0 /* return sign(x)*0 */=0A= -=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0 return Zero[(uint32_t)sx>>31];=0A= -=A0=A0=A0=A0=A0=A0 while(hx<0x00800000) {=A0=A0=A0=A0=A0=A0=A0=A0=A0 /* no= rmalize x */=0A= -=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0 hx =3D hx+hx;=0A= -=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0 iy -=3D 1;=0A= -=A0=A0=A0=A0=A0=A0 }=0A= -=A0=A0=A0=A0=A0=A0 if(iy>=3D -126) {=A0=A0=A0=A0=A0=A0=A0=A0 /* normalize = output */=0A= -=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0 hx =3D ((hx-0x00800000)|((iy+127)<<23));=0A= -=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0 SET_FLOAT_WORD(x,hx|sx);=0A= -=A0=A0=A0=A0=A0=A0 } else {=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0 /= * subnormal output */=0A= -=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0 n =3D -126 - iy;=0A= -=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0 hx >>=3D n;=0A= -=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0 SET_FLOAT_WORD(x,hx|sx);=0A= -=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0 x *=3D one;=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0= =A0 /* create necessary signal */=0A= -=A0=A0=A0=A0=A0=A0 }=0A= -=A0=A0=A0=A0=A0=A0 return x;=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0 /* = exact output */=0A= +=A0 uint32_t hx =3D asuint (x);=0A= +=A0 uint32_t hy =3D asuint (y);=0A= +=0A= +=A0 uint32_t sx =3D hx & SIGN_MASK;=0A= +=A0 /* Get |x| and |y|.=A0 */=0A= +=A0 hx ^=3D sx;=0A= +=A0 hy &=3D ~SIGN_MASK;=0A= +=0A= +=A0 /* Special cases:=0A= +=A0=A0=A0=A0 - If x or y is a Nan, NaN is returned.=0A= +=A0=A0=A0=A0 - If x is an inifinity, a NaN is returned.=0A= +=A0=A0=A0=A0 - If y is zero, Nan is returned.=0A= +=A0=A0=A0=A0 - If x is +0/-0, and y is not zero, +0/-0 is returned.=A0 */= =0A= +=A0 if (__glibc_unlikely (hy =3D=3D 0=A0=A0=A0=A0=A0=A0=A0 || hx >=3D EXPO= NENT_MASK || hy > EXPONENT_MASK))=0A= +=A0=A0=A0 return (x * y) / (x * y);=0A= +=0A= +=A0 if (__glibc_unlikely (hx <=3D hy))=0A= +=A0=A0=A0 {=0A= +=A0=A0=A0=A0=A0 if (hx < hy)=0A= +=A0=A0=A0=A0=A0=A0 return x;=0A= +=A0=A0=A0=A0=A0 return asfloat (sx);=0A= +=A0=A0=A0 }=0A= +=0A= +=A0 int ex =3D hx >> MANTISSA_WIDTH;=0A= +=A0 int ey =3D hy >> MANTISSA_WIDTH;=0A= +=0A= +=A0 /* Common case where exponents are close: ey >=3D -103 and |x/y| < 2^8= ,=A0 */=0A= +=A0 if (__glibc_likely (ey > MANTISSA_WIDTH && ex - ey <=3D EXPONENT_WIDTH= ))=0A= +=A0=A0=A0 {=0A= +=A0=A0=A0=A0=A0 uint64_t mx =3D (hx & MANTISSA_MASK) | (MANTISSA_MASK + 1)= ;=0A= +=A0=A0=A0=A0=A0 uint64_t my =3D (hy & MANTISSA_MASK) | (MANTISSA_MASK + 1)= ;=0A= +=0A= +=A0=A0=A0=A0=A0 uint32_t d =3D (ex =3D=3D ey) ? (mx - my) : (mx << (ex - e= y)) % my;=0A= +=A0=A0=A0=A0=A0 return make_float (d, ey - 1, sx);=0A= +=A0=A0=A0 }=0A= +=0A= +=A0 /* Special case, both x and y are subnormal.=A0 */=0A= +=A0 if (__glibc_unlikely (ex =3D=3D 0 && ey =3D=3D 0))=0A= +=A0=A0=A0 return asfloat (sx | hx % hy);=0A= =0A= OK=0A= =0A= +=A0 /* Convert |x| and |y| to 'mx + 2^ex' and 'my + 2^ey'.=A0 Assume that = hx is=0A= +=A0=A0=A0=A0 not subnormal by conditions above.=A0 */=0A= +=A0 uint32_t mx =3D get_mantissa (hx) | (MANTISSA_MASK + 1);=0A= +=A0 ex--;=0A= +=0A= +=A0 uint32_t my =3D get_mantissa (hy) | (MANTISSA_MASK + 1);=0A= +=A0 int lead_zeros_my =3D EXPONENT_WIDTH;=0A= +=A0 if (__glibc_likely (ey > 0))=0A= +=A0=A0=A0 ey--;=0A= +=A0 else=0A= +=A0=A0=A0 {=0A= +=A0=A0=A0=A0=A0 my =3D hy;=0A= +=A0=A0=A0=A0=A0 lead_zeros_my =3D __builtin_clz (my);=0A= +=A0=A0=A0 }=0A= =0A= OK=0A= =0A= +=A0 int tail_zeros_my =3D __builtin_ctz (my);=0A= +=A0 int sides_zeroes =3D lead_zeros_my + tail_zeros_my;=0A= +=A0 int exp_diff =3D ex - ey;=0A= +=0A= +=A0 int right_shift =3D exp_diff < tail_zeros_my ? exp_diff : tail_zeros_m= y;=0A= +=A0 my >>=3D right_shift;=0A= +=A0 exp_diff -=3D right_shift;=0A= +=A0 ey +=3D right_shift;=0A= +=0A= +=A0 int left_shift =3D exp_diff < EXPONENT_WIDTH ? exp_diff : EXPONENT_WID= TH;=0A= +=A0 mx <<=3D left_shift;=0A= +=A0 exp_diff -=3D left_shift;=0A= +=0A= +=A0 mx %=3D my;=0A= +=0A= +=A0 if (__glibc_unlikely (mx =3D=3D 0))=0A= +=A0=A0=A0 return asfloat (sx);=0A= =0A= =0A= OK=0A= =0A= +=A0 if (exp_diff =3D=3D 0)=0A= +=A0=A0=A0 return make_float (mx, ey, sx);=0A= +=0A= +=A0 /* Assume modulo/divide operation is slow, so use multiplication with = invert=0A= +=A0=A0=A0=A0 values.=A0 */=0A= +=A0 uint32_t inv_hy =3D UINT32_MAX / my;=0A= +=A0 while (exp_diff > sides_zeroes) {=0A= +=A0=A0=A0 exp_diff -=3D sides_zeroes;=0A= +=A0=A0=A0 uint32_t hd =3D (mx * inv_hy) >> (BIT_WIDTH - sides_zeroes);=0A= +=A0=A0=A0 mx <<=3D sides_zeroes;=0A= +=A0=A0=A0 mx -=3D hd * my;=0A= +=A0=A0=A0 while (__glibc_unlikely (mx > my))=0A= +=A0=A0=A0=A0=A0 mx -=3D my;=0A= +=A0 }=0A= +=A0 uint32_t hd =3D (mx * inv_hy) >> (BIT_WIDTH - exp_diff);=0A= +=A0 mx <<=3D exp_diff;=0A= +=A0 mx -=3D hd * my;=0A= +=A0 while (__glibc_unlikely (mx > my))=0A= +=A0=A0=A0 mx -=3D my;=0A= +=0A= +=A0 return make_float (mx, ey, sx);=0A= =A0}=0A= =A0libm_alias_finite (__ieee754_fmodf, __fmodf)=0A= =0A= OK=0A= =0A= diff --git a/sysdeps/ieee754/flt-32/math_config.h b/sysdeps/ieee754/flt-32/= math_config.h=0A= index 23045f59d6..829430ea28 100644=0A= --- a/sysdeps/ieee754/flt-32/math_config.h=0A= +++ b/sysdeps/ieee754/flt-32/math_config.h=0A= @@ -110,6 +110,47 @@ issignalingf_inline (float x)=0A= =A0=A0 return 2 * (ix ^ 0x00400000) > 2 * 0x7fc00000UL;=0A= =A0}=0A= =A0=0A= +#define BIT_WIDTH=A0=A0=A0=A0=A0=A0 32=0A= +#define MANTISSA_WIDTH=A0 23=0A= +#define EXPONENT_WIDTH=A0 8=0A= +#define MANTISSA_MASK=A0=A0 0x007fffff=0A= +#define EXPONENT_MASK=A0=A0 0x7f800000=0A= +#define EXP_MANT_MASK=A0=A0 0x7fffffff=0A= +#define QUIET_NAN_MASK=A0 0x00400000=0A= +#define SIGN_MASK=A0=A0=A0=A0=A0=A0 0x80000000=0A= +=0A= +static inline bool=0A= +is_nan (uint32_t x)=0A= +{=0A= +=A0 return (x & EXP_MANT_MASK) > EXPONENT_MASK;=0A= +}=0A= +=0A= +static inline uint32_t=0A= +get_mantissa (uint32_t x)=0A= +{=0A= +=A0 return x & MANTISSA_MASK;=0A= +}=0A= +=0A= +/* Convert integer number X, unbiased exponent EP, and sign S to double:= =0A= +=0A= +=A0=A0 result =3D X * 2^(EP+1 - exponent_bias)=0A= +=0A= +=A0=A0 NB: zero is not supported.=A0 */=0A= +static inline double=0A= +make_float (uint32_t x, int ep, uint32_t s)=0A= +{=0A= +=A0 int lz =3D __builtin_clz (x) - EXPONENT_WIDTH;=0A= +=A0 x <<=3D lz;=0A= +=A0 ep -=3D lz;=0A= +=0A= +=A0 if (__glibc_unlikely (ep < 0 || x =3D=3D 0))=0A= +=A0=A0=A0 {=0A= +=A0=A0=A0=A0=A0 x >>=3D -ep;=0A= +=A0=A0=A0=A0=A0 ep =3D 0;=0A= +=A0=A0=A0 }=0A= +=A0 return asfloat (s + x + (ep << MANTISSA_WIDTH));=0A= +}=0A= =0A= OK=0A=