From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Received: from EUR05-AM6-obe.outbound.protection.outlook.com (mail-am6eur05on2070.outbound.protection.outlook.com [40.107.22.70]) by sourceware.org (Postfix) with ESMTPS id 87AA53858C5F for ; Mon, 3 Apr 2023 13:29:36 +0000 (GMT) DMARC-Filter: OpenDMARC Filter v1.4.2 sourceware.org 87AA53858C5F 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=EFJXGtZ1OyLaa63aoqj76/vwukLH/VqjKMVtlI9BpRw=; b=7mNzi8RpK9x/sj/GaKrcXMKvnffcZTdX84T6J5xF3w2iJhi96Dhtgciocfh9FTA1kGr7GmtK4O4yoRbZIJLu4PAzle27Jj0RmKDgEmUh1S1mfYvIIAX6OEQBJ6fKWuBZDFt/RKFr9fHpRAIi3Tbv2aic66S6y/CwfwtlMzJteX0= Received: from AS9PR0301CA0059.eurprd03.prod.outlook.com (2603:10a6:20b:469::21) by VI1PR08MB10174.eurprd08.prod.outlook.com (2603:10a6:800:1ca::6) 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:29:21 +0000 Received: from AM7EUR03FT058.eop-EUR03.prod.protection.outlook.com (2603:10a6:20b:469:cafe::f9) by AS9PR0301CA0059.outlook.office365.com (2603:10a6:20b:469::21) 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:29:21 +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 AM7EUR03FT058.mail.protection.outlook.com (100.127.140.247) with Microsoft SMTP Server (version=TLS1_2, cipher=TLS_ECDHE_RSA_WITH_AES_256_GCM_SHA384) id 15.20.6277.16 via Frontend Transport; Mon, 3 Apr 2023 13:29:20 +0000 Received: ("Tessian outbound 5bb4c51d5a1f:v136"); Mon, 03 Apr 2023 13:29:20 +0000 X-CheckRecipientChecked: true X-CR-MTA-CID: 844cc586876e06f6 X-CR-MTA-TID: 64aa7808 Received: from f8952af500cf.1 by 64aa7808-outbound-1.mta.getcheckrecipient.com id 3B48A2BC-220B-49FC-A902-7F2AD1847520.1; Mon, 03 Apr 2023 13:29:13 +0000 Received: from EUR04-DB3-obe.outbound.protection.outlook.com by 64aa7808-outbound-1.mta.getcheckrecipient.com with ESMTPS id f8952af500cf.1 (version=TLSv1.2 cipher=ECDHE-RSA-AES256-GCM-SHA384); Mon, 03 Apr 2023 13:29:13 +0000 ARC-Seal: i=1; a=rsa-sha256; s=arcselector9901; d=microsoft.com; cv=none; b=JTo7ibLbU3eUM3JGNahjXqZhqrzcvfQrYk1c5jzN0p11Rxj1jf6j/w5py4pAo1aXw31lrzvhApPcrxqHCw9zGE06LAjn5zHmRv6aTZ8b1BoOsuZ1B6UIow60wuPj/IlzDDHrg7QrXb1XwCSUMRBl4ZgNqFytSDkKXB61VschViCcR2pAGYgbDL4g2miHF2LpYBOCiwWt4db9U6RrLRVXaS2ohF5b8rp3sU+Wa05pH3tSS7DeW/JBxO7ZG50Re3lt8gN5QD821F1p+wYkCunBBFgjNB/qyuFSBzfyixhuPtHVG2wbgxIwJ6HpPdbXN1I+UHmBwv40a6SH4xxHrsc/hw== 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=EFJXGtZ1OyLaa63aoqj76/vwukLH/VqjKMVtlI9BpRw=; b=TrfvrVmf3y4Fvx9C4o76T/WmCXUfwQ5OLvvbceXpVMxcQkykpeuRIYB2Lu1WoEsoC7YMUggkNu2gETlR40I042vzEXCHhsVSnHHaIg7j3vIKwV+OPbmmFvEkHO5nI8BnZH5X7MGCHebOqe1MzwhZVh2Fjb2nCF/FApaE0xYkMseQKxftlHDClo6ZFxXix1wfm+x8/Ytj/PAvdY+5pupkPIhQS0vk8p4gas5ozIcV5Ffl6Y9L5tnJZhVDD7y7Hy3+UA0sti3VO4JQ9FUrX4ujlyh5I3m4Rv5uJW/SEeN924SHzr7BDD+MoN+a7ZLMjEV8edCHH2CMohBBPkjtIsH/og== 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=EFJXGtZ1OyLaa63aoqj76/vwukLH/VqjKMVtlI9BpRw=; b=7mNzi8RpK9x/sj/GaKrcXMKvnffcZTdX84T6J5xF3w2iJhi96Dhtgciocfh9FTA1kGr7GmtK4O4yoRbZIJLu4PAzle27Jj0RmKDgEmUh1S1mfYvIIAX6OEQBJ6fKWuBZDFt/RKFr9fHpRAIi3Tbv2aic66S6y/CwfwtlMzJteX0= Received: from PAWPR08MB8982.eurprd08.prod.outlook.com (2603:10a6:102:33f::20) by DU0PR08MB9348.eurprd08.prod.outlook.com (2603:10a6:10:420::7) 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:29:11 +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:29:10 +0000 From: Wilco Dijkstra To: Adhemerval Zanella , "libc-alpha@sourceware.org" , "H . J . Lu" CC: kirill Subject: Re: [PATCH v4 3/5] math: Improve fmod Thread-Topic: [PATCH v4 3/5] math: Improve fmod Thread-Index: AQHZW0VODHcYYNb/+0qeX4RZaWWL0a8Zpppf Date: Mon, 3 Apr 2023 13:29:10 +0000 Message-ID: References: <20230320160118.352206-1-adhemerval.zanella@linaro.org> <20230320160118.352206-4-adhemerval.zanella@linaro.org> In-Reply-To: <20230320160118.352206-4-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_|DU0PR08MB9348:EE_|AM7EUR03FT058:EE_|VI1PR08MB10174:EE_ X-MS-Office365-Filtering-Correlation-Id: 7b1f447c-440e-46b8-708f-08db34476ecf 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: h2KbFQKVtEc0ePPoG/+1Xr4QD3mcXuZTMQmaiah4jjLYsgV1l+XiTPnmcrLWOVK3Jvkv/tnPCD8qTwn6lyoHEq3klTYC47JR50IgzT/vAPvv3ifGPI488qxBx1sO1XYumeVQ27uPay+fgNFr6ndZziPi7n2huZBx0pe6cPD/ml0Wr2Yrg+JzOhNxTfKfF6bbHk+n3BDDp9fqS5qOIeXS5wzzr4egHu+UD8PYaCZYN1WIDXNgy+fj2T8DXKfQCkU4U3V05hrtSkERNo96sDXgXUOEeKuCLWJrUpZlxEUopqP8Bi1d36nlOuRcwhJp1kkIxOiROvcYO8TEItWsTYdCqhbGe2Eab0LdkqdOKWVsHcZT5AwbqHeedxVt+NvyV8UHmZh400mNZsK5XzrEi2wFLSV7q6/ODBkpy1jQRFVqWQ48fGiI43Q+r6l8Je7FiOZ7zmoS3Vh3cIPw8CyfcP4wxvPkXoQCm/Csile03yXf6lCl/9MyQ3rtaD3LQ6P7Qp7gLMzPO3KMJ5l5stiQVOdC3iOkC9uJW/BYhXIsknWsq4m6n9PtSw4m4U/AWnrkDSQUdh2Y/2bVFPWFsvHZhPFsNpFKtqRBV7w62erFUcRAfm9akGe1prF6bNnBouQl45QB 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)(39860400002)(396003)(366004)(376002)(346002)(136003)(451199021)(86362001)(38070700005)(33656002)(2906002)(7696005)(55016003)(71200400001)(83380400001)(186003)(26005)(6506007)(9686003)(4326008)(91956017)(8676002)(66476007)(76116006)(478600001)(66446008)(66556008)(66946007)(64756008)(41300700001)(38100700002)(122000001)(5660300002)(316002)(110136005)(8936002)(30864003)(52536014)(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: DU0PR08MB9348 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: AM7EUR03FT058.eop-EUR03.prod.protection.outlook.com X-MS-PublicTrafficType: Email X-MS-Office365-Filtering-Correlation-Id-Prvs: c0a42683-c7e6-4148-831c-08db344768a3 X-Microsoft-Antispam: BCL:0; X-Microsoft-Antispam-Message-Info: YB16mRZDK2+kexnuBIwj5TopNtFAtFEAr+onMbqW+kbNQmIx4n6UBG8Vp62Zo2HmOkO7N3iyChZwPQnLJBAyVkLwYxzx5ZaE3fzktdfwfi1gT8IjQlelfDGrrsyH5KaufuGKJq6vdcElcTzQXCyNZ69hApRjxYWyWfh+s+V9vL9TuONrUO0/1eKvl/nPATZKrPoylQrTfu4eIdK2y9KBp1eX/iqMzuGFKGlS3hciUoMfKJT3IZs/PKrv/ZGW0uUTwyCcqPnS582gikLe433MhFXBqcb+KqNDENJMbnwTSSUs6pO616SN1C0E6XMCCe0QbNKAdIM+/YdS22v3rXdWq4WXEo8YOshYE1zFOl01HEHDcEo6YqNeAZlAt3frHCPnPFArHXu6iKUWt57EEs4QZ0l+ZshyYix4DCUr9D64aFmN+Kyj4gzECeftseqGS6Z0pjVgyt7b8KEAC3R/De6zKJzp0jCyMPdLOPVv6FbiWUTdoR+S6kGu8/7FK3k6XcCb3yjmKDsJ5jBqtxBsWZ7VAMo0sp1wV7d0/Om/bk88QQZvQtDxiP34UZfoHYjdfAFDDDLnG+562iX75RL4JtAwpisJ3IEVrzcHPaaoj6GeIhgyFlSF/n9aUZY0QY6/mBL/u41P4UGzeb6j/ZZTu4vgI4tRizgtjaC/On/P4OGCMYcFw9iizyhbZ17oTz/dB2ZG 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)(346002)(376002)(396003)(39850400004)(136003)(451199021)(46966006)(36840700001)(356005)(316002)(110136005)(30864003)(8676002)(82310400005)(52536014)(70586007)(41300700001)(4326008)(70206006)(81166007)(5660300002)(8936002)(82740400003)(478600001)(86362001)(2906002)(7696005)(36860700001)(186003)(47076005)(55016003)(40480700001)(26005)(107886003)(9686003)(6506007)(33656002)(83380400001)(336012)(2004002);DIR:OUT;SFP:1101; X-OriginatorOrg: arm.com X-MS-Exchange-CrossTenant-OriginalArrivalTime: 03 Apr 2023 13:29:20.8686 (UTC) X-MS-Exchange-CrossTenant-Network-Message-Id: 7b1f447c-440e-46b8-708f-08db34476ecf 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: AM7EUR03FT058.eop-EUR03.prod.protection.outlook.com X-MS-Exchange-CrossTenant-AuthAs: Anonymous X-MS-Exchange-CrossTenant-FromEntityHeader: HybridOnPrem X-MS-Exchange-Transport-CrossTenantHeadersStamped: VI1PR08MB10174 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= =0A= diff --git a/math/libm-test-fmod.inc b/math/libm-test-fmod.inc=0A= index 39fd02f9ef..8841c1332c 100644=0A= --- a/math/libm-test-fmod.inc=0A= +++ b/math/libm-test-fmod.inc=0A= @@ -152,6 +152,16 @@ static const struct test_ff_f_data fmod_test_data[] = =3D=0A= =A0=A0=A0=A0 TEST_ff_f (fmod, 6.5, -2.25L, 2.0L, NO_INEXACT_EXCEPTION|ERRNO= _UNCHANGED),=0A= =A0=A0=A0=A0 TEST_ff_f (fmod, -6.5, -2.25L, -2.0L, NO_INEXACT_EXCEPTION|ERR= NO_UNCHANGED),=0A= =A0=0A= +=A0=A0=A0 TEST_ff_f (fmod, 0x1.e848p+18, 0x1.6p+3, 0x1.8p+2, NO_INEXACT_EX= CEPTION|ERRNO_UNCHANGED),=0A= +=A0=A0=A0 TEST_ff_f (fmod, 0x1.e848p+18, -0x1.6p+3, 0x1.8p+2, NO_INEXACT_E= XCEPTION|ERRNO_UNCHANGED),=0A= +=A0=A0=A0 TEST_ff_f (fmod, -0x1.e848p+18, 0x1.6p+3, -0x1.8p+2, NO_INEXACT_= EXCEPTION|ERRNO_UNCHANGED),=0A= +=A0=A0=A0 TEST_ff_f (fmod, -0x1.e848p+18, -0x1.6p+3, -0x1.8p+2, NO_INEXACT= _EXCEPTION|ERRNO_UNCHANGED),=0A= +=0A= +=A0=A0=A0 TEST_ff_f (fmod, 0x1.4p+2, 0x1p+0, 0x0p+0, NO_INEXACT_EXCEPTION|= ERRNO_UNCHANGED),=0A= +=A0=A0=A0 TEST_ff_f (fmod, 0x1.4p+2, -0x1p+0, 0x0p+0, NO_INEXACT_EXCEPTION= |ERRNO_UNCHANGED),=0A= +=A0=A0=A0 TEST_ff_f (fmod, -0x1.4p+2, 0x1p+0, -0x0p+0, NO_INEXACT_EXCEPTIO= N|ERRNO_UNCHANGED),=0A= +=A0=A0=A0 TEST_ff_f (fmod, -0x1.4p+2, -0x1p+0, -0x0p+0, NO_INEXACT_EXCEPTI= ON|ERRNO_UNCHANGED),=0A= +=0A= =A0=A0=A0=A0 TEST_ff_f (fmod, max_value, max_value, 0, NO_INEXACT_EXCEPTION= |ERRNO_UNCHANGED),=0A= =A0=A0=A0=A0 TEST_ff_f (fmod, max_value, -max_value, 0, NO_INEXACT_EXCEPTIO= N|ERRNO_UNCHANGED),=0A= =A0=A0=A0=A0 TEST_ff_f (fmod, max_value, min_value, 0, NO_INEXACT_EXCEPTION= |ERRNO_UNCHANGED),=0A= @@ -216,6 +226,10 @@ static const struct test_ff_f_data fmod_test_data[] = =3D=0A= =A0=A0=A0=A0 TEST_ff_f (fmod, -0x1p1023L, -0x3p-1073L, -0x1p-1073L, NO_INEX= ACT_EXCEPTION|ERRNO_UNCHANGED),=0A= =A0=A0=A0=A0 TEST_ff_f (fmod, -0x1p1023L, 0x3p-1022L, -0x1p-1021L, NO_INEXA= CT_EXCEPTION|ERRNO_UNCHANGED),=0A= =A0=A0=A0=A0 TEST_ff_f (fmod, -0x1p1023L, -0x3p-1022L, -0x1p-1021L, NO_INEX= ACT_EXCEPTION|ERRNO_UNCHANGED),=0A= +=A0=A0=A0 TEST_ff_f (fmod, 0x0.cded0a47373e9p-1022, 0x0.3212f5b8c8c16p-102= 2, 0x0.05a1336414391p-1022, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),=0A= +=A0=A0=A0 TEST_ff_f (fmod, 0x0.cded0a47373e9p-1022, -0x0.3212f5b8c8c16p-10= 22, 0x0.05a1336414391p-1022, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),=0A= +=A0=A0=A0 TEST_ff_f (fmod, -0x0.cded0a47373e9p-1022, 0x0.3212f5b8c8c16p-10= 22, -0x0.05a1336414391p-1022, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),=0A= +=A0=A0=A0 TEST_ff_f (fmod, -0x0.cded0a47373e9p-1022, -0x0.3212f5b8c8c16p-1= 022, -0x0.05a1336414391p-1022, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),=0A= =A0#endif=0A= =A0#if MIN_EXP <=3D -16381=0A= =A0=A0=A0=A0 TEST_ff_f (fmod, 0x1p16383L, 0x3p-16445L, 0x1p-16445L, NO_INEX= ACT_EXCEPTION|ERRNO_UNCHANGED),=0A= =0A= OK=0A= =0A= diff --git a/sysdeps/ieee754/dbl-64/e_fmod.c b/sysdeps/ieee754/dbl-64/e_fmo= d.c=0A= index 60b8bbb9d2..d20466fd5d 100644=0A= --- a/sysdeps/ieee754/dbl-64/e_fmod.c=0A= +++ b/sysdeps/ieee754/dbl-64/e_fmod.c=0A= @@ -1,105 +1,145 @@=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_fmod(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= -#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 double 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= =A0double=0A= =A0__ieee754_fmod (double x, double y)=0A= =A0{=0A= -=A0=A0=A0=A0=A0=A0 int32_t n,ix,iy;=0A= -=A0=A0=A0=A0=A0=A0 int64_t hx,hy,hz,sx,i;=0A= -=0A= -=A0=A0=A0=A0=A0=A0 EXTRACT_WORDS64(hx,x);=0A= -=A0=A0=A0=A0=A0=A0 EXTRACT_WORDS64(hy,y);=0A= -=A0=A0=A0=A0=A0=A0 sx =3D hx&UINT64_C(0x8000000000000000);=A0=A0 /* sign o= f 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= =A0=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 UINT64_C(0x7fffffffffffffff);=A0=A0=A0=A0 /* |y= | */=0A= -=0A= -=A0=A0=A0 /* purge off exception values */=0A= -=A0=A0=A0=A0=A0=A0 if(__builtin_expect(hy=3D=3D0=0A= -=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0= =A0=A0 || hx >=3D UINT64_C(0x7ff0000000000000)=0A= -=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0= =A0=A0 || hy > UINT64_C(0x7ff0000000000000), 0))=0A= -=A0=A0=A0=A0=A0=A0=A0=A0 /* y=3D0,or x not finite 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(__builtin_expect(hx<=3Dhy, 0)) {=0A= -=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0 if(hx>63];=A0=A0=A0=A0= =A0=A0 /* |x|=3D|y| return x*0*/=0A= -=A0=A0=A0=A0=A0=A0 }=0A= -=0A= -=A0=A0=A0 /* determine ix =3D ilogb(x) */=0A= -=A0=A0=A0=A0=A0=A0 if(__builtin_expect(hx0; i<<=3D1) ix = -=3D1;=0A= -=A0=A0=A0=A0=A0=A0 } else ix =3D (hx>>52)-1023;=0A= -=0A= -=A0=A0=A0 /* determine iy =3D ilogb(y) */=0A= -=A0=A0=A0=A0=A0=A0 if(__builtin_expect(hy0; i<<=3D1) iy = -=3D1;=0A= -=A0=A0=A0=A0=A0=A0 } else iy =3D (hy>>52)-1023;=0A= -=0A= -=A0=A0=A0 /* set up hx, hy and align y to x */=0A= -=A0=A0=A0=A0=A0=A0 if(__builtin_expect(ix >=3D -1022, 1))=0A= -=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0 hx =3D UINT64_C(0x0010000000000000)|(UINT64= _C(0x000fffffffffffff)&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 -1022-ix;=0A= -=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0 hx<<=3Dn;=0A= -=A0=A0=A0=A0=A0=A0 }=0A= -=A0=A0=A0=A0=A0=A0 if(__builtin_expect(iy >=3D -1022, 1))=0A= -=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0 hy =3D UINT64_C(0x0010000000000000)|(UINT64= _C(0x000fffffffffffff)&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 -1022-iy;=0A= -=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0 hy<<=3Dn;=0A= -=A0=A0=A0=A0=A0=A0 }=0A= -=0A= -=A0=A0=A0 /* fix point fmod */=0A= -=A0=A0=A0=A0=A0=A0 n =3D ix - iy;=0A= -=A0=A0=A0=A0=A0=A0 while(n--) {=0A= -=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0 hz=3Dhx-hy;=0A= -=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0 if(hz<0){hx =3D hx+hx;}=0A= -=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0 else {=0A= -=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0 if(hz=3D=3D0)=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=A0=A0=A0=A0=A0=A0=A0=A0 return Zero[(uint64= _t)sx>>63];=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[(uint64_t)sx>>63];=0A= -=A0=A0=A0=A0=A0=A0 while(hx=3D -1022, 1)) {=A0=A0 /* normal= ize output */=0A= -=A0=A0=A0=A0=A0=A0=A0=A0 hx =3D ((hx-UINT64_C(0x0010000000000000))|((uint6= 4_t)(iy+1023)<<52));=0A= -=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0 INSERT_WORDS64(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 -1022 - iy;=0A= -=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0 hx>>=3Dn;=0A= -=A0=A0=A0=A0=A0=A0=A0=A0=A0=A0 INSERT_WORDS64(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 uint64_t hx =3D asuint64 (x);=0A= +=A0 uint64_t hy =3D asuint64 (y);=0A= +=0A= +=A0 uint64_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 asdouble (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 -907 and |x/y| < 2^5= 2,=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 uint64_t d =3D (ex =3D=3D ey) ? (mx - my) : (mx << (ex - e= y)) % my;=0A= +=A0=A0=A0=A0=A0 return make_double (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 asdouble (sx | hx % hy);=0A= =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 uint64_t mx =3D get_mantissa (hx) | (MANTISSA_MASK + 1);=0A= +=A0 ex--;=0A= +=A0 uint64_t my =3D get_mantissa (hy) | (MANTISSA_MASK + 1);=0A= +=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 clz_uint64 (my);=0A= +=A0=A0=A0 }=0A= =0A= OK=0A= =0A= +=A0 /* Assume hy !=3D 0=A0 */=0A= +=A0 int tail_zeros_my =3D ctz_uint64 (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 asdouble (sx);=0A= =0A= =0A= OK=0A= =0A= +=A0 if (exp_diff =3D=3D 0)=0A= +=A0=A0=A0 return make_double (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 uint64_t inv_hy =3D UINT64_MAX / my;=0A= +=A0 while (exp_diff > sides_zeroes) {=0A= +=A0=A0=A0 exp_diff -=3D sides_zeroes;=0A= +=A0=A0=A0 uint64_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 uint64_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_double (mx, ey, sx);=0A= =A0}=0A= =A0libm_alias_finite (__ieee754_fmod, __fmod)=0A= diff --git a/sysdeps/ieee754/dbl-64/math_config.h b/sysdeps/ieee754/dbl-64/= math_config.h=0A= index 3cbaeede64..2049cea3f7 100644=0A= --- a/sysdeps/ieee754/dbl-64/math_config.h=0A= +++ b/sysdeps/ieee754/dbl-64/math_config.h=0A= @@ -43,6 +43,24 @@=0A= =A0# define TOINT_INTRINSICS 0=0A= =A0#endif=0A= =A0=0A= +static inline int=0A= +clz_uint64 (uint64_t x)=0A= +{=0A= +=A0 if (sizeof (uint64_t) =3D=3D sizeof (unsigned long))=0A= +=A0=A0=A0 return __builtin_clzl (x);=0A= +=A0 else=0A= +=A0=A0=A0 return __builtin_clzll (x);=0A= +}=0A= +=0A= +static inline int=0A= +ctz_uint64 (uint64_t x)=0A= +{=0A= +=A0 if (sizeof (uint64_t) =3D=3D sizeof (unsigned long))=0A= +=A0=A0=A0 return __builtin_ctzl (x);=0A= +=A0 else=0A= +=A0=A0=A0 return __builtin_ctzll (x);=0A= +}=0A= +=0A= =A0#if TOINT_INTRINSICS=0A= =A0/* Round x to nearest int in all rounding modes, ties have to be rounded= =0A= =A0=A0=A0 consistently with converttoint so the results match.=A0 If the re= sult=0A= @@ -88,6 +106,49 @@ issignaling_inline (double x)=0A= =A0=A0 return 2 * (ix ^ 0x0008000000000000) > 2 * 0x7ff8000000000000ULL;=0A= =A0}=0A= =A0=0A= +#define BIT_WIDTH=A0=A0=A0=A0=A0=A0 64=0A= +#define MANTISSA_WIDTH=A0 52=0A= +#define EXPONENT_WIDTH=A0 11=0A= +#define MANTISSA_MASK=A0=A0 UINT64_C(0x000fffffffffffff)=0A= +#define EXPONENT_MASK=A0=A0 UINT64_C(0x7ff0000000000000)=0A= +#define EXP_MANT_MASK=A0=A0 UINT64_C(0x7fffffffffffffff)=0A= +#define QUIET_NAN_MASK=A0 UINT64_C(0x0008000000000000)=0A= +#define SIGN_MASK=A0=A0=A0=A0=A0 UINT64_C(0x8000000000000000)=0A= +=0A= +static inline bool=0A= +is_nan (uint64_t x)=0A= +{=0A= +=A0 return (x & EXP_MANT_MASK) > EXPONENT_MASK;=0A= +}=0A= +=0A= +static inline uint64_t=0A= +get_mantissa (uint64_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_double (uint64_t x, int64_t ep, uint64_t s)=0A= +{=0A= +=A0 int lz =3D clz_uint64 (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= +=0A= +=A0 return asdouble (s + x + (ep << MANTISSA_WIDTH));=0A= +}=0A= =0A= OK=0A=