From patchwork Thu Jun 18 16:44:47 2020 Content-Type: text/plain; charset="utf-8" MIME-Version: 1.0 Content-Transfer-Encoding: 7bit X-Patchwork-Submitter: Wilco Dijkstra X-Patchwork-Id: 39673 Return-Path: X-Original-To: patchwork@sourceware.org Delivered-To: patchwork@sourceware.org Received: from server2.sourceware.org (localhost [IPv6:::1]) by sourceware.org (Postfix) with ESMTP id 408D13870914; Thu, 18 Jun 2020 16:45:03 +0000 (GMT) X-Original-To: libc-alpha@sourceware.org Delivered-To: libc-alpha@sourceware.org Received: from EUR02-VE1-obe.outbound.protection.outlook.com (mail-eopbgr20043.outbound.protection.outlook.com [40.107.2.43]) by sourceware.org (Postfix) with ESMTPS id 016973870852 for ; Thu, 18 Jun 2020 16:44:59 +0000 (GMT) DMARC-Filter: OpenDMARC Filter v1.3.2 sourceware.org 016973870852 Authentication-Results: sourceware.org; dmarc=none (p=none dis=none) header.from=arm.com Authentication-Results: sourceware.org; spf=pass smtp.mailfrom=Wilco.Dijkstra@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=0EUXjoShY/iLrFYH7dNVotR6xmNGBBy0jsrN5oJQelw=; b=4KzKUihzbgBVpRHY68YKlVMiYeqkGKkAShO6LMrLP5i09bATzWFNwxxR30G/ixkyxV5QKArtPWVkqj+1rexBBiNY/dMy3ar8neZ0ADounCKRBwFye44tL35mnglXl/+pz39Y4HTXYXiG4Akaqj0w6s6YPtW4a1msAWlLeWvnyPk= Received: from MR2P264CA0002.FRAP264.PROD.OUTLOOK.COM (2603:10a6:500:1::14) by AM0PR08MB4276.eurprd08.prod.outlook.com (2603:10a6:208:13a::19) with Microsoft SMTP Server (version=TLS1_2, cipher=TLS_ECDHE_RSA_WITH_AES_256_GCM_SHA384) id 15.20.3109.22; Thu, 18 Jun 2020 16:44:56 +0000 Received: from VE1EUR03FT004.eop-EUR03.prod.protection.outlook.com (2603:10a6:500:1:cafe::aa) by MR2P264CA0002.outlook.office365.com (2603:10a6:500:1::14) with Microsoft SMTP Server (version=TLS1_2, cipher=TLS_ECDHE_RSA_WITH_AES_256_GCM_SHA384) id 15.20.3109.21 via Frontend Transport; Thu, 18 Jun 2020 16:44:56 +0000 X-MS-Exchange-Authentication-Results: spf=pass (sender IP is 63.35.35.123) smtp.mailfrom=arm.com; sourceware.org; dkim=pass (signature was verified) header.d=armh.onmicrosoft.com; sourceware.org; dmarc=bestguesspass 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; Received: from 64aa7808-outbound-1.mta.getcheckrecipient.com (63.35.35.123) by VE1EUR03FT004.mail.protection.outlook.com (10.152.18.106) with Microsoft SMTP Server (version=TLS1_2, cipher=TLS_ECDHE_RSA_WITH_AES_256_GCM_SHA384) id 15.20.3109.22 via Frontend Transport; Thu, 18 Jun 2020 16:44:56 +0000 Received: ("Tessian outbound cdc313860a50:v59"); Thu, 18 Jun 2020 16:44:56 +0000 X-CheckRecipientChecked: true X-CR-MTA-CID: 0e7714b202a9304e X-CR-MTA-TID: 64aa7808 Received: from 3ff9626143e8.1 by 64aa7808-outbound-1.mta.getcheckrecipient.com id 2824936B-27C8-4C85-B2C3-F99BEECFE3E9.1; Thu, 18 Jun 2020 16:44:50 +0000 Received: from EUR04-VI1-obe.outbound.protection.outlook.com by 64aa7808-outbound-1.mta.getcheckrecipient.com with ESMTPS id 3ff9626143e8.1 (version=TLSv1.2 cipher=ECDHE-RSA-AES256-GCM-SHA384); Thu, 18 Jun 2020 16:44:50 +0000 ARC-Seal: i=1; a=rsa-sha256; s=arcselector9901; d=microsoft.com; cv=none; b=nqf8BbSAi/1g6a7jLSqov+x/SrnMPV9D/XFKxlMhsIK7T+PU7FmY+AJASbhltjx53telWf/WcZf0j/EXehfCbaNVk7x5TZXERkjGpAesM7jFtC5YrIKeCIbmabtbTh78iEzIfkoKKsSYox1hCuReghJSw1RWlMYrrbq6JoTsPIuDU/yMB9hOs0e4xB8ClCooIxwWxzCSBbFiC0Fm7Oz0+jsk4N5oUpCLjKKQ0k+fQjD/Q80ACbl4IKcN6eJMJO6dR0PV8vrlG6VQPEScWZcJSbzZso5wWHfr4Z7yUNWqaz94p/nxBZ97Hah73/QV7oRkbN40BBUFJ+KyhZ6maqy1rw== 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-SenderADCheck; bh=0EUXjoShY/iLrFYH7dNVotR6xmNGBBy0jsrN5oJQelw=; b=W9mekcOifBtkmkqPnsZ48VxAqEkl0Wm3devpw5pT5ILGKSQR+5AqvhpBStUgLhj8iD3MHBdUbikKnSUVLBJE0nLgTGUP/UOVMnLHqhxfu6SOuwyXSYBZAwbuMNSvpEPhLBB23oQFBWZz7b2nNTzO0TbfMbFyzEgvgqIFmfCeJzEjS/y9RebkcomUonm28p3fWtI/jgpx/tf4YuRosGBdG85hse1iptSaDmmp2LurZ1eGMPIBLAZTZrblB6Msp9NOvIr0ynxOOVjXTJyhqL0WtQ4iiAzv3uHKPvn40cIBYX53aJ+N4bF5ru8ZN9fCg53h/40bVIm6yU3h+Uixlv8xNg== 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=0EUXjoShY/iLrFYH7dNVotR6xmNGBBy0jsrN5oJQelw=; b=4KzKUihzbgBVpRHY68YKlVMiYeqkGKkAShO6LMrLP5i09bATzWFNwxxR30G/ixkyxV5QKArtPWVkqj+1rexBBiNY/dMy3ar8neZ0ADounCKRBwFye44tL35mnglXl/+pz39Y4HTXYXiG4Akaqj0w6s6YPtW4a1msAWlLeWvnyPk= Received: from DB8PR08MB5036.eurprd08.prod.outlook.com (2603:10a6:10:ed::20) by DB8PR08MB4202.eurprd08.prod.outlook.com (2603:10a6:10:ae::17) with Microsoft SMTP Server (version=TLS1_2, cipher=TLS_ECDHE_RSA_WITH_AES_256_GCM_SHA384) id 15.20.3088.23; Thu, 18 Jun 2020 16:44:47 +0000 Received: from DB8PR08MB5036.eurprd08.prod.outlook.com ([fe80::40e3:3b43:9af2:d4ff]) by DB8PR08MB5036.eurprd08.prod.outlook.com ([fe80::40e3:3b43:9af2:d4ff%3]) with mapi id 15.20.3088.029; Thu, 18 Jun 2020 16:44:47 +0000 From: Wilco Dijkstra To: Adhemerval Zanella Subject: [PATCH 2/5] math: Optimized generic exp10f with wrappers Thread-Topic: [PATCH 2/5] math: Optimized generic exp10f with wrappers Thread-Index: AQHWRY1nPnWW6yR6sEa1W40++nPLXA== Date: Thu, 18 Jun 2020 16:44:47 +0000 Message-ID: Accept-Language: en-GB, en-US Content-Language: en-GB X-MS-Has-Attach: X-MS-TNEF-Correlator: Authentication-Results-Original: linaro.org; dkim=none (message not signed) header.d=none;linaro.org; dmarc=none action=none header.from=arm.com; x-originating-ip: [82.24.199.97] x-ms-publictraffictype: Email X-MS-Office365-Filtering-HT: Tenant X-MS-Office365-Filtering-Correlation-Id: c0a3e300-5768-4c9d-a753-08d813a6eea4 x-ms-traffictypediagnostic: DB8PR08MB4202:|AM0PR08MB4276: X-Microsoft-Antispam-PRVS: x-checkrecipientrouted: true nodisclaimer: true x-ms-oob-tlc-oobclassifiers: OLM:2657;OLM:2657; x-forefront-prvs: 0438F90F17 X-MS-Exchange-SenderADCheck: 1 X-Microsoft-Antispam-Untrusted: BCL:0; X-Microsoft-Antispam-Message-Info-Original: G48xoSmbjgcnKqQH+LopdUdRf882NwITyhH8lt7qGP4tBIqIYkiFSuqPynMegxwugUkK9meWKptnoJd+jDhb300GdJraicrzRzLR+jUUydhMCcJbiCn5CqdXjUYPGFjBNjBWXw+BcCS5bBpBSxhvt8MWfPLVXfZiTReAZGLb4K9/oVCrofDQOANvwMxZIAtqWYGft+9w1AS49SQlmImt1K/P8AI1GWdRqoJgl39rQ0Mlckb4adKZNKWunCev39ggaTL58c3Kpr7wGZcwIwAdzOc9a7nz8tQpYPtWO8MY0fb+Ns2v33OoVxDX+IG/0A9lirhrEMy5jnyLKUaTQAiPUFhL86zR21jyj8eLg0vyP575w6ocYFTEnTTsSDX+NX89okjJ6lxA5MhbxFU5uDcH9hbw86GachxHOlBfRYy2qiGoH08OOTA6dV241e3GupHyWv/Vryk8DkG2MaXh/baCoWLMUvzdrC32o3cQ+FdALa0= X-Forefront-Antispam-Report-Untrusted: CIP:255.255.255.255; CTRY:; LANG:en; SCL:1; SRV:; IPV:NLI; SFV:NSPM; H:DB8PR08MB5036.eurprd08.prod.outlook.com; PTR:; CAT:NONE; SFTY:; SFS:(4636009)(39860400002)(346002)(396003)(136003)(366004)(376002)(7696005)(86362001)(316002)(6506007)(478600001)(186003)(4326008)(26005)(54906003)(9686003)(66476007)(6916009)(52536014)(2906002)(76116006)(66946007)(83380400001)(8676002)(5660300002)(71200400001)(55016002)(64756008)(66556008)(66446008)(33656002)(8936002)(41533002)(2004002); DIR:OUT; SFP:1101; x-ms-exchange-antispam-messagedata: J4DfglNGbWLeuKyEF83EBnzHGbVOSq8hdkrQcUw0SYtpFVtVozgVMwIXUgx73DPpdLuzIBa1LhdMoheHgarial4ts4ii8ZEFiF+rbJbwlpdQKPsSEZJ95fy4oeZIl5MWk7n9/WrPy0ycx26Znbj8rbyOCuKDUDiNME2i8HO5kgI1AWGLjW+ptV5cqqd2La3QQHF6HptobJ5l/jlyzULAmJ94wP1oLCj3xahkExEGm3KpCQejZkz1Saec5+hvcEgx7BoTjndvc/TRfjrvYfm5/zcS7hwiMwQkOnGKTvZfPruOjZtDWqGTvBry+iBIfqHe5I2uZRhLaeSaTKOltFO+BFz1rJFaRI5QqURRpQGzpUVPwHvAoQZSRfLuP3CZEAn9Mo05EnmUAeo+lqyIfnkZLliPWYqMpwFo1Wk5pq3gDS0uduHluFPM8CrNgHQAWac/p6qExZkx9uVRkb979NOekAxCbnjs9ZDiIu6fNru2AfE= x-ms-exchange-transport-forked: True MIME-Version: 1.0 X-MS-Exchange-Transport-CrossTenantHeadersStamped: DB8PR08MB4202 Original-Authentication-Results: linaro.org; dkim=none (message not signed) header.d=none;linaro.org; dmarc=none action=none header.from=arm.com; X-EOPAttributedMessage: 0 X-MS-Exchange-Transport-CrossTenantHeadersStripped: VE1EUR03FT004.eop-EUR03.prod.protection.outlook.com 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; SFTY:; SFS:(4636009)(136003)(346002)(396003)(39860400002)(376002)(46966005)(8676002)(5660300002)(2906002)(52536014)(478600001)(8936002)(36906005)(4326008)(9686003)(55016002)(83380400001)(81166007)(107886003)(6862004)(316002)(54906003)(336012)(70586007)(70206006)(26005)(186003)(33656002)(7696005)(6506007)(86362001)(82740400003)(82310400002)(47076004)(356005)(41533002)(2004002); DIR:OUT; SFP:1101; X-MS-Office365-Filtering-Correlation-Id-Prvs: 7e1113bc-e813-4470-3dc7-08d813a6e99e X-Forefront-PRVS: 0438F90F17 X-Microsoft-Antispam: BCL:0; X-Microsoft-Antispam-Message-Info: Wtj0bvtpO/g+/FefTLpjmLsa9WEGaTKXlXNlez1zrcGpnix9IYqfkA7NuUU3vSN4J003EbEEyjBpwu5iHrn6L/ETn6JY+6sPiwow/EYzp0l7ZGZJGnOgzXLorajFr6fX9lfTHsTwEMWf1mI3RDeIYKSEaAlNqmFb7OhBPjc0qYdBsWAf3OqInP7FQnSgeqvA/fpexyBmzUDyGGrXSdToHHPSROLjhBH3BlRAKxOKgegO5t9DL6asmXscAppEbWtVOMiedKbMi6EEpsgdAzCg5RD0hbRQ1LogVcrGsO2kTQOO0QDqKKNazptKIUl3zHKR1leaKCHu7MsLUZqaBBslSbULUfWyx10UNAc+UAENQdejypOYNNdmGufXXee/klgVDwkaGMK5OwIPltTU87bAIU9U/fit6wTO2lLc5QMHVdirB1n7t+XzeuN5n8fP9eCmxujLln49ie9yTmEinMo/TAE9CXEFFWZvmRbTeBsJ87TSpq/Up4+7LzVQrpze8C0lGSolSbH2HHEUd3B8MJA6yQ== X-OriginatorOrg: arm.com X-MS-Exchange-CrossTenant-OriginalArrivalTime: 18 Jun 2020 16:44:56.0722 (UTC) X-MS-Exchange-CrossTenant-Network-Message-Id: c0a3e300-5768-4c9d-a753-08d813a6eea4 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-FromEntityHeader: HybridOnPrem X-MS-Exchange-Transport-CrossTenantHeadersStamped: AM0PR08MB4276 X-Spam-Status: No, score=-12.5 required=5.0 tests=BAYES_00, DKIM_SIGNED, DKIM_VALID, GIT_PATCH_0, KAM_SHORT, RCVD_IN_BARRACUDACENTRAL, RCVD_IN_DNSWL_LOW, RCVD_IN_MSPIKE_H2, SPF_HELO_PASS, SPF_PASS, TXREP, UNPARSEABLE_RELAY 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: libc-alpha@sourceware.org X-Mailman-Version: 2.1.29 Precedence: list List-Id: Libc-alpha mailing list List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , Cc: 'GNU C Library' Errors-To: libc-alpha-bounces@sourceware.org Sender: "Libc-alpha" Hi, This looks good (I like the extra comments with the bounds), OK with adding the missing 'f' on the float constant below. Wilco math/e_exp10f.c | 32 ----- sysdeps/ieee754/flt-32/e_exp10f.c | 198 +++++++++++++++++++++++++++ sysdeps/ieee754/flt-32/math_config.h | 2 +- 3 files changed, 199 insertions(+), 33 deletions(-) delete mode 100644 math/e_exp10f.c create mode 100644 sysdeps/ieee754/flt-32/e_exp10f.c OK diff --git a/math/e_exp10f.c b/math/e_exp10f.c deleted file mode 100644 index 93c41d00a6..0000000000 --- a/math/e_exp10f.c +++ /dev/null @@ -1,32 +0,0 @@ -/* Copyright (C) 1998-2020 Free Software Foundation, Inc. - This file is part of the GNU C Library. - Contributed by Ulrich Drepper , 1998. - - The GNU C Library is free software; you can redistribute it and/or - modify it under the terms of the GNU Lesser General Public - License as published by the Free Software Foundation; either - version 2.1 of the License, or (at your option) any later version. - - The GNU C Library is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - Lesser General Public License for more details. - - You should have received a copy of the GNU Lesser General Public - License along with the GNU C Library; if not, see - . */ - -#include -#include -#include - -float -__ieee754_exp10f (float arg) -{ - /* The argument to exp needs to be calculated in enough precision - that the fractional part has as much precision as float, in - addition to the bits in the integer part; using double ensures - this. */ - return __ieee754_exp (M_LN10 * arg); -} -libm_alias_finite (__ieee754_exp10f, __exp10f) OK diff --git a/sysdeps/ieee754/flt-32/e_exp10f.c b/sysdeps/ieee754/flt-32/e_exp10f.c new file mode 100644 index 0000000000..034a9e364a --- /dev/null +++ b/sysdeps/ieee754/flt-32/e_exp10f.c @@ -0,0 +1,198 @@ +/* Single-precision 10^x function. + Copyright (C) 2020 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with the GNU C Library; if not, see + . */ + +#include +#include +#include +#include +#include +#include "math_config.h" + +/* + EXP2F_TABLE_BITS 5 + EXP2F_POLY_ORDER 3 + + Max. ULP error: 0.502 (normal range, nearest rounding). + Max. relative error: 2^-33.240 (before rounding to float). + Wrong count: 169839. + Non-nearest ULP error: 1 (rounded ULP error). + + Detailed error analysis (for EXP2F_TABLE_BITS=3 thus N=32): + + - We first compute z = RN(InvLn10N * x) where + + InvLn10N = N*log(10)/log(2) * (1 + theta1) with |theta1| < 2^-54.150 + + since z is rounded to nearest: + + z = InvLn10N * x * (1 + theta2) with |theta2| < 2^-53 + + thus z = N*log(10)/log(2) * x * (1 + theta3) with |theta3| < 2^-52.463 + + - Since |x| < 38 in the main branch, we deduce: + + z = N*log(10)/log(2) * x + theta4 with |theta4| < 2^-40.483 + + - We then write z = k + r where k is an integer and |r| <= 0.5 (exact) + + - We thus have + + x = z*log(2)/(N*log(10)) - theta4*log(2)/(N*log(10)) + = z*log(2)/(N*log(10)) + theta5 with |theta5| < 2^-47.215 + + 10^x = 2^(k/N) * 2^(r/N) * 10^theta5 + = 2^(k/N) * 2^(r/N) * (1 + theta6) with |theta6| < 2^-46.011 + + - Then 2^(k/N) is approximated by table lookup, the maximal relative error + is for (k%N) = 5, with + + s = 2^(i/N) * (1 + theta7) with |theta7| < 2^-53.249 + + - 2^(r/N) is approximated by a degree-3 polynomial, where the maximal + mathematical relative error is 2^-33.243. + + - For C[0] * r + C[1], assuming no FMA is used, since |r| <= 0.5 and + |C[0]| < 1.694e-6, |C[0] * r| < 8.47e-7, and the absolute error on + C[0] * r is bounded by 1/2*ulp(8.47e-7) = 2^-74. Then we add C[1] with + |C[1]| < 0.000235, thus the absolute error on C[0] * r + C[1] is bounded + by 2^-65.994 (z is bounded by 0.000236). + + - For r2 = r * r, since |r| <= 0.5, we have |r2| <= 0.25 and the absolute + error is bounded by 1/4*ulp(0.25) = 2^-56 (the factor 1/4 is because |r2| + cannot exceed 1/4, and on the left of 1/4 the distance between two + consecutive numbers is 1/4*ulp(1/4)). + + - For y = C[2] * r + 1, assuming no FMA is used, since |r| <= 0.5 and + |C[2]| < 0.0217, the absolute error on C[2] * r is bounded by 2^-60, + and thus the absolute error on C[2] * r + 1 is bounded by 1/2*ulp(1)+2^60 + < 2^-52.988, and |y| < 1.01085 (the error bound is better if a FMA is + used). + + - for z * r2 + y: the absolute error on z is bounded by 2^-65.994, with + |z| < 0.000236, and the absolute error on r2 is bounded by 2^-56, with + r2 < 0.25, thus |z*r2| < 0.000059, and the absolute error on z * r2 + (including the rounding error) is bounded by: + + 2^-65.994 * 0.25 + 0.000236 * 2^-56 + 1/2*ulp(0.000059) < 2^-66.429. + + Now we add y, with |y| < 1.01085, and error on y bounded by 2^-52.988, + thus the absolute error is bounded by: + + 2^-66.429 + 2^-52.988 + 1/2*ulp(1.01085) < 2^-51.993 + + - Now we convert the error on y into relative error. Recall that y + approximates 2^(r/N), for |r| <= 0.5 and N=32. Thus 2^(-0.5/32) <= y, + and the relative error on y is bounded by + + 2^-51.993/2^(-0.5/32) < 2^-51.977 + + - Taking into account the mathematical relative error of 2^-33.243 we have: + + y = 2^(r/N) * (1 + theta8) with |theta8| < 2^-33.242 + + - Since we had s = 2^(k/N) * (1 + theta7) with |theta7| < 2^-53.249, + after y = y * s we get y = 2^(k/N) * 2^(r/N) * (1 + theta9) + with |theta9| < 2^-33.241 + + - Finally, taking into account the error theta6 due to the rounding error on + InvLn10N, and when multiplying InvLn10N * x, we get: + + y = 10^x * (1 + theta10) with |theta10| < 2^-33.240 + + - Converting into binary64 ulps, since y < 2^53*ulp(y), the error is at most + 2^19.76 ulp(y) + + - If the result is a binary32 value in the normal range (i.e., >= 2^-126), + then the error is at most 2^-9.24 ulps, i.e., less than 0.00166 (in the + subnormal range, the error in ulps might be larger). + + Note that this bound is tight, since for x=-0x25.54ac0p0 the final value of + y (before conversion to float) differs from 879853 ulps from the correctly + rounded value, and 879853 ~ 2^19.7469. */ + +#define N (1 << EXP2F_TABLE_BITS) + +#define InvLn10N (0x3.5269e12f346e2p0 * N) /* log(10)/log(2) to nearest */ +#define T __exp2f_data.tab +#define C __exp2f_data.poly_scaled +#define SHIFT __exp2f_data.shift + +static inline uint32_t +top13 (float x) +{ + return asuint (x) >> 19; +} + +float +__ieee754_exp10f (float x) +{ + uint32_t abstop; + uint64_t ki, t; + double kd, xd, z, r, r2, y, s; + + xd = (double) x; + abstop = top13 (x) & 0xfff; /* Ignore sign. */ + if (__glibc_unlikely (abstop >= top13 (38.0f))) + { + /* |x| >= 38 or x is nan. */ + if (asuint (x) == asuint (-INFINITY)) + return 0.0f; + if (abstop >= top13 (INFINITY)) + return x + x; + /* 0x26.8826ap0 is the largest value such that 10^x < 2^128. */ + if (x > 0x26.8826ap0f) + return __math_oflowf (0); + /* -0x2d.278d4p0 is the smallest value such that 10^x > 2^-150. */ + if (x < -0x2d.278d4p0f) + return __math_uflowf (0); +#if WANT_ERRNO_UFLOW + if (x < -0x2c.da7cfp0) This is missing an 'f' on the float constant. + return __math_may_uflowf (0); +#endif + /* the smallest value such that 10^x >= 2^-126 (normal range) + is x = -0x25.ee060p0 */ + /* we go through here for 2014929 values out of 2060451840 + (not counting NaN and infinities, i.e., about 0.1% */ + } + + /* x*N*Ln10/Ln2 = k + r with r in [-1/2, 1/2] and int k. */ + z = InvLn10N * xd; + /* |xd| < 38 thus |z| < 1216 */ +#if TOINT_INTRINSICS + kd = roundtoint (z); + ki = converttoint (z); +#else +# define SHIFT __exp2f_data.shift + kd = math_narrow_eval ((double) (z + SHIFT)); /* Needs to be double. */ + ki = asuint64 (kd); + kd -= SHIFT; +#endif + r = z - kd; + + /* 10^x = 10^(k/N) * 10^(r/N) ~= s * (C0*r^3 + C1*r^2 + C2*r + 1) */ + t = T[ki % N]; + t += ki << (52 - EXP2F_TABLE_BITS); + s = asdouble (t); + z = C[0] * r + C[1]; + r2 = r * r; + y = C[2] * r + 1; + y = z * r2 + y; + y = y * s; + return (float) y; +} +libm_alias_finite (__ieee754_exp10f, __exp10f) OK diff --git a/sysdeps/ieee754/flt-32/math_config.h b/sysdeps/ieee754/flt-32/math_config.h index bf79274ce7..4817e500e1 100644 --- a/sysdeps/ieee754/flt-32/math_config.h +++ b/sysdeps/ieee754/flt-32/math_config.h @@ -109,7 +109,7 @@ attribute_hidden float __math_may_uflowf (uint32_t); attribute_hidden float __math_divzerof (uint32_t); attribute_hidden float __math_invalidf (float); -/* Shared between expf, exp2f and powf. */ +/* Shared between expf, exp2f, exp10f, and powf. */ #define EXP2F_TABLE_BITS 5 #define EXP2F_POLY_ORDER 3 extern const struct exp2f_data