From patchwork Mon Sep 25 11:04:34 2017 Content-Type: text/plain; charset="utf-8" MIME-Version: 1.0 Content-Transfer-Encoding: 7bit X-Patchwork-Submitter: Szabolcs Nagy X-Patchwork-Id: 23129 Received: (qmail 22500 invoked by alias); 25 Sep 2017 11:04:48 -0000 Mailing-List: contact libc-alpha-help@sourceware.org; run by ezmlm Precedence: bulk List-Id: List-Unsubscribe: List-Subscribe: List-Archive: List-Post: List-Help: , Sender: libc-alpha-owner@sourceware.org Delivered-To: mailing list libc-alpha@sourceware.org Received: (qmail 21028 invoked by uid 89); 25 Sep 2017 11:04:45 -0000 Authentication-Results: sourceware.org; auth=none X-Virus-Found: No X-Spam-SWARE-Status: No, score=-24.0 required=5.0 tests=AWL, BAYES_00, GIT_PATCH_0, GIT_PATCH_1, GIT_PATCH_2, GIT_PATCH_3, KAM_ASCII_DIVIDERS, RCVD_IN_DNSWL_NONE, SPF_HELO_PASS, SPF_PASS autolearn=ham version=3.3.2 spammy=zj X-HELO: EUR02-HE1-obe.outbound.protection.outlook.com Authentication-Results: spf=none (sender IP is ) smtp.mailfrom=Szabolcs.Nagy@arm.com; Message-ID: <59C8E2C2.6000509@arm.com> Date: Mon, 25 Sep 2017 12:04:34 +0100 From: Szabolcs Nagy User-Agent: Mozilla/5.0 (X11; Linux x86_64; rv:31.0) Gecko/20100101 Thunderbird/31.8.0 MIME-Version: 1.0 To: GNU C Library CC: nd@arm.com Subject: [PATCH 4/7 v2] New generic powf References: <59C8E136.6070606@arm.com> In-Reply-To: <59C8E136.6070606@arm.com> X-ClientProxiedBy: DB6PR07CA0006.eurprd07.prod.outlook.com (2603:10a6:6:2d::16) To HE1PR0802MB2489.eurprd08.prod.outlook.com (2603:10a6:3:d8::23) X-MS-PublicTrafficType: Email X-MS-Office365-Filtering-Correlation-Id: cbcf463e-aeb6-4096-86a5-08d504053694 X-MS-Office365-Filtering-HT: Tenant X-Microsoft-Antispam: UriScan:; BCL:0; PCL:0; RULEID:(22001)(2017030254152)(48565401081)(2017052603199)(49563074)(201703131423075)(201703031133081)(201702281549075); SRVR:HE1PR0802MB2489; X-Microsoft-Exchange-Diagnostics: 1; HE1PR0802MB2489; 3:RlDU/sHedX/IZkz3N8hOcVllStZml09OMgnTR8CHzh8TRlptBe+uHnOI0BhG6jo8A8vJdUxQMlEqIRJKO/HBPo5eIbtTav1mfmUvqodnzbPHwUYcTTr7SgJGTJh4fEZwO8BefPoNbUQMEbnSi+ZohFboOMlTDflFs+2D9xtMoWHT7/YblMwHuYA+YIUfDLxepbupavfTXIczH9yxk4f7L9NzmX3MdECltQoym51NnMMkK8uGojaPzKYdLIKOtHMN; 25:2nmTgB89lSqYWxGePvYEIzMUJS1iOo4KZFCLO4M0X+q1txH85+BANIYchKSSlAHpQxlKTZwp2uhH3xjLHSC1XU21jmxz404/UBdyihGDsuEyopO0lKt1d1Cb3foso7vsymzcf7LjNIfzqKeAtXe9MzZ+GGInuZ50hOWaz6OGMU1ecU6jaZSZ0+4XSGHDihJbiehDkI4FZxoC4D5UbIwCQ7J0o4q+zQGO7x3PEt+qWRpMuBlTaChLiOfBeXNOIWbly9ufrXkM3eM0qfMoDLpS2NGIayPZWdGrF1HjYzzi9aSITbU6L5U8mp7tv6ooKqunGK9fKrVgzjSZMaGk6SNUCQ==; 31:ZPyy3425ycXt2jo9SzHagbSBS92Qv1mPnqE/9d8h5xMfTGX9EXRbqKbBKvda8w3MHyJjSbEVYp1in8ODzflJjD1UjBdf2Oo6CSdjjo09UGKsyy4nz8cSRP+/EhL6g/CobKcyIrAchk+15KT8pNS9RciaZ+0fgVW2CScsh6uxwX1RFjY3mI/qmAror4HzLpcyYYbq0QjGcF1P7IkxQBQDL3QzqDo+b3Rrgof4NqPKTJo= X-MS-TrafficTypeDiagnostic: HE1PR0802MB2489: NoDisclaimer: True X-Microsoft-Exchange-Diagnostics: 1; HE1PR0802MB2489; 20:O/H9f+1qQmQGSKZlSVeTflLCOde2Ha/9Joo8T+r0ETZq2FSPHLX5cqyIk7y/QLd6qJMwq5fmoN4BRjZuEC1V7U7vESBbp5IH6Q2w5ObKuEN8yK6juojd5AHtFcJM6NAq+/dBUqsqm1MuFrRAc99A7RNfgPs3bFwG6EJ0m6nr6wg=; 4:Yf2XUOrcVRjEDWwInNntXjOoWAZUckwdUKVyhjltCqmtNKvMF5IHakCtv8v9c4QhqHgXYZtEHjaqCLREpDMjLZZg3NiYobTl2fMMEUZOTAoQiHGRw1iglIUc1qvmZWLIBkfqHGUWCCgd7xTORhikUXrlU54YF82+9+ZQdqXAcPKEOoZqlf56KBS0x0dyQPZPoDLSrCNRkwY5rIojlOOyvPpiNrpxi4qyEJ8UyJCRvHrSL6uEzhhsT9lOUad8rM/e X-Exchange-Antispam-Report-Test: UriScan:; X-Microsoft-Antispam-PRVS: X-Exchange-Antispam-Report-CFA-Test: BCL:0; PCL:0; RULEID:(100000700101)(100105000095)(100000701101)(100105300095)(100000702101)(100105100095)(102415395)(6040450)(2401047)(8121501046)(5005006)(100000703101)(100105400095)(10201501046)(93006095)(93001095)(3002001)(6055026)(6041248)(20161123558100)(20161123560025)(20161123555025)(20161123562025)(20161123564025)(201703131423075)(201702281528075)(201703061421075)(201703061406153)(6072148)(201708071742011)(100000704101)(100105200095)(100000705101)(100105500095); SRVR:HE1PR0802MB2489; BCL:0; PCL:0; RULEID:(100000800101)(100110000095)(100000801101)(100110300095)(100000802101)(100110100095)(100000803101)(100110400095)(100000804101)(100110200095)(100000805101)(100110500095); SRVR:HE1PR0802MB2489; X-Forefront-PRVS: 04410E544A X-Forefront-Antispam-Report: SFV:NSPM; SFS:(10009020)(6009001)(6049001)(346002)(39860400002)(376002)(189002)(199003)(97736004)(16586007)(86362001)(68736007)(65956001)(66066001)(305945005)(65806001)(6916009)(33656002)(2950100002)(478600001)(36756003)(7736002)(72206003)(4610100001)(25786009)(54356999)(76176999)(65816999)(87266999)(5000100001)(59896002)(50986999)(8936002)(81166006)(81156014)(4326008)(64126003)(80316001)(2906002)(84326002)(83506001)(53936002)(58126008)(101416001)(8676002)(5660300001)(3846002)(270700001)(6116002)(568964002)(5890100001)(16576012)(16526017)(2476003)(316002)(6666003)(189998001)(105586002)(564344004)(106356001)(77096006)(6486002)(41533002); DIR:OUT; SFP:1101; SCL:1; SRVR:HE1PR0802MB2489; H:[10.2.206.69]; FPR:; SPF:None; PTR:InfoNoRecords; MX:1; A:1; LANG:en; Received-SPF: None (protection.outlook.com: arm.com does not designate permitted sender hosts) X-Microsoft-Exchange-Diagnostics: =?us-ascii?Q?1; HE1PR0802MB2489; 23:+GWHW4r06rEU6I4lGNHBJJGly77KTA2Nz6J+49B?= =?us-ascii?Q?SB1MN0Qf1Uh3jPfLvfZCMoE/RTMoO3tItccrR/awCWcVGZFv15olVF7Z2qG+?= =?us-ascii?Q?ml1P8NdiJHc1c5bAyr+TKi+Ek6M9agoOOOIY2dpGVTnNGzfmf04Y2SM/KQjI?= =?us-ascii?Q?c3gvQ9Y8A0yCkgb8OW38NsFSxYb6R02BLnxPCD8wwQccjq2rSBwJZQQ4z+0P?= =?us-ascii?Q?BUuuTtP1EFw8+vaqb2O8RDpeBMNUUhy67t0juY+Jv3rEUgwEIBqLoDZixUj9?= =?us-ascii?Q?HBi50nJipVOaSA6J6aDfKnCF8E3g1dZ3Fu2+URfrjx6APYGoHRtgaGlCwCFN?= =?us-ascii?Q?91dS3gtBYnv5BZgBzzJZ31xGLd+ynMy0KR7H2A1qFyMjyvtJzV4dBNphQDuD?= =?us-ascii?Q?mdMKlz7EDqkacg+KQ66MiA+VoCfW8EDw1uoEpb/kpuHnN3uuZZ7RZ4uCKo0q?= =?us-ascii?Q?WWa52IodbSHGIcBjwvf2EzSIK5mHgwfasfOXuP/Dr6JPYxp2YDbLv0joORsY?= =?us-ascii?Q?uHZ+ZYUU3YIJvOnlz/nh4flejFBlk85hSYRZ2uR1bfNPpsedCwkC0NO7Huw3?= =?us-ascii?Q?9ewc2Zzs9w6/GHTtaghk/d3TncxzlF5lqbahW2XlanZaHXFUe1x/pHTDMq7X?= =?us-ascii?Q?HF9ollKub0LTh4C2E9JDCDt6ySP52FEsfxsq3DdFhAKWpW34/VYpL/4NQPKS?= =?us-ascii?Q?ndooSWhQaMYgoH6zsEowbYfWBioKyet6cHnW75ccgmpqKWffYhRugIC5a4iI?= =?us-ascii?Q?sgDsRgHWxb5cSpg5wETSrLm7yKitior4WaSUczsdXVJ0lOcCfqWW27fIfT2W?= =?us-ascii?Q?IxavqHyHgeU8dwM94XIfXX6s9mk4+uoc2lamE489yMMvaawxJOUuNrrFTepb?= =?us-ascii?Q?lvEyPBKLHuYQRbPoYlsHvldeSEtj4VzcquHpuDSaVjPqemAMWHtQst1Rb3+I?= =?us-ascii?Q?VW41VCkApAL50+dJHan/vUixc5cjk9V98nygka+UxXdVBE+os+AI752fnros?= =?us-ascii?Q?xL5c1TFNJoSFdKoKXz4QqY273x7YXYm6iczQbRxhDatOb/XRy/gTLiwRlCMb?= =?us-ascii?Q?fCtDZ1LvII6mD6XQAU1u6KNJEEyfmFL/jjd6IzOsOqmIBq5WgIeveLNCmTEG?= =?us-ascii?Q?NMJRes0A/3I+qrzs8I3y2QFyu4eH2OBoxTqeQnF90JD6klZ/HGTQtoYAU80a?= =?us-ascii?Q?eTCv8IEvbWyQ8ZfpE3eViMcZq1n4po+xfPmlCntdafPIIE07PsYN1nncdCCe?= =?us-ascii?Q?YPgnSGniKFN0wvf6Bh19k/7WPEgTdYOl5IOCM3fHihG621c1e6xaVgW4zC0X?= =?us-ascii?Q?/8najnIkaESZx9OEhUYrqF8AA1gkqXBT5XAOEC0nKWVUJA9VUCPGCtf2ksvV?= =?us-ascii?Q?J8fcHon0kvEivzjXiOPZGDPrzaECh6zRcLE1+0o/ImOdMdN6I600Hf5RW2f2?= =?us-ascii?Q?XsxqAifgvye/OQPyV01tuQotCb9+fSbM=3D?= X-Microsoft-Exchange-Diagnostics: 1; HE1PR0802MB2489; 6:OoXY8R1PZiWifToQzOvd17RNM8wCTPAYxSTWUe/ZlQE4fAFCi10P4W/LG+XwSGBcpQ7sAyIhoyGzzAZw78p+3HUYGrFW063obnMVY7Zz2B7avgVtgmBhl/f2ARhb4dlkUTFwZ/ELWaMtFr1Mo+G3gr6x3BmGjm/PdmxYJQnI7ctAyrtZnQfo0afcHECALwcCBLS1/RSBY8OEzKPstddqlE27xaaKU1kFh3fa4rOD6gusZx4UDiTMNilFy577bVzz6afkeL0UWORlzVXwlLOs0zrTzC/fSmyDf/loLrFnadoG+4td5yiBqpsg14sMHC6KY4g/YQ3iqqvaEKxVazSKTQ==; 5:GcdjcyeJRcMeJUIaman6yP56Hdi9lRcxF7XPwLAaFxCELgQP7HzeKdIeYDyTFLb8VrsQFhTwL/Yu0rzb/QpLU47hjcnL+8mTB6mpLO2w1KzCR/N4Lsn++kO3DlNgYBGh0QXb4P7JxJr3QR1zDABNjw==; 24:5nilMl7dMdvJ1PBgAYosxcOhpWBdIbqd0b0+0/SjIIgHv0NATofOGijIJLrz63bRQ62lZ7mD11KlYLUxAkkeHjuVs69oS67f+EvoVP5WF6M=; 7:RWbYZ2IuzsyDwrJ+nhBMLdSV8MVI3AhFzIeVeCMkKWvG2GftS2SJjEOdmnWohvkBmi1rjlqCJKkX1Co991n64jbjHsbSwdEr05Iu7TTbP/dyS9Oa+gvr+He/7TA7ju5E6Pl437UtF322uQVKQSdocV2i9jKcekyGz3STjYiEM71pJP0YVHduRjzPIIpy/uk36NvI9EVyTAgKnqf1KrYRQwTlmocMruns+wYM8r2l8m4= SpamDiagnosticOutput: 1:99 SpamDiagnosticMetadata: NSPM X-OriginatorOrg: arm.com X-MS-Exchange-CrossTenant-OriginalArrivalTime: 25 Sep 2017 11:04:37.6250 (UTC) X-MS-Exchange-CrossTenant-FromEntityHeader: Hosted X-MS-Exchange-CrossTenant-Id: f34e5979-57d9-4aaa-ad4d-b122a662184d X-MS-Exchange-Transport-CrossTenantHeadersStamped: HE1PR0802MB2489 v2: - __powf_data is attribute_hidden From 22ad3727cd7ff290cee95423f8f74c265229e52d Mon Sep 17 00:00:00 2001 From: Szabolcs Nagy Date: Mon, 4 Sep 2017 17:55:33 +0100 Subject: [PATCH 4/7] New generic powf without wrapper on aarch64: powf reciprocal-throughput: 4.2x faster powf latency: 2.6x faster old worst-case error: 1.11 ulp new worst-case error: 0.82 ulp aarch64 .text size: -780 bytes aarch64 .rodata size: +144 bytes powf(x,y) is implemented as exp2(y*log2(x)) with the same algorithms that are used in exp2f and log2f, except that the log2f polynomial is larger for extra precision and its output (and exp2f input) may be scaled by a power of 2 (POWF_SCALE) to simplify the argument reduction step of exp2 (possible when efficient round and convert toint operation is available). The special case handling tries to minimize the checks in the hot path. When the input of exp2_inline is checked, int arithmetics is used as that was faster on the tested aarch64 cores. 2017-09-19 Szabolcs Nagy * math/Makefile (type-float-routines): Add e_powf_log2_data. * sysdeps/ieee754/flt-32/e_powf.c: New implementation. * sysdeps/ieee754/flt-32/e_powf_log2_data.c: New file. * sysdeps/ieee754/flt-32/math_config.h (__powf_data): Define. (issignalingf_inline): Likewise. (POWF_LOG2_TABLE_BITS): Likewise. (POWF_LOG2_POLY_ORDER): Likewise. (POWF_SCALE_BITS): Likewise. (POWF_SCALE): Likewise. * sysdeps/i386/fpu/e_powf_log2_data.c: New file. * sysdeps/ia64/fpu/e_powf_log2_data.c: New file. * sysdeps/m68k/m680x0/fpu/e_powf_log2_data.c: New file. --- math/Makefile | 2 +- sysdeps/i386/fpu/e_powf_log2_data.c | 1 + sysdeps/ia64/fpu/e_powf_log2_data.c | 1 + sysdeps/ieee754/flt-32/e_powf.c | 388 ++++++++++++++--------------- sysdeps/ieee754/flt-32/e_powf_log2_data.c | 45 ++++ sysdeps/ieee754/flt-32/math_config.h | 27 ++ sysdeps/m68k/m680x0/fpu/e_powf_log2_data.c | 1 + 7 files changed, 265 insertions(+), 200 deletions(-) create mode 100644 sysdeps/i386/fpu/e_powf_log2_data.c create mode 100644 sysdeps/ia64/fpu/e_powf_log2_data.c create mode 100644 sysdeps/ieee754/flt-32/e_powf_log2_data.c create mode 100644 sysdeps/m68k/m680x0/fpu/e_powf_log2_data.c diff --git a/math/Makefile b/math/Makefile index b4b3101592..6c8aa3e413 100644 --- a/math/Makefile +++ b/math/Makefile @@ -116,7 +116,7 @@ type-double-routines := branred doasin dosincos halfulp mpa mpatan2 \ # float support type-float-suffix := f type-float-routines := k_rem_pio2f math_errf e_exp2f_data e_logf_data \ - e_log2f_data + e_log2f_data e_powf_log2_data # _Float128 support type-float128-suffix := f128 diff --git a/sysdeps/i386/fpu/e_powf_log2_data.c b/sysdeps/i386/fpu/e_powf_log2_data.c new file mode 100644 index 0000000000..1cc8931700 --- /dev/null +++ b/sysdeps/i386/fpu/e_powf_log2_data.c @@ -0,0 +1 @@ +/* Not needed. */ diff --git a/sysdeps/ia64/fpu/e_powf_log2_data.c b/sysdeps/ia64/fpu/e_powf_log2_data.c new file mode 100644 index 0000000000..1cc8931700 --- /dev/null +++ b/sysdeps/ia64/fpu/e_powf_log2_data.c @@ -0,0 +1 @@ +/* Not needed. */ diff --git a/sysdeps/ieee754/flt-32/e_powf.c b/sysdeps/ieee754/flt-32/e_powf.c index ce8e11f1ea..644a18d05e 100644 --- a/sysdeps/ieee754/flt-32/e_powf.c +++ b/sysdeps/ieee754/flt-32/e_powf.c @@ -1,7 +1,5 @@ -/* e_powf.c -- float version of e_pow.c. - * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. - */ -/* Copyright (C) 2017 Free Software Foundation, Inc. +/* Single-precision pow function. + Copyright (C) 2017 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 @@ -18,210 +16,202 @@ License along with the GNU C Library; if not, see . */ -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - #include -#include - -static const float huge = 1.0e+30, tiny = 1.0e-30; - -static const float -bp[] = {1.0, 1.5,}, -zero = 0.0, -one = 1.0, -two = 2.0, -two24 = 16777216.0, /* 0x4b800000 */ - /* poly coefs for (3/2)*(log(x)-2s-2/3*s**3 */ -L1 = 6.0000002384e-01, /* 0x3f19999a */ -L2 = 4.2857143283e-01, /* 0x3edb6db7 */ -L3 = 3.3333334327e-01, /* 0x3eaaaaab */ -L4 = 2.7272811532e-01, /* 0x3e8ba305 */ -L5 = 2.3066075146e-01, /* 0x3e6c3255 */ -L6 = 2.0697501302e-01, /* 0x3e53f142 */ -P1 = 1.6666667163e-01, /* 0x3e2aaaab */ -P2 = -2.7777778450e-03, /* 0xbb360b61 */ -P3 = 6.6137559770e-05, /* 0x388ab355 */ -P4 = -1.6533901999e-06, /* 0xb5ddea0e */ -P5 = 4.1381369442e-08, /* 0x3331bb4c */ -ovt = 4.2995665694e-08; /* -(128-log2(ovfl+.5ulp)) */ - -static const double - dp[] = { 0.0, 0x1.2b803473f7ad1p-1, }, /* log2(1.5) */ - lg2 = M_LN2, - cp = 2.0/3.0/M_LN2, - invln2 = 1.0/M_LN2; +#include +#include "math_config.h" -float -__ieee754_powf(float x, float y) +/* +POWF_LOG2_POLY_ORDER = 5 +EXP2F_TABLE_BITS = 5 + +ULP error: 0.82 (~ 0.5 + relerr*2^24) +relerr: 1.27 * 2^-26 (Relative error ~= 128*Ln2*relerr_log2 + relerr_exp2) +relerr_log2: 1.83 * 2^-33 (Relative error of logx.) +relerr_exp2: 1.69 * 2^-34 (Relative error of exp2(ylogx).) +*/ + +#define N (1 << POWF_LOG2_TABLE_BITS) +#define T __powf_log2_data.tab +#define A __powf_log2_data.poly +#define OFF 0x3f330000 + +/* Subnormal input is normalized so ix has negative biased exponent. + Output is multiplied by N (POWF_SCALE) if TOINT_INTRINICS is set. */ +static inline double_t +log2_inline (uint32_t ix) { - float z, ax, s; - double d1, d2; - int32_t i,j,k,yisint,n; - int32_t hx,hy,ix,iy; - - GET_FLOAT_WORD(hy,y); - iy = hy&0x7fffffff; - - /* y==zero: x**0 = 1 */ - if(iy==0 && !issignaling (x)) return one; - - /* x==+-1 */ - if(x == 1.0 && !issignaling (y)) return one; - if(x == -1.0 && isinf(y)) return one; - - GET_FLOAT_WORD(hx,x); - ix = hx&0x7fffffff; - - /* +-NaN return x+y */ - if(__builtin_expect(ix > 0x7f800000 || - iy > 0x7f800000, 0)) - return x+y; - - /* special value of y */ - if (__builtin_expect(iy==0x7f800000, 0)) { /* y is +-inf */ - if (ix==0x3f800000) - return y - y; /* inf**+-1 is NaN */ - else if (ix > 0x3f800000)/* (|x|>1)**+-inf = inf,0 */ - return (hy>=0)? y: zero; - else /* (|x|<1)**-,+inf = inf,0 */ - return (hy<0)?-y: zero; - } - if(iy==0x3f800000) { /* y is +-1 */ - if(hy<0) return one/x; else return x; - } - if(hy==0x40000000) return x*x; /* y is 2 */ - if(hy==0x3f000000) { /* y is 0.5 */ - if(__builtin_expect(hx>=0, 1)) /* x >= +0 */ - return __ieee754_sqrtf(x); - } + /* double_t for better performance on targets with FLT_EVAL_METHOD==2. */ + double_t z, r, r2, r4, p, q, y, y0, invc, logc; + uint32_t iz, top, tmp; + int k, i; + + /* x = 2^k z; where z is in range [OFF,2*OFF] and exact. + The range is split into N subintervals. + The ith subinterval contains z and c is near its center. */ + tmp = ix - OFF; + i = (tmp >> (23 - POWF_LOG2_TABLE_BITS)) % N; + top = tmp & 0xff800000; + iz = ix - top; + k = (int32_t) top >> (23 - POWF_SCALE_BITS); /* arithmetic shift */ + invc = T[i].invc; + logc = T[i].logc; + z = (double_t) asfloat (iz); + + /* log2(x) = log1p(z/c-1)/ln2 + log2(c) + k */ + r = z * invc - 1; + y0 = logc + (double_t) k; + + /* Pipelined polynomial evaluation to approximate log1p(r)/ln2. */ + r2 = r * r; + y = A[0] * r + A[1]; + p = A[2] * r + A[3]; + r4 = r2 * r2; + q = A[4] * r + y0; + q = p * r2 + q; + y = y * r4 + q; + return y; +} - /* determine if y is an odd int when x < 0 - * yisint = 0 ... y is not an integer - * yisint = 1 ... y is an odd int - * yisint = 2 ... y is an even int - */ - yisint = 0; - if(hx<0) { - if(iy>=0x4b800000) yisint = 2; /* even integer y */ - else if(iy>=0x3f800000) { - k = (iy>>23)-0x7f; /* exponent */ - j = iy>>(23-k); - if((j<<(23-k))==iy) yisint = 2-(j&1); - } - } +#undef N +#undef T +#define N (1 << EXP2F_TABLE_BITS) +#define T __exp2f_data.tab +#define SIGN_BIAS (1 << (EXP2F_TABLE_BITS + 11)) + +/* The output of log2 and thus the input of exp2 is either scaled by N + (in case of fast toint intrinsics) or not. The unscaled xd must be + in [-1021,1023], sign_bias sets the sign of the result. */ +static inline double_t +exp2_inline (double_t xd, unsigned long sign_bias) +{ + uint64_t ki, ski, t; + /* double_t for better performance on targets with FLT_EVAL_METHOD==2. */ + double_t kd, z, r, r2, y, s; + +#if TOINT_INTRINSICS +# define C __exp2f_data.poly_scaled + /* N*x = k + r with r in [-1/2, 1/2] */ + kd = roundtoint (xd); /* k */ + ki = converttoint (xd); +#else +# define C __exp2f_data.poly +# define SHIFT __exp2f_data.shift_scaled + /* x = k/N + r with r in [-1/(2N), 1/(2N)] */ + kd = (double) (xd + SHIFT); /* Rounding to double precision is required. */ + ki = asuint64 (kd); + kd -= SHIFT; /* k/N */ +#endif + r = xd - kd; + + /* exp2(x) = 2^(k/N) * 2^r ~= s * (C0*r^3 + C1*r^2 + C2*r + 1) */ + t = T[ki % N]; + ski = ki + sign_bias; + t += ski << (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 y; +} - ax = fabsf(x); - /* special value of x */ - if(__builtin_expect(ix==0x7f800000||ix==0||ix==0x3f800000, 0)){ - z = ax; /*x is +-0,+-inf,+-1*/ - if(hy<0) z = one/z; /* z = (1/|x|) */ - if(hx<0) { - if(((ix-0x3f800000)|yisint)==0) { - z = (z-z)/(z-z); /* (-1)**non-int is NaN */ - } else if(yisint==1) - z = -z; /* (x<0)**odd = -(|x|**odd) */ - } - return z; - } +/* Returns 0 if not int, 1 if odd int, 2 if even int. */ +static inline int +checkint (uint32_t iy) +{ + int e = iy >> 23 & 0xff; + if (e < 0x7f) + return 0; + if (e > 0x7f + 23) + return 2; + if (iy & ((1 << (0x7f + 23 - e)) - 1)) + return 0; + if (iy & (1 << (0x7f + 23 - e))) + return 1; + return 2; +} - /* (x<0)**(non-int) is NaN */ - if(__builtin_expect(((((uint32_t)hx>>31)-1)|yisint)==0, 0)) - return (x-x)/(x-x); - - /* |y| is huge */ - if(__builtin_expect(iy>0x4d000000, 0)) { /* if |y| > 2**27 */ - /* over/underflow if x is not close to one */ - if(ix<0x3f7ffff8) return (hy<0)? huge*huge:tiny*tiny; - if(ix>0x3f800007) return (hy>0)? huge*huge:tiny*tiny; - /* now |1-x| is tiny <= 2**-20, suffice to compute - log(x) by x-x^2/2+x^3/3-x^4/4 */ - d2 = ax-1; /* d2 has 20 trailing zeros. */ - d2 = d2 * invln2 - - (d2 * d2) * (0.5 - d2 * (0.333333333333 - d2 * 0.25)) * invln2; - } else { - /* Avoid internal underflow for tiny y. The exact value - of y does not matter if |y| <= 2**-32. */ - if (iy < 0x2f800000) - SET_FLOAT_WORD (y, (hy & 0x80000000) | 0x2f800000); - n = 0; - /* take care subnormal number */ - if(ix<0x00800000) - {ax *= two24; n -= 24; GET_FLOAT_WORD(ix,ax); } - n += ((ix)>>23)-0x7f; - j = ix&0x007fffff; - /* determine interval */ - ix = j|0x3f800000; /* normalize ix */ - if(j<=0x1cc471) k=0; /* |x|= 2u * 0x7f800000 - 1; +} - s = one; /* s (sign of result -ve**odd) = -1 else = 1 */ - if(((((uint32_t)hx>>31)-1)|(yisint-1))==0) - s = -one; /* (-ve)**(odd int) */ - - /* compute y * d2 */ - d1 = y * d2; - z = d1; - GET_FLOAT_WORD(j,z); - if (__builtin_expect(j>0x43000000, 0)) /* if z > 128 */ - return s*huge*huge; /* overflow */ - else if (__builtin_expect(j==0x43000000, 0)) { /* if z == 128 */ - if(ovt>(z-d1)) return s*huge*huge; /* overflow */ +float +__ieee754_powf (float x, float y) +{ + unsigned long sign_bias = 0; + uint32_t ix, iy; + + ix = asuint (x); + iy = asuint (y); + if (__glibc_unlikely (ix - 0x00800000 >= 0x7f800000 - 0x00800000 + || zeroinfnan (iy))) + { + /* Either (x < 0x1p-126 or inf or nan) or (y is 0 or inf or nan). */ + if (__glibc_unlikely (zeroinfnan (iy))) + { + if (2 * iy == 0) + return issignalingf_inline (x) ? x + y : 1.0f; + if (ix == 0x3f800000) + return issignalingf_inline (y) ? x + y : 1.0f; + if (2 * ix > 2u * 0x7f800000 || 2 * iy > 2u * 0x7f800000) + return x + y; + if (2 * ix == 2 * 0x3f800000) + return 1.0f; + if ((2 * ix < 2 * 0x3f800000) == !(iy & 0x80000000)) + return 0.0f; /* |x|<1 && y==inf or |x|>1 && y==-inf. */ + return y * y; + } + if (__glibc_unlikely (zeroinfnan (ix))) + { + float_t x2 = x * x; + if (ix & 0x80000000 && checkint (iy) == 1) + { + x2 = -x2; + sign_bias = 1; + } +#if WANT_ERRNO + if (2 * ix == 0 && iy & 0x80000000) + return __math_divzerof (sign_bias); +#endif + return iy & 0x80000000 ? 1 / x2 : x2; } - else if (__builtin_expect((j&0x7fffffff)>0x43160000, 0))/* z <= -150 */ - return s*tiny*tiny; /* underflow */ - else if (__builtin_expect((uint32_t) j==0xc3160000, 0)){/* z == -150*/ - if(0.0<=(z-d1)) return s*tiny*tiny; /* underflow */ + /* x and y are non-zero finite. */ + if (ix & 0x80000000) + { + /* Finite x < 0. */ + int yint = checkint (iy); + if (yint == 0) + return __math_invalidf (x); + if (yint == 1) + sign_bias = SIGN_BIAS; + ix &= 0x7fffffff; } - /* - * compute 2**d1 - */ - i = j&0x7fffffff; - k = (i>>23)-0x7f; - n = 0; - if(i>0x3f000000) { /* if |z| > 0.5, set n = [z+0.5] */ - n = j+(0x00800000>>(k+1)); - k = ((n&0x7fffffff)>>23)-0x7f; /* new k for n */ - SET_FLOAT_WORD(z,n&~(0x007fffff>>k)); - n = ((n&0x007fffff)|0x00800000)>>(23-k); - if(j<0) n = -n; - d1 -= z; + if (ix < 0x00800000) + { + /* Normalize subnormal x so exponent becomes negative. */ + ix = asuint (x * 0x1p23f); + ix &= 0x7fffffff; + ix -= 23 << 23; } - d1 = d1 * lg2; - d2 = d1*d1; - d2 = d1 - d2*(P1+d2*(P2+d2*(P3+d2*(P4+d2*P5)))); - d2 = (d1*d2)/(d2-two); - z = one - (d2-d1); - GET_FLOAT_WORD(j,z); - j += (n<<23); - if((j>>23)<=0) /* subnormal output */ - { - z = __scalbnf (z, n); - float force_underflow = z * z; - math_force_eval (force_underflow); - } - else SET_FLOAT_WORD(z,j); - return s*z; + } + double_t logx = log2_inline (ix); + double_t ylogx = y * logx; /* Note: cannot overflow, y is single prec. */ + if (__glibc_unlikely ((asuint64 (ylogx) >> 47 & 0xffff) + >= asuint64 (126.0 * POWF_SCALE) >> 47)) + { + /* |y*log(x)| >= 126. */ + if (ylogx > 0x1.fffffffd1d571p+6 * POWF_SCALE) + return __math_oflowf (sign_bias); + if (ylogx <= -150.0 * POWF_SCALE) + return __math_uflowf (sign_bias); +#if WANT_ERRNO_UFLOW + if (ylogx < -149.0 * POWF_SCALE) + return __math_may_uflowf (sign_bias); +#endif + } + return (float) exp2_inline (ylogx, sign_bias); } strong_alias (__ieee754_powf, __powf_finite) diff --git a/sysdeps/ieee754/flt-32/e_powf_log2_data.c b/sysdeps/ieee754/flt-32/e_powf_log2_data.c new file mode 100644 index 0000000000..7cff06f59b --- /dev/null +++ b/sysdeps/ieee754/flt-32/e_powf_log2_data.c @@ -0,0 +1,45 @@ +/* Data definition for powf. + Copyright (C) 2017 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 "math_config.h" + +const struct powf_log2_data __powf_log2_data = { + .tab = { + { 0x1.661ec79f8f3bep+0, -0x1.efec65b963019p-2 * POWF_SCALE }, + { 0x1.571ed4aaf883dp+0, -0x1.b0b6832d4fca4p-2 * POWF_SCALE }, + { 0x1.49539f0f010bp+0, -0x1.7418b0a1fb77bp-2 * POWF_SCALE }, + { 0x1.3c995b0b80385p+0, -0x1.39de91a6dcf7bp-2 * POWF_SCALE }, + { 0x1.30d190c8864a5p+0, -0x1.01d9bf3f2b631p-2 * POWF_SCALE }, + { 0x1.25e227b0b8eap+0, -0x1.97c1d1b3b7afp-3 * POWF_SCALE }, + { 0x1.1bb4a4a1a343fp+0, -0x1.2f9e393af3c9fp-3 * POWF_SCALE }, + { 0x1.12358f08ae5bap+0, -0x1.960cbbf788d5cp-4 * POWF_SCALE }, + { 0x1.0953f419900a7p+0, -0x1.a6f9db6475fcep-5 * POWF_SCALE }, + { 0x1p+0, 0x0p+0 * POWF_SCALE }, + { 0x1.e608cfd9a47acp-1, 0x1.338ca9f24f53dp-4 * POWF_SCALE }, + { 0x1.ca4b31f026aap-1, 0x1.476a9543891bap-3 * POWF_SCALE }, + { 0x1.b2036576afce6p-1, 0x1.e840b4ac4e4d2p-3 * POWF_SCALE }, + { 0x1.9c2d163a1aa2dp-1, 0x1.40645f0c6651cp-2 * POWF_SCALE }, + { 0x1.886e6037841edp-1, 0x1.88e9c2c1b9ff8p-2 * POWF_SCALE }, + { 0x1.767dcf5534862p-1, 0x1.ce0a44eb17bccp-2 * POWF_SCALE }, + }, + .poly = { + 0x1.27616c9496e0bp-2 * POWF_SCALE, -0x1.71969a075c67ap-2 * POWF_SCALE, + 0x1.ec70a6ca7baddp-2 * POWF_SCALE, -0x1.7154748bef6c8p-1 * POWF_SCALE, + 0x1.71547652ab82bp0 * POWF_SCALE, + } +}; diff --git a/sysdeps/ieee754/flt-32/math_config.h b/sysdeps/ieee754/flt-32/math_config.h index cca0b53d36..aeff4808c6 100644 --- a/sysdeps/ieee754/flt-32/math_config.h +++ b/sysdeps/ieee754/flt-32/math_config.h @@ -21,6 +21,7 @@ #include #include +#include #include #ifndef WANT_ROUNDING @@ -90,6 +91,15 @@ asdouble (uint64_t i) return u.f; } +static inline int +issignalingf_inline (float x) +{ + uint32_t ix = asuint (x); + if (HIGH_ORDER_BIT_IS_SET_FOR_SNAN) + return (ix & 0x7fc00000) == 0x7fc00000; + return 2 * (ix ^ 0x00400000) > 2u * 0x7fc00000; +} + #define NOINLINE __attribute__ ((noinline)) attribute_hidden float __math_oflowf (unsigned long); @@ -134,4 +144,21 @@ extern const struct log2f_data double poly[LOG2F_POLY_ORDER]; } __log2f_data attribute_hidden; +#define POWF_LOG2_TABLE_BITS 4 +#define POWF_LOG2_POLY_ORDER 5 +#if TOINT_INTRINSICS +#define POWF_SCALE_BITS EXP2F_TABLE_BITS +#else +#define POWF_SCALE_BITS 0 +#endif +#define POWF_SCALE ((double) (1 << POWF_SCALE_BITS)) +extern const struct powf_log2_data +{ + struct + { + double invc, logc; + } tab[1 << POWF_LOG2_TABLE_BITS]; + double poly[POWF_LOG2_POLY_ORDER]; +} __powf_log2_data attribute_hidden; + #endif diff --git a/sysdeps/m68k/m680x0/fpu/e_powf_log2_data.c b/sysdeps/m68k/m680x0/fpu/e_powf_log2_data.c new file mode 100644 index 0000000000..1cc8931700 --- /dev/null +++ b/sysdeps/m68k/m680x0/fpu/e_powf_log2_data.c @@ -0,0 +1 @@ +/* Not needed. */ -- 2.11.0