1
0
mirror of https://sourceware.org/git/glibc.git synced 2025-12-06 12:01:08 +03:00

math: Remove the SVID error handling wrapper from sqrt

i386 and m68k architectures should use math-use-builtins-sqrt.h rather
than relying on architecture-specific or inline assembly implementations.

The PowerPC optimization for PPC 601/603 (30 years old) is removed.

Tested on x86_64-linux-gnu and i686-linux-gnu.

Reviewed-by: Wilco Dijkstra  <Wilco.Dijkstra@arm.com>
This commit is contained in:
Adhemerval Zanella
2025-10-31 13:08:56 -03:00
parent f27a146409
commit 0dfc849eff
36 changed files with 123 additions and 210 deletions

View File

@@ -164,7 +164,7 @@ fabsf128 (_Float128 x)
# define MATH_REDIRECT_UNARY_ARGS(TYPE) TYPE
# define MATH_REDIRECT_BINARY_ARGS(TYPE) TYPE, TYPE
# define MATH_REDIRECT_TERNARY_ARGS(TYPE) TYPE, TYPE, TYPE
MATH_REDIRECT (sqrt, "__ieee754_", MATH_REDIRECT_UNARY_ARGS)
MATH_REDIRECT (sqrt, "__", MATH_REDIRECT_UNARY_ARGS)
MATH_REDIRECT (ceil, "__", MATH_REDIRECT_UNARY_ARGS)
MATH_REDIRECT (floor, "__", MATH_REDIRECT_UNARY_ARGS)
MATH_REDIRECT (roundeven, "__", MATH_REDIRECT_UNARY_ARGS)

View File

@@ -700,6 +700,7 @@ libm {
remainder;
remainderf;
sinhf;
sqrtf;
y0f;
y1f;
ynf;

View File

@@ -22,15 +22,15 @@
#include <libm-alias-float.h>
#if LIBM_SVID_COMPAT
#if LIBM_SVID_COMPAT && SHLIB_COMPAT (libm, GLIBC_2_0, GLIBC_2_43)
/* wrapper sqrtf */
float
__sqrtf (float x)
__sqrtf_svid (float x)
{
if (__builtin_expect (isless (x, 0.0f), 0) && _LIB_VERSION != _IEEE_)
return __kernel_standard_f (x, x, 126); /* sqrt(negative) */
return __ieee754_sqrtf (x);
return __sqrtf (x);
}
libm_alias_float (__sqrt, sqrt)
compat_symbol (libm, __sqrtf_svid, sqrtf, GLIBC_2_0);
#endif

View File

@@ -3,6 +3,6 @@
#include <sysdeps/ieee754/flt-32/e_sqrtf.c>
#if SHLIB_COMPAT (libm, GLIBC_2_18, GLIBC_2_31)
strong_alias (__ieee754_sqrtf, __sqrtf_finite_2_18)
strong_alias (__sqrtf, __sqrtf_finite_2_18)
compat_symbol (libm, __sqrtf_finite_2_18, __sqrtf_finite, GLIBC_2_18);
#endif

View File

@@ -1,13 +0,0 @@
/*
* Public domain.
*/
#include <machine/asm.h>
#include <libm-alias-finite.h>
ENTRY(__ieee754_sqrtf)
flds 4(%esp)
fsqrt
ret
END (__ieee754_sqrtf)
libm_alias_finite (__ieee754_sqrtf, __sqrtf)

View File

@@ -0,0 +1,4 @@
#define USE_SQRT_BUILTIN 1
#define USE_SQRTF_BUILTIN 1
#define USE_SQRTL_BUILTIN 0
#define USE_SQRTF128_BUILTIN 0

View File

@@ -12,81 +12,103 @@
* ====================================================
*/
/* The internal alias to avoid PLT calls interfere with the default
symbol alias for !LIBM_SVID_COMPAT. */
#define sqrtf __redirect_sqrtf
#include <math.h>
#include <math_private.h>
#undef sqrtf
#include <libm-alias-finite.h>
#include <libm-alias-float.h>
#include <math-svid-compat.h>
#include <math-use-builtins.h>
#include "math_config.h"
float
__ieee754_sqrtf(float x)
__sqrtf (float x)
{
#if USE_SQRTF_BUILTIN
return __builtin_sqrtf (x);
if (__glibc_unlikely (isless (x, 0.0f)))
return __math_invalidf (x);
return __builtin_sqrtf (x);
#else
/* Use generic implementation. */
float z;
int32_t sign = (int)0x80000000;
int32_t ix,s,q,m,t,i;
uint32_t r;
/* Use generic implementation. */
float z;
int32_t sign = (int) 0x80000000;
int32_t ix, s, q, m, t, i;
uint32_t r;
GET_FLOAT_WORD(ix,x);
ix = asuint (x);
/* take care of Inf and NaN */
if((ix&0x7f800000)==0x7f800000) {
return x*x+x; /* sqrt(NaN)=NaN, sqrt(+inf)=+inf
sqrt(-inf)=sNaN */
}
/* take care of zero */
if(ix<=0) {
if((ix&(~sign))==0) return x;/* sqrt(+-0) = +-0 */
else if(ix<0)
return (x-x)/(x-x); /* sqrt(-ve) = sNaN */
}
/* normalize x */
m = (ix>>23);
if(m==0) { /* subnormal x */
for(i=0;(ix&0x00800000)==0;i++) ix<<=1;
m -= i-1;
}
m -= 127; /* unbias exponent */
ix = (ix&0x007fffff)|0x00800000;
if(m&1) /* odd m, double x to make it even */
ix += ix;
m >>= 1; /* m = [m/2] */
/* take care of Inf and NaN */
if ((ix & 0x7f800000) == 0x7f800000)
{
if (ix == 0xff800000)
return __math_invalidf (0.0f);
return x * x + x; /* sqrt(NaN)=NaN, sqrt(+inf)=+inf
sqrt(-inf)=sNaN */
}
/* take care of zero */
if (ix <= 0)
{
if ((ix & (~sign)) == 0)
return x; /* sqrt(+-0) = +-0 */
else if (ix < 0)
return __math_invalidf (0.0f); /* sqrt(-ve) = sNaN */
}
/* normalize x */
m = (ix >> 23);
if (m == 0)
{ /* subnormal x */
for (i = 0; (ix & 0x00800000) == 0; i++)
ix <<= 1;
m -= i - 1;
}
m -= 127; /* unbias exponent */
ix = (ix & 0x007fffff) | 0x00800000;
if (m & 1) /* odd m, double x to make it even */
ix += ix;
m >>= 1; /* m = [m/2] */
/* generate sqrt(x) bit by bit */
ix += ix;
q = s = 0; /* q = sqrt(x) */
r = 0x01000000; /* r = moving bit from right to left */
/* generate sqrt(x) bit by bit */
ix += ix;
q = s = 0; /* q = sqrt(x) */
r = 0x01000000; /* r = moving bit from right to left */
while(r!=0) {
t = s+r;
if(t<=ix) {
s = t+r;
ix -= t;
q += r;
}
ix += ix;
r>>=1;
while (r != 0)
{
t = s + r;
if (t <= ix)
{
s = t + r;
ix -= t;
q += r;
}
ix += ix;
r >>= 1;
}
/* use floating add to find out rounding direction */
if(ix!=0) {
z = 0x1p0f - 0x1.4484cp-100f; /* trigger inexact flag. */
if (z >= 0x1p0f) { /* rounding to nearest or upward */
z = 0x1p0f + 0x1.4484cp-100f;
if (z > 0x1p0f) /* rounding upward */
q += 2;
else
q += (q&1);
}
/* use floating add to find out rounding direction */
if (ix != 0)
{
z = 0x1p0 - 0x1.4484cp-100; /* trigger inexact flag. */
if (z >= 0x1p0)
{
z = 0x1p0 + 0x1.4484cp-100;
if (z > 0x1p0)
q += 2;
else
q += (q & 1);
}
ix = (q>>1)+0x3f000000;
ix += (m <<23);
SET_FLOAT_WORD(z,ix);
return z;
}
ix = (q >> 1) + 0x3f000000;
ix += (m << 23);
return asfloat (ix);
#endif /* ! USE_SQRTF_BUILTIN */
}
#ifndef __ieee754_sqrtf
libm_alias_finite (__ieee754_sqrtf, __sqrtf)
libm_alias_finite (__sqrtf, __sqrtf)
#if LIBM_SVID_COMPAT
versioned_symbol (libm, __sqrtf, sqrtf, GLIBC_2_43);
libm_alias_float_other (__sqrt, sqrt)
#else
libm_alias_float (__sqrt, sqrt)
#endif

View File

@@ -0,0 +1 @@
/* Not needed */

View File

@@ -0,0 +1,4 @@
#define USE_SQRT_BUILTIN 1
#define USE_SQRTF_BUILTIN 1
#define USE_SQRTL_BUILTIN 0
#define USE_SQRTF128_BUILTIN 0

View File

@@ -1,3 +0,0 @@
#define FUNC __ieee754_sqrtf
#define FUNC_FINITE __sqrtf
#include <e_acosf.c>

View File

@@ -1331,6 +1331,7 @@ GLIBC_2.43 log10f F
GLIBC_2.43 remainder F
GLIBC_2.43 remainderf F
GLIBC_2.43 sinhf F
GLIBC_2.43 sqrtf F
GLIBC_2.43 y0f F
GLIBC_2.43 y1f F
GLIBC_2.43 ynf F

View File

@@ -1,128 +0,0 @@
/* Single-precision floating point square root.
Copyright (C) 1997-2025 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
<https://www.gnu.org/licenses/>. */
#include <math.h>
#include <math_private.h>
#include <fenv_libc.h>
#include <libm-alias-finite.h>
#include <math-use-builtins.h>
float
__ieee754_sqrtf (float x)
{
#if USE_SQRTF_BUILTIN
return __builtin_sqrtf (x);
#else
/* The method is based on a description in
Computation of elementary functions on the IBM RISC System/6000 processor,
P. W. Markstein, IBM J. Res. Develop, 34(1) 1990.
Basically, it consists of two interleaved Newton-Raphson approximations,
one to find the actual square root, and one to find its reciprocal
without the expense of a division operation. The tricky bit here
is the use of the POWER/PowerPC multiply-add operation to get the
required accuracy with high speed.
The argument reduction works by a combination of table lookup to
obtain the initial guesses, and some careful modification of the
generated guesses (which mostly runs on the integer unit, while the
Newton-Raphson is running on the FPU). */
extern const float __t_sqrt[1024];
if (x > 0)
{
if (x != INFINITY)
{
/* Variables named starting with 's' exist in the
argument-reduced space, so that 2 > sx >= 0.5,
1.41... > sg >= 0.70.., 0.70.. >= sy > 0.35... .
Variables named ending with 'i' are integer versions of
floating-point values. */
float sx; /* The value of which we're trying to find the square
root. */
float sg, g; /* Guess of the square root of x. */
float sd, d; /* Difference between the square of the guess and x. */
float sy; /* Estimate of 1/2g (overestimated by 1ulp). */
float sy2; /* 2*sy */
float e; /* Difference between y*g and 1/2 (note that e==se). */
float shx; /* == sx * fsg */
float fsg; /* sg*fsg == g. */
fenv_t fe; /* Saved floating-point environment (stores rounding
mode and whether the inexact exception is
enabled). */
uint32_t xi, sxi, fsgi;
const float *t_sqrt;
GET_FLOAT_WORD (xi, x);
fe = fegetenv_register ();
relax_fenv_state ();
sxi = (xi & 0x3fffffff) | 0x3f000000;
SET_FLOAT_WORD (sx, sxi);
t_sqrt = __t_sqrt + (xi >> (23 - 8 - 1) & 0x3fe);
sg = t_sqrt[0];
sy = t_sqrt[1];
/* Here we have three Newton-Raphson iterations each of a
division and a square root and the remainder of the
argument reduction, all interleaved. */
sd = -__builtin_fmaf (sg, sg, -sx);
fsgi = (xi + 0x40000000) >> 1 & 0x7f800000;
sy2 = sy + sy;
sg = __builtin_fmaf (sy, sd, sg); /* 16-bit approximation to
sqrt(sx). */
e = -__builtin_fmaf (sy, sg, -0x1.0000020365653p-1);
SET_FLOAT_WORD (fsg, fsgi);
sd = -__builtin_fmaf (sg, sg, -sx);
sy = __builtin_fmaf (e, sy2, sy);
if ((xi & 0x7f800000) == 0)
goto denorm;
shx = sx * fsg;
sg = __builtin_fmaf (sy, sd, sg); /* 32-bit approximation to
sqrt(sx), but perhaps
rounded incorrectly. */
sy2 = sy + sy;
g = sg * fsg;
e = -__builtin_fmaf (sy, sg, -0x1.0000020365653p-1);
d = -__builtin_fmaf (g, sg, -shx);
sy = __builtin_fmaf (e, sy2, sy);
fesetenv_register (fe);
return __builtin_fmaf (sy, d, g);
denorm:
/* For denormalised numbers, we normalise, calculate the
square root, and return an adjusted result. */
fesetenv_register (fe);
return __ieee754_sqrtf (x * 0x1p+48) * 0x1p-24;
}
}
else if (x < 0)
{
/* For some reason, some PowerPC32 processors don't implement
FE_INVALID_SQRT. */
# ifdef FE_INVALID_SQRT
feraiseexcept (FE_INVALID_SQRT);
fenv_union_t u = { .fenv = fegetenv_register () };
if ((u.l & FE_INVALID) == 0)
# endif
feraiseexcept (FE_INVALID);
x = NAN;
}
return f_washf (x);
#endif /* USE_SQRTF_BUILTIN */
}
libm_alias_finite (__ieee754_sqrtf, __sqrtf)

View File

@@ -1297,6 +1297,7 @@ GLIBC_2.43 log10f F
GLIBC_2.43 remainder F
GLIBC_2.43 remainderf F
GLIBC_2.43 sinhf F
GLIBC_2.43 sqrtf F
GLIBC_2.43 y0f F
GLIBC_2.43 y1f F
GLIBC_2.43 ynf F

View File

@@ -1456,6 +1456,7 @@ GLIBC_2.43 log10f F
GLIBC_2.43 remainder F
GLIBC_2.43 remainderf F
GLIBC_2.43 sinhf F
GLIBC_2.43 sqrtf F
GLIBC_2.43 y0f F
GLIBC_2.43 y1f F
GLIBC_2.43 ynf F

View File

@@ -962,6 +962,7 @@ GLIBC_2.43 log10f F
GLIBC_2.43 remainder F
GLIBC_2.43 remainderf F
GLIBC_2.43 sinhf F
GLIBC_2.43 sqrtf F
GLIBC_2.43 y0f F
GLIBC_2.43 y1f F
GLIBC_2.43 ynf F

View File

@@ -962,6 +962,7 @@ GLIBC_2.43 log10f F
GLIBC_2.43 remainder F
GLIBC_2.43 remainderf F
GLIBC_2.43 sinhf F
GLIBC_2.43 sqrtf F
GLIBC_2.43 y0f F
GLIBC_2.43 y1f F
GLIBC_2.43 ynf F

View File

@@ -962,6 +962,7 @@ GLIBC_2.43 log10f F
GLIBC_2.43 remainder F
GLIBC_2.43 remainderf F
GLIBC_2.43 sinhf F
GLIBC_2.43 sqrtf F
GLIBC_2.43 y0f F
GLIBC_2.43 y1f F
GLIBC_2.43 ynf F

View File

@@ -1338,6 +1338,7 @@ GLIBC_2.43 log10f F
GLIBC_2.43 remainder F
GLIBC_2.43 remainderf F
GLIBC_2.43 sinhf F
GLIBC_2.43 sqrtf F
GLIBC_2.43 y0f F
GLIBC_2.43 y1f F
GLIBC_2.43 ynf F

View File

@@ -962,6 +962,7 @@ GLIBC_2.43 log10f F
GLIBC_2.43 remainder F
GLIBC_2.43 remainderf F
GLIBC_2.43 sinhf F
GLIBC_2.43 sqrtf F
GLIBC_2.43 y0f F
GLIBC_2.43 y1f F
GLIBC_2.43 ynf F

View File

@@ -994,6 +994,7 @@ GLIBC_2.43 j1f F
GLIBC_2.43 jnf F
GLIBC_2.43 remainder F
GLIBC_2.43 remainderf F
GLIBC_2.43 sqrtf F
GLIBC_2.43 y0f F
GLIBC_2.43 y1f F
GLIBC_2.43 ynf F

View File

@@ -962,6 +962,7 @@ GLIBC_2.43 log10f F
GLIBC_2.43 remainder F
GLIBC_2.43 remainderf F
GLIBC_2.43 sinhf F
GLIBC_2.43 sqrtf F
GLIBC_2.43 y0f F
GLIBC_2.43 y1f F
GLIBC_2.43 ynf F

View File

@@ -962,6 +962,7 @@ GLIBC_2.43 log10f F
GLIBC_2.43 remainder F
GLIBC_2.43 remainderf F
GLIBC_2.43 sinhf F
GLIBC_2.43 sqrtf F
GLIBC_2.43 y0f F
GLIBC_2.43 y1f F
GLIBC_2.43 ynf F

View File

@@ -962,6 +962,7 @@ GLIBC_2.43 log10f F
GLIBC_2.43 remainder F
GLIBC_2.43 remainderf F
GLIBC_2.43 sinhf F
GLIBC_2.43 sqrtf F
GLIBC_2.43 y0f F
GLIBC_2.43 y1f F
GLIBC_2.43 ynf F

View File

@@ -1297,6 +1297,7 @@ GLIBC_2.43 log10f F
GLIBC_2.43 remainder F
GLIBC_2.43 remainderf F
GLIBC_2.43 sinhf F
GLIBC_2.43 sqrtf F
GLIBC_2.43 y0f F
GLIBC_2.43 y1f F
GLIBC_2.43 ynf F

View File

@@ -1109,6 +1109,7 @@ GLIBC_2.43 log10f F
GLIBC_2.43 remainder F
GLIBC_2.43 remainderf F
GLIBC_2.43 sinhf F
GLIBC_2.43 sqrtf F
GLIBC_2.43 y0f F
GLIBC_2.43 y1f F
GLIBC_2.43 ynf F

View File

@@ -1108,6 +1108,7 @@ GLIBC_2.43 log10f F
GLIBC_2.43 remainder F
GLIBC_2.43 remainderf F
GLIBC_2.43 sinhf F
GLIBC_2.43 sqrtf F
GLIBC_2.43 y0f F
GLIBC_2.43 y1f F
GLIBC_2.43 ynf F

View File

@@ -1102,6 +1102,7 @@ GLIBC_2.43 log10f F
GLIBC_2.43 remainder F
GLIBC_2.43 remainderf F
GLIBC_2.43 sinhf F
GLIBC_2.43 sqrtf F
GLIBC_2.43 y0f F
GLIBC_2.43 y1f F
GLIBC_2.43 ynf F

View File

@@ -1486,6 +1486,7 @@ GLIBC_2.43 log10f F
GLIBC_2.43 remainder F
GLIBC_2.43 remainderf F
GLIBC_2.43 sinhf F
GLIBC_2.43 sqrtf F
GLIBC_2.43 y0f F
GLIBC_2.43 y1f F
GLIBC_2.43 ynf F

View File

@@ -1400,6 +1400,7 @@ GLIBC_2.43 log10f F
GLIBC_2.43 remainder F
GLIBC_2.43 remainderf F
GLIBC_2.43 sinhf F
GLIBC_2.43 sqrtf F
GLIBC_2.43 y0f F
GLIBC_2.43 y1f F
GLIBC_2.43 ynf F

View File

@@ -1400,6 +1400,7 @@ GLIBC_2.43 log10f F
GLIBC_2.43 remainder F
GLIBC_2.43 remainderf F
GLIBC_2.43 sinhf F
GLIBC_2.43 sqrtf F
GLIBC_2.43 y0f F
GLIBC_2.43 y1f F
GLIBC_2.43 ynf F

View File

@@ -962,6 +962,7 @@ GLIBC_2.43 log10f F
GLIBC_2.43 remainder F
GLIBC_2.43 remainderf F
GLIBC_2.43 sinhf F
GLIBC_2.43 sqrtf F
GLIBC_2.43 y0f F
GLIBC_2.43 y1f F
GLIBC_2.43 ynf F

View File

@@ -962,6 +962,7 @@ GLIBC_2.43 log10f F
GLIBC_2.43 remainder F
GLIBC_2.43 remainderf F
GLIBC_2.43 sinhf F
GLIBC_2.43 sqrtf F
GLIBC_2.43 y0f F
GLIBC_2.43 y1f F
GLIBC_2.43 ynf F

View File

@@ -1407,6 +1407,7 @@ GLIBC_2.43 log10f F
GLIBC_2.43 remainder F
GLIBC_2.43 remainderf F
GLIBC_2.43 sinhf F
GLIBC_2.43 sqrtf F
GLIBC_2.43 y0f F
GLIBC_2.43 y1f F
GLIBC_2.43 ynf F

View File

@@ -1297,6 +1297,7 @@ GLIBC_2.43 log10f F
GLIBC_2.43 remainder F
GLIBC_2.43 remainderf F
GLIBC_2.43 sinhf F
GLIBC_2.43 sqrtf F
GLIBC_2.43 y0f F
GLIBC_2.43 y1f F
GLIBC_2.43 ynf F

View File

@@ -1330,6 +1330,7 @@ GLIBC_2.43 log10f F
GLIBC_2.43 remainder F
GLIBC_2.43 remainderf F
GLIBC_2.43 sinhf F
GLIBC_2.43 sqrtf F
GLIBC_2.43 y0f F
GLIBC_2.43 y1f F
GLIBC_2.43 ynf F

View File

@@ -1330,6 +1330,7 @@ GLIBC_2.43 log10f F
GLIBC_2.43 remainder F
GLIBC_2.43 remainderf F
GLIBC_2.43 sinhf F
GLIBC_2.43 sqrtf F
GLIBC_2.43 y0f F
GLIBC_2.43 y1f F
GLIBC_2.43 ynf F