diff options
author | Gabriel F. T. Gomes <gabriel@inconstante.eti.br> | 2017-09-04 14:34:04 -0300 |
---|---|---|
committer | Gabriel F. T. Gomes <gabriel@inconstante.eti.br> | 2017-09-21 17:37:39 -0300 |
commit | c5c2e667bf834951ea95fc5124da4ef093711ebe (patch) | |
tree | ed57ab53a769264ba7762935ed73fdb53e61d2f7 /sysdeps/ieee754/ldbl-128ibm/s_cbrtl.c | |
parent | e010deb231113e92f5f351315645a8197857dffd (diff) | |
download | glibc-c5c2e667bf834951ea95fc5124da4ef093711ebe.tar.gz glibc-c5c2e667bf834951ea95fc5124da4ef093711ebe.tar.xz glibc-c5c2e667bf834951ea95fc5124da4ef093711ebe.zip |
ldbl-128ibm: Copy implementations from ldbl-128 instead of including them
Some files under sysdeps/ieee754/ldbl-128ibm/ are able to reuse the implementation in sysdeps/ieee754/ldbl-128/ by defining _Float128 to long double. This relied on compiler support for _Float128 being disabled. On powerpc, such support was disabled by default, however, it got enabled by default [1] in GCC 8. This patch copies the implementations from ldbl-128 to ldbl-128ibm. The uses of _Float128 and L() are kept intact in this patch and are replaced with a script in a subsequent patch. [1] https://gcc.gnu.org/ml/gcc-patches/2017-08/msg01028.html Tested for powerpc64 and powerpc64le. * sysdeps/ieee754/ldbl-128ibm/e_expl.c: Include tables from sysdeps/ieee754/ldbl-128ibm. * sysdeps/ieee754/ldbl-128ibm/e_j0l.c: Copy contents from the equivalent implementation in sysdeps/ieee754/ldbl-128/ instead of including it. Keep _Float128 and L() intact. These will be reviewed by a separate patch. * sysdeps/ieee754/ldbl-128ibm/e_j1l.c: Likewise. * sysdeps/ieee754/ldbl-128ibm/e_lgammal_r.c: Likewise. * sysdeps/ieee754/ldbl-128ibm/s_cbrtl.c: Likewise. * sysdeps/ieee754/ldbl-128ibm/t_expl.h: Likewise.
Diffstat (limited to 'sysdeps/ieee754/ldbl-128ibm/s_cbrtl.c')
-rw-r--r-- | sysdeps/ieee754/ldbl-128ibm/s_cbrtl.c | 104 |
1 files changed, 99 insertions, 5 deletions
diff --git a/sysdeps/ieee754/ldbl-128ibm/s_cbrtl.c b/sysdeps/ieee754/ldbl-128ibm/s_cbrtl.c index 64bfc46414..615707ca98 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_cbrtl.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_cbrtl.c @@ -1,10 +1,104 @@ -/* Looks like we can use ieee854 s_cbrtl.c as is for IBM extended format. */ -#include <math_ldbl_opt.h> -#undef weak_alias -#define weak_alias(n,a) +/* Implementation of cbrtl. IBM Extended Precision version. + Cephes Math Library Release 2.2: January, 1991 + Copyright 1984, 1991 by Stephen L. Moshier + Adapted for glibc October, 2001. + + This 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. + + This 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 this library; if not, see + <http://www.gnu.org/licenses/>. */ + +/* This file was copied from sysdeps/ieee754/ldbl-128/e_j0l.c. */ #define _Float128 long double #define L(x) x ## L -#include <sysdeps/ieee754/ldbl-128/s_cbrtl.c> +#include <math_ldbl_opt.h> +#include <math.h> +#include <math_private.h> + +static const _Float128 CBRT2 = L(1.259921049894873164767210607278228350570251); +static const _Float128 CBRT4 = L(1.587401051968199474751705639272308260391493); +static const _Float128 CBRT2I = L(0.7937005259840997373758528196361541301957467); +static const _Float128 CBRT4I = L(0.6299605249474365823836053036391141752851257); + + +_Float128 +__cbrtl (_Float128 x) +{ + int e, rem, sign; + _Float128 z; + + if (!isfinite (x)) + return x + x; + + if (x == 0) + return (x); + + if (x > 0) + sign = 1; + else + { + sign = -1; + x = -x; + } + + z = x; + /* extract power of 2, leaving mantissa between 0.5 and 1 */ + x = __frexpl (x, &e); + + /* Approximate cube root of number between .5 and 1, + peak relative error = 1.2e-6 */ + x = ((((L(1.3584464340920900529734e-1) * x + - L(6.3986917220457538402318e-1)) * x + + L(1.2875551670318751538055e0)) * x + - L(1.4897083391357284957891e0)) * x + + L(1.3304961236013647092521e0)) * x + L(3.7568280825958912391243e-1); + + /* exponent divided by 3 */ + if (e >= 0) + { + rem = e; + e /= 3; + rem -= 3 * e; + if (rem == 1) + x *= CBRT2; + else if (rem == 2) + x *= CBRT4; + } + else + { /* argument less than 1 */ + e = -e; + rem = e; + e /= 3; + rem -= 3 * e; + if (rem == 1) + x *= CBRT2I; + else if (rem == 2) + x *= CBRT4I; + e = -e; + } + + /* multiply by power of 2 */ + x = __ldexpl (x, e); + + /* Newton iteration */ + x -= (x - (z / (x * x))) * L(0.3333333333333333333333333333333333333333); + x -= (x - (z / (x * x))) * L(0.3333333333333333333333333333333333333333); + x -= (x - (z / (x * x))) * L(0.3333333333333333333333333333333333333333); + + if (sign < 0) + x = -x; + return (x); +} + long_double_symbol (libm, __cbrtl, cbrtl); |