mirror of
https://sourceware.org/git/glibc.git
synced 2024-11-23 13:30:06 +00:00
60f435bb0c
Various sysdeps/ieee754/dbl-64 functions use double constants defined using a union between a double and two ints, with separate big-endian and little-endian definitions of the constants. With modern C, this is unnecessary complication; hex float constants (or __builtin_inf etc.) suffice to specify the exact value desired, and so can avoid separate versions for each endianness. Having this complication also complicates cleanups such as removing slow paths from these library functions, as they need to make sure to remove both copies of variables that are no longer used after such a cleanup (and in at least one case, proper removal of a slow path will also involve removing slow-path-only values from the middle of an array - an array with both big-endian and little-endian copies - and adjusting other references to that array). So it makes sense to clean up the code to define these constants using hex floats and so eliminate the endianness conditional. This patch does so in the case of sqrt, where the two constants are such that it makes sense just to put them directly in the code using them and eliminate the names for them altogether. Tested for arm (the code generated for sqrt does change, though not in any significant way). * sysdeps/ieee754/dbl-64/e_sqrt.c: Do not include uroot.h. (__ieee754_sqrt): Use hex float constants instead of tm256.x and t512.x. * sysdeps/ieee754/dbl-64/uroot.h: Remove file.
140 lines
4.8 KiB
C
140 lines
4.8 KiB
C
/*
|
|
* IBM Accurate Mathematical Library
|
|
* written by International Business Machines Corp.
|
|
* Copyright (C) 2001-2015 Free Software Foundation, Inc.
|
|
*
|
|
* This program 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 program 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 program; if not, see <http://www.gnu.org/licenses/>.
|
|
*/
|
|
/*********************************************************************/
|
|
/* MODULE_NAME: uroot.c */
|
|
/* */
|
|
/* FUNCTION: usqrt */
|
|
/* */
|
|
/* FILES NEEDED: dla.h endian.h mydefs.h */
|
|
/* uroot.tbl */
|
|
/* */
|
|
/* An ultimate sqrt routine. Given an IEEE double machine number x */
|
|
/* it computes the correctly rounded (to nearest) value of square */
|
|
/* root of x. */
|
|
/* Assumption: Machine arithmetic operations are performed in */
|
|
/* round to nearest mode of IEEE 754 standard. */
|
|
/* */
|
|
/*********************************************************************/
|
|
|
|
#include "endian.h"
|
|
#include "mydefs.h"
|
|
#include <dla.h>
|
|
#include "MathLib.h"
|
|
#include "root.tbl"
|
|
#include <math_private.h>
|
|
|
|
/*********************************************************************/
|
|
/* An ultimate sqrt routine. Given an IEEE double machine number x */
|
|
/* it computes the correctly rounded (to nearest) value of square */
|
|
/* root of x. */
|
|
/*********************************************************************/
|
|
double
|
|
__ieee754_sqrt (double x)
|
|
{
|
|
static const double
|
|
rt0 = 9.99999999859990725855365213134618E-01,
|
|
rt1 = 4.99999999495955425917856814202739E-01,
|
|
rt2 = 3.75017500867345182581453026130850E-01,
|
|
rt3 = 3.12523626554518656309172508769531E-01;
|
|
static const double big = 134217728.0;
|
|
double y, t, del, res, res1, hy, z, zz, p, hx, tx, ty, s;
|
|
mynumber a, c = { { 0, 0 } };
|
|
int4 k;
|
|
|
|
a.x = x;
|
|
k = a.i[HIGH_HALF];
|
|
a.i[HIGH_HALF] = (k & 0x001fffff) | 0x3fe00000;
|
|
t = inroot[(k & 0x001fffff) >> 14];
|
|
s = a.x;
|
|
/*----------------- 2^-1022 <= | x |< 2^1024 -----------------*/
|
|
if (k > 0x000fffff && k < 0x7ff00000)
|
|
{
|
|
int rm = __fegetround ();
|
|
fenv_t env;
|
|
libc_feholdexcept_setround (&env, FE_TONEAREST);
|
|
double ret;
|
|
y = 1.0 - t * (t * s);
|
|
t = t * (rt0 + y * (rt1 + y * (rt2 + y * rt3)));
|
|
c.i[HIGH_HALF] = 0x20000000 + ((k & 0x7fe00000) >> 1);
|
|
y = t * s;
|
|
hy = (y + big) - big;
|
|
del = 0.5 * t * ((s - hy * hy) - (y - hy) * (y + hy));
|
|
res = y + del;
|
|
if (res == (res + 1.002 * ((y - res) + del)))
|
|
ret = res * c.x;
|
|
else
|
|
{
|
|
res1 = res + 1.5 * ((y - res) + del);
|
|
EMULV (res, res1, z, zz, p, hx, tx, hy, ty); /* (z+zz)=res*res1 */
|
|
res = ((((z - s) + zz) < 0) ? max (res, res1) :
|
|
min (res, res1));
|
|
ret = res * c.x;
|
|
}
|
|
math_force_eval (ret);
|
|
libc_fesetenv (&env);
|
|
double dret = x / ret;
|
|
if (dret != ret)
|
|
{
|
|
double force_inexact = 1.0 / 3.0;
|
|
math_force_eval (force_inexact);
|
|
/* The square root is inexact, ret is the round-to-nearest
|
|
value which may need adjusting for other rounding
|
|
modes. */
|
|
switch (rm)
|
|
{
|
|
#ifdef FE_UPWARD
|
|
case FE_UPWARD:
|
|
if (dret > ret)
|
|
ret = (res + 0x1p-1022) * c.x;
|
|
break;
|
|
#endif
|
|
|
|
#ifdef FE_DOWNWARD
|
|
case FE_DOWNWARD:
|
|
#endif
|
|
#ifdef FE_TOWARDZERO
|
|
case FE_TOWARDZERO:
|
|
#endif
|
|
#if defined FE_DOWNWARD || defined FE_TOWARDZERO
|
|
if (dret < ret)
|
|
ret = (res - 0x1p-1022) * c.x;
|
|
break;
|
|
#endif
|
|
|
|
default:
|
|
break;
|
|
}
|
|
}
|
|
/* Otherwise (x / ret == ret), either the square root was exact or
|
|
the division was inexact. */
|
|
return ret;
|
|
}
|
|
else
|
|
{
|
|
if ((k & 0x7ff00000) == 0x7ff00000)
|
|
return x * x + x; /* sqrt(NaN)=NaN, sqrt(+inf)=+inf, sqrt(-inf)=sNaN */
|
|
if (x == 0)
|
|
return x; /* sqrt(+0)=+0, sqrt(-0)=-0 */
|
|
if (k < 0)
|
|
return (x - x) / (x - x); /* sqrt(-ve)=sNaN */
|
|
return 0x1p-256 * __ieee754_sqrt (x * 0x1p512);
|
|
}
|
|
}
|
|
strong_alias (__ieee754_sqrt, __sqrt_finite)
|