diff options
Diffstat (limited to 'sysdeps/ieee754')
27 files changed, 446 insertions, 229 deletions
diff --git a/sysdeps/ieee754/dbl-64/e_ilogb.c b/sysdeps/ieee754/dbl-64/e_ilogb.c index 1e338a5..1ea2f23 100644 --- a/sysdeps/ieee754/dbl-64/e_ilogb.c +++ b/sysdeps/ieee754/dbl-64/e_ilogb.c @@ -1,63 +1 @@ -/* @(#)s_ilogb.c 5.1 93/09/24 */ -/* - * ==================================================== - * 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. - * ==================================================== - */ - -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: s_ilogb.c,v 1.9 1995/05/10 20:47:28 jtc Exp $"; -#endif - -/* ilogb(double x) - * return the binary exponent of non-zero x - * ilogb(0) = FP_ILOGB0 - * ilogb(NaN) = FP_ILOGBNAN (no signal is raised) - * ilogb(+-Inf) = INT_MAX (no signal is raised) - */ - -#include <limits.h> -#include <math.h> -#include <math_private.h> - -int -__ieee754_ilogb (double x) -{ - int32_t hx, lx, ix; - - GET_HIGH_WORD (hx, x); - hx &= 0x7fffffff; - if (hx < 0x00100000) - { - GET_LOW_WORD (lx, x); - if ((hx | lx) == 0) - return FP_ILOGB0; /* ilogb(0) = FP_ILOGB0 */ - else /* subnormal x */ - if (hx == 0) - { - for (ix = -1043; lx > 0; lx <<= 1) - ix -= 1; - } - else - { - for (ix = -1022, hx <<= 11; hx > 0; hx <<= 1) - ix -= 1; - } - return ix; - } - else if (hx < 0x7ff00000) - return (hx >> 20) - 1023; - else if (FP_ILOGBNAN != INT_MAX) - { - /* ISO C99 requires ilogb(+-Inf) == INT_MAX. */ - GET_LOW_WORD (lx, x); - if (((hx ^ 0x7ff00000) | lx) == 0) - return INT_MAX; - } - return FP_ILOGBNAN; -} +/* ilogb is implemented at w_ilogb.c */ diff --git a/sysdeps/ieee754/dbl-64/math_config.h b/sysdeps/ieee754/dbl-64/math_config.h index 3382e38..d9288c4 100644 --- a/sysdeps/ieee754/dbl-64/math_config.h +++ b/sysdeps/ieee754/dbl-64/math_config.h @@ -109,6 +109,7 @@ issignaling_inline (double x) #define BIT_WIDTH 64 #define MANTISSA_WIDTH 52 #define EXPONENT_WIDTH 11 +#define EXPONENT_BIAS 1023 #define MANTISSA_MASK UINT64_C(0x000fffffffffffff) #define EXPONENT_MASK UINT64_C(0x7ff0000000000000) #define EXP_MANT_MASK UINT64_C(0x7fffffffffffffff) @@ -121,12 +122,24 @@ is_nan (uint64_t x) return (x & EXP_MANT_MASK) > EXPONENT_MASK; } +static inline bool +is_inf (uint64_t x) +{ + return (x << 1) == (EXPONENT_MASK << 1); +} + static inline uint64_t get_mantissa (uint64_t x) { return x & MANTISSA_MASK; } +static inline int +get_exponent (uint64_t x) +{ + return (int)((x >> MANTISSA_WIDTH & 0x7ff) - EXPONENT_BIAS); +} + /* Convert integer number X, unbiased exponent EP, and sign S to double: result = X * 2^(EP+1 - exponent_bias) @@ -164,6 +177,8 @@ attribute_hidden double __math_divzero (uint32_t); /* Invalid input unless it is a quiet NaN. */ attribute_hidden double __math_invalid (double); +attribute_hidden int __math_invalid_i (int); +attribute_hidden long int __math_invalid_li (long int); /* Error handling using output checking, only for errno setting. */ diff --git a/sysdeps/ieee754/dbl-64/math_err.c b/sysdeps/ieee754/dbl-64/math_err.c index 4a07fd5..b8c645a 100644 --- a/sysdeps/ieee754/dbl-64/math_err.c +++ b/sysdeps/ieee754/dbl-64/math_err.c @@ -29,8 +29,24 @@ with_errno (double y, int e) errno = e; return y; } + +NOINLINE static int +with_errno_i (int y, int e) +{ + errno = e; + return y; +} + +NOINLINE static long int +with_errno_li (long int y, int e) +{ + errno = e; + return y; +} #else #define with_errno(x, e) (x) +#define with_errno_i(x, e) (x) +#define with_errno_li(x, e) (x) #endif attribute_hidden double @@ -83,6 +99,22 @@ __math_invalid (double x) return isnan (x) ? y : with_errno (y, EDOM); } +attribute_hidden int +__math_invalid_i (int r) +{ + double y = 0.0 / 0.0; + math_force_eval (y); + return with_errno_i (r, EDOM); +} + +attribute_hidden long int +__math_invalid_li (long int r) +{ + double y = 0.0 / 0.0; + math_force_eval (y); + return with_errno_li (r, EDOM); +} + /* Check result and set errno if necessary. */ attribute_hidden double diff --git a/sysdeps/ieee754/dbl-64/s_modf.c b/sysdeps/ieee754/dbl-64/s_modf.c index 0de2084..90cd8e8 100644 --- a/sysdeps/ieee754/dbl-64/s_modf.c +++ b/sysdeps/ieee754/dbl-64/s_modf.c @@ -1,63 +1,68 @@ -/* - * ==================================================== - * 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. - * ==================================================== - */ +/* Extract signed integral and fractional values. + Copyright (C) 1993-2025 Free Software Foundation, Inc. + This file is part of the GNU C Library. -/* - * modf(double x, double *iptr) - * return fraction part of x, and return x's integral part in *iptr. - * Method: - * Bit twiddling. - * - * Exception: - * No exception. - */ + 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 <libm-alias-double.h> -#include <stdint.h> - -static const double one = 1.0; +#include "math_config.h" +#include <math-use-builtins-trunc.h> double -__modf(double x, double *iptr) +__modf (double x, double *iptr) { - int64_t i0; - int32_t j0; - EXTRACT_WORDS64(i0,x); - j0 = ((i0>>52)&0x7ff)-0x3ff; /* exponent of x */ - if(j0<52) { /* integer part in x */ - if(j0<0) { /* |x|<1 */ - /* *iptr = +-0 */ - INSERT_WORDS64(*iptr,i0&UINT64_C(0x8000000000000000)); - return x; - } else { - uint64_t i = UINT64_C(0x000fffffffffffff)>>j0; - if((i0&i)==0) { /* x is integral */ - *iptr = x; - /* return +-0 */ - INSERT_WORDS64(x,i0&UINT64_C(0x8000000000000000)); - return x; - } else { - INSERT_WORDS64(*iptr,i0&(~i)); - return x - *iptr; - } - } - } else { /* no fraction part */ - *iptr = x*one; - /* We must handle NaNs separately. */ - if (j0 == 0x400 && (i0 & UINT64_C(0xfffffffffffff))) - return x*one; - INSERT_WORDS64(x,i0&UINT64_C(0x8000000000000000)); /* return +-0 */ - return x; + uint64_t t = asuint64 (x); +#if USE_TRUNC_BUILTIN + if (is_inf (t)) + { + *iptr = x; + return copysign (0.0, x); + } + *iptr = trunc (x); + return copysign (x - *iptr, x); +#else + int e = get_exponent (t); + /* No fraction part. */ + if (e < MANTISSA_WIDTH) + { + if (e < 0) + { + /* |x|<1 -> *iptr = +-0 */ + *iptr = asdouble (t & SIGN_MASK); + return x; + } + + uint64_t i = MANTISSA_MASK >> e; + if ((t & i) == 0) + { + /* x in integral, return +-0 */ + *iptr = x; + return asdouble (t & SIGN_MASK); } + + *iptr = asdouble (t & ~i); + return x - *iptr; + } + + /* Set invalid operation for sNaN. */ + *iptr = x * 1.0; + if ((e == 0x400) && (t & MANTISSA_MASK)) + return *iptr; + return asdouble (t & SIGN_MASK); +#endif } #ifndef __modf libm_alias_double (__modf, modf) diff --git a/sysdeps/ieee754/dbl-64/w_ilogb-impl.h b/sysdeps/ieee754/dbl-64/w_ilogb-impl.h new file mode 100644 index 0000000..c919735 --- /dev/null +++ b/sysdeps/ieee754/dbl-64/w_ilogb-impl.h @@ -0,0 +1,37 @@ +/* Get integer exponent of a floating-point value. + Copyright (C) 1999-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/>. */ + +static inline RET_TYPE +IMPL_NAME (double x) +{ + uint64_t ux = asuint64 (x); + int ex = (ux & ~SIGN_MASK) >> MANTISSA_WIDTH; + if (__glibc_unlikely (ex == 0)) /* zero or subnormal */ + { + /* Clear sign and exponent */ + ux <<= 12; + if (ux == 0) + return RET_INVALID (RET_LOGB0); + /* subnormal */ + return (RET_TYPE)-1023 - stdc_leading_zeros (ux); + } + if (__glibc_unlikely (ex == EXPONENT_MASK >> MANTISSA_WIDTH)) + /* NaN or Inf */ + return RET_INVALID (ux << 12 ? RET_LOGBNAN : RET_LOGMAX); + return ex - 1023; +} diff --git a/sysdeps/ieee754/dbl-64/w_ilogb.c b/sysdeps/ieee754/dbl-64/w_ilogb.c new file mode 100644 index 0000000..e460f14 --- /dev/null +++ b/sysdeps/ieee754/dbl-64/w_ilogb.c @@ -0,0 +1,52 @@ +/* Get integer exponent of a floating-point value. + Copyright (C) 1999-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 <limits.h> +#include <math.h> +#include <stdbit.h> +#include <libm-alias-double.h> +#include "math_config.h" + +#ifdef DEF_AS_LLOGB +# define DECL_NAME __llogb +# define FUNC_NAME llogb +# define RET_TYPE long int +# define RET_LOGB0 FP_LLOGB0 +# define RET_LOGBNAN FP_LLOGBNAN +# define RET_LOGMAX LONG_MAX +# define RET_INVALID __math_invalid_li +#else +# define DECL_NAME __ilogb +# define FUNC_NAME ilogb +# define RET_TYPE int +# define RET_LOGB0 FP_ILOGB0 +# define RET_LOGBNAN FP_ILOGBNAN +# define RET_LOGMAX INT_MAX +# define RET_INVALID __math_invalid_i +#endif +#define __IMPL_NAME(x,y) x ## _ ## y +#define _IMPL_NAME(x,y) __IMPL_NAME(x,y) +#define IMPL_NAME _IMPL_NAME(FUNC_NAME, impl) +#include <w_ilogb-impl.h> + +RET_TYPE +DECL_NAME (double x) +{ + return IMPL_NAME (x); +} +libm_alias_double (DECL_NAME, FUNC_NAME) diff --git a/sysdeps/ieee754/dbl-64/w_llogb.c b/sysdeps/ieee754/dbl-64/w_llogb.c new file mode 100644 index 0000000..c984cd15 --- /dev/null +++ b/sysdeps/ieee754/dbl-64/w_llogb.c @@ -0,0 +1,2 @@ +#define DEF_AS_LLOGB +#include "w_ilogb.c" diff --git a/sysdeps/ieee754/flt-32/e_atanhf.c b/sysdeps/ieee754/flt-32/e_atanhf.c index 5138408..eeb1aae 100644 --- a/sysdeps/ieee754/flt-32/e_atanhf.c +++ b/sysdeps/ieee754/flt-32/e_atanhf.c @@ -3,7 +3,7 @@ Copyright (c) 2023-2024 Alexei Sibidanov. The original version of this file was copied from the CORE-MATH -project (file src/binary32/acosh/acoshf.c, revision bc385c2). +project (file src/binary32/acosh/acoshf.c, revision 4d6192d2). Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal @@ -135,7 +135,7 @@ __ieee754_atanhf (float x) } double sgn = s[ux >> 31]; unsigned int e = ax >> 24; - unsigned int md = ((ux << 8) | 1 << 31) >> (126 - e); + unsigned int md = ((ux << 8) | 1U << 31) >> (126 - e); unsigned int mn = -md; int nz = __builtin_clz (mn) + 1; mn <<= nz; diff --git a/sysdeps/ieee754/flt-32/e_coshf.c b/sysdeps/ieee754/flt-32/e_coshf.c index 5f6ff8c..382cd55 100644 --- a/sysdeps/ieee754/flt-32/e_coshf.c +++ b/sysdeps/ieee754/flt-32/e_coshf.c @@ -3,7 +3,7 @@ Copyright (c) 2022-2024 Alexei Sibidanov. The original version of this file was copied from the CORE-MATH -project (file src/binary32/cosh/coshf.c, revision 5c58ea1). +project (file src/binary32/cosh/coshf.c, revision de59ecfb). Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal @@ -69,8 +69,8 @@ __ieee754_coshf (float x) double h2 = h * h; int64_t jp = asuint64 (ia + 0x1.8p52); int64_t jm = -jp; - double sp = asdouble (TB[jp & 31] + ((jp >> 5) << 52)); - double sm = asdouble (TB[jm & 31] + ((jm >> 5) << 52)); + double sp = asdouble (TB[jp & 31] + ((uint64_t)(jp >> 5) << 52)); + double sm = asdouble (TB[jm & 31] + ((uint64_t)(jm >> 5) << 52)); double te = C[0] + h2 * C[2]; double to = (C[1] + h2 * C[3]); double rp = sp * (te + h * to); diff --git a/sysdeps/ieee754/flt-32/e_ilogbf.c b/sysdeps/ieee754/flt-32/e_ilogbf.c index db24012..a27fb94 100644 --- a/sysdeps/ieee754/flt-32/e_ilogbf.c +++ b/sysdeps/ieee754/flt-32/e_ilogbf.c @@ -1,43 +1 @@ -/* s_ilogbf.c -- float version of s_ilogb.c. - */ - -/* - * ==================================================== - * 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. - * ==================================================== - */ - -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: s_ilogbf.c,v 1.4 1995/05/10 20:47:31 jtc Exp $"; -#endif - -#include <limits.h> -#include <math.h> -#include <math_private.h> - -int __ieee754_ilogbf(float x) -{ - int32_t hx,ix; - - GET_FLOAT_WORD(hx,x); - hx &= 0x7fffffff; - if(hx<0x00800000) { - if(hx==0) - return FP_ILOGB0; /* ilogb(0) = FP_ILOGB0 */ - else /* subnormal x */ - for (ix = -126,hx<<=8; hx>0; hx<<=1) ix -=1; - return ix; - } - else if (hx<0x7f800000) return (hx>>23)-127; - else if (FP_ILOGBNAN != INT_MAX) { - /* ISO C99 requires ilogbf(+-Inf) == INT_MAX. */ - if (hx==0x7f800000) - return INT_MAX; - } - return FP_ILOGBNAN; -} +/* ilogbf is implemented at w_ilogbf.c */ diff --git a/sysdeps/ieee754/flt-32/e_logf.c b/sysdeps/ieee754/flt-32/e_logf.c index 6a595cf..207151c 100644 --- a/sysdeps/ieee754/flt-32/e_logf.c +++ b/sysdeps/ieee754/flt-32/e_logf.c @@ -70,7 +70,7 @@ __logf (float x) tmp = ix - OFF; i = (tmp >> (23 - LOGF_TABLE_BITS)) % N; k = (int32_t) tmp >> 23; /* arithmetic shift */ - iz = ix - (tmp & 0x1ff << 23); + iz = ix - (tmp & 0xff800000); invc = T[i].invc; logc = T[i].logc; z = (double_t) asfloat (iz); diff --git a/sysdeps/ieee754/flt-32/e_sinhf.c b/sysdeps/ieee754/flt-32/e_sinhf.c index 754b84a..6c8c1db 100644 --- a/sysdeps/ieee754/flt-32/e_sinhf.c +++ b/sysdeps/ieee754/flt-32/e_sinhf.c @@ -3,7 +3,7 @@ Copyright (c) 2022-2024 Alexei Sibidanov. The original version of this file was copied from the CORE-MATH -project (file src/binary32/sinh/sinhf.c, revision 572ecec). +project (file src/binary32/sinh/sinhf.c, revision bbfabd99). Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal @@ -83,8 +83,8 @@ __ieee754_sinhf (float x) double h2 = h * h; int64_t jp = asuint64 (ia + 0x1.8p52); int64_t jm = -jp; - double sp = asdouble (TB[jp & 31] + ((jp >> 5) << 52)); - double sm = asdouble (TB[jm & 31] + ((jm >> 5) << 52)); + double sp = asdouble (TB[jp & 31] + ((uint64_t)(jp >> 5) << 52)); + double sm = asdouble (TB[jm & 31] + ((uint64_t)(jm >> 5) << 52)); double te = C[0] + h2 * C[2]; double to = (C[1] + h2 * C[3]); double rp = sp * (te + h * to); diff --git a/sysdeps/ieee754/flt-32/math_config.h b/sysdeps/ieee754/flt-32/math_config.h index 8d9c8ee..33ea631 100644 --- a/sysdeps/ieee754/flt-32/math_config.h +++ b/sysdeps/ieee754/flt-32/math_config.h @@ -165,6 +165,7 @@ issignalingf_inline (float x) #define BIT_WIDTH 32 #define MANTISSA_WIDTH 23 #define EXPONENT_WIDTH 8 +#define EXPONENT_BIAS 127 #define MANTISSA_MASK 0x007fffff #define EXPONENT_MASK 0x7f800000 #define EXP_MANT_MASK 0x7fffffff @@ -177,12 +178,24 @@ is_nan (uint32_t x) return (x & EXP_MANT_MASK) > EXPONENT_MASK; } +static inline bool +is_inf (uint32_t x) +{ + return (x << 1) == (EXPONENT_MASK << 1); +} + static inline uint32_t get_mantissa (uint32_t x) { return x & MANTISSA_MASK; } +static inline int +get_exponent (uint32_t x) +{ + return (int)((x >> MANTISSA_WIDTH & 0xff) - EXPONENT_BIAS); +} + /* Convert integer number X, unbiased exponent EP, and sign S to double: result = X * 2^(EP+1 - exponent_bias) @@ -208,6 +221,8 @@ attribute_hidden float __math_uflowf (uint32_t); attribute_hidden float __math_may_uflowf (uint32_t); attribute_hidden float __math_divzerof (uint32_t); attribute_hidden float __math_invalidf (float); +attribute_hidden int __math_invalidf_i (int); +attribute_hidden long int __math_invalidf_li (long int); attribute_hidden float __math_edomf (float x); /* Shared between expf, exp2f, exp10f, and powf. */ diff --git a/sysdeps/ieee754/flt-32/math_errf.c b/sysdeps/ieee754/flt-32/math_errf.c index edcc4c0..244e38a 100644 --- a/sysdeps/ieee754/flt-32/math_errf.c +++ b/sysdeps/ieee754/flt-32/math_errf.c @@ -16,6 +16,7 @@ License along with the GNU C Library; if not, see <https://www.gnu.org/licenses/>. */ +#include <math-barriers.h> #include "math_config.h" #if WANT_ERRNO @@ -27,8 +28,24 @@ with_errnof (float y, int e) errno = e; return y; } + +NOINLINE static int +with_errnof_i (int y, int e) +{ + errno = e; + return y; +} + +NOINLINE static long int +with_errnof_li (long int y, int e) +{ + errno = e; + return y; +} #else # define with_errnof(x, e) (x) +# define with_errnof_i(x, x) (x) +# define with_errnof_li(x, x) (x) #endif attribute_hidden float @@ -80,3 +97,19 @@ __math_invalidf (float x) float y = (x - x) / (x - x); return isnan (x) ? y : with_errnof (y, EDOM); } + +attribute_hidden int +__math_invalidf_i (int x) +{ + float y = 0.0f / 0.0f; + math_force_eval (y); + return with_errnof_i (x, EDOM); +} + +attribute_hidden long int +__math_invalidf_li (long int x) +{ + float y = 0.0f / 0.0f; + math_force_eval (y); + return with_errnof_li (x, EDOM); +} diff --git a/sysdeps/ieee754/flt-32/s_cbrtf.c b/sysdeps/ieee754/flt-32/s_cbrtf.c index 5a7a9a9..df9e888 100644 --- a/sysdeps/ieee754/flt-32/s_cbrtf.c +++ b/sysdeps/ieee754/flt-32/s_cbrtf.c @@ -3,7 +3,7 @@ Copyright (c) 2023, 2024 Alexei Sibidanov. The original version of this file was copied from the CORE-MATH -project (file src/binary32/cbrt/cbrtf.c, revision bc385c2). +project (file src/binary32/cbrt/cbrtf.c, revision f7c7408d). Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal @@ -61,8 +61,8 @@ __cbrtf (float x) e += 899; uint32_t et = e / 3, it = e % 3; uint64_t isc = escale[it].u; - isc += (int64_t) (et - 342) << 52; - isc |= (int64_t) sgn << 63; + isc += (uint64_t) (et - 342) << 52; + isc |= (uint64_t) sgn << 63; double cvt2 = asdouble (isc); static const double c[] = { diff --git a/sysdeps/ieee754/flt-32/s_cospif.c b/sysdeps/ieee754/flt-32/s_cospif.c index 1e83803..eb4a10f 100644 --- a/sysdeps/ieee754/flt-32/s_cospif.c +++ b/sysdeps/ieee754/flt-32/s_cospif.c @@ -3,7 +3,7 @@ Copyright (c) 2022-2025 Alexei Sibidanov. The original version of this file was copied from the CORE-MATH -project (src/binary32/cospi/cospif.c, revision f786e13). +project (src/binary32/cospi/cospif.c, revision bbfabd99). Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal @@ -58,10 +58,10 @@ __cospif (float x) { if (__glibc_unlikely (p > 63)) return 1.0f; - int32_t iq = m << (p - 32); + int32_t iq = (uint32_t)m << (p - 32); return S[(iq + 32) & 127]; } - int32_t k = m << p; + int32_t k = (uint32_t)m << p; if (__glibc_unlikely (k == 0)) { int32_t iq = m >> (32 - p); diff --git a/sysdeps/ieee754/flt-32/s_erfcf.c b/sysdeps/ieee754/flt-32/s_erfcf.c index 3dae2a0..955f129 100644 --- a/sysdeps/ieee754/flt-32/s_erfcf.c +++ b/sysdeps/ieee754/flt-32/s_erfcf.c @@ -3,7 +3,7 @@ Copyright (c) 2023, 2024 Alexei Sibidanov. This file is part of the CORE-MATH project -project (file src/binary32/erfc/erfcf.c revision bc385c2). +project (file src/binary32/erfc/erfcf.c revision d0a2be20). Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal @@ -142,7 +142,7 @@ __erfcf (float xf) const double ln2l = 0x1.cf79abd6f5dc8p-47; uint64_t jt = asuint64 (fma (x2, iln2, -(1024 + 0x1p-8))); int64_t j = (int64_t) (jt << 12) >> 48; - double S = asdouble (((j >> 7) + (0x3ff | sgn << 11)) << 52); + double S = asdouble ((uint64_t)((j >> 7) + (0x3ff | sgn << 11)) << 52); static const double ch[] = { -0x1.ffffffffff333p-2, 0x1.5555555556a14p-3, -0x1.55556666659b4p-5, diff --git a/sysdeps/ieee754/flt-32/s_modff.c b/sysdeps/ieee754/flt-32/s_modff.c index ad2e91d..965136b 100644 --- a/sysdeps/ieee754/flt-32/s_modff.c +++ b/sysdeps/ieee754/flt-32/s_modff.c @@ -1,54 +1,69 @@ -/* s_modff.c -- float version of s_modf.c. - */ - -/* - * ==================================================== - * 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. - * ==================================================== - */ +/* Extract signed integral and fractional values. + Copyright (C) 1993-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 <libm-alias-float.h> - -static const float one = 1.0; +#include "math_config.h" +#include <math-use-builtins-trunc.h> float -__modff(float x, float *iptr) +__modff (float x, float *iptr) { - int32_t i0,j0; - uint32_t i; - GET_FLOAT_WORD(i0,x); - j0 = ((i0>>23)&0xff)-0x7f; /* exponent of x */ - if(__builtin_expect(j0<23, 1)) { /* integer part in x */ - if(j0<0) { /* |x|<1 */ - SET_FLOAT_WORD(*iptr,i0&0x80000000); /* *iptr = +-0 */ - return x; - } else { - i = (0x007fffff)>>j0; - if((i0&i)==0) { /* x is integral */ - uint32_t ix; - *iptr = x; - GET_FLOAT_WORD(ix,x); - SET_FLOAT_WORD(x,ix&0x80000000); /* return +-0 */ - return x; - } else { - SET_FLOAT_WORD(*iptr,i0&(~i)); - return x - *iptr; - } - } - } else { /* no fraction part */ - *iptr = x*one; - /* We must handle NaNs separately. */ - if (j0 == 0x80 && (i0 & 0x7fffff)) - return x*one; - SET_FLOAT_WORD(x,i0&0x80000000); /* return +-0 */ - return x; + uint32_t t = asuint (x); +#if USE_TRUNCF_BUILTIN + if (is_inf (t)) + { + *iptr = x; + return copysignf (0.0, x); + } + *iptr = truncf (x); + return copysignf (x - *iptr, x); +#else + int e = get_exponent (t); + /* No fraction part. */ + if (e < MANTISSA_WIDTH) + { + if (e < 0) + { + /* |x|<1 -> *iptr = +-0 */ + *iptr = asfloat (t & SIGN_MASK); + return x; } + + uint32_t i = MANTISSA_MASK >> e; + if ((t & i) == 0) + { + /* x in integral, return +-0 */ + *iptr = x; + return asfloat (t & SIGN_MASK); + } + + *iptr = asfloat (t & ~i); + return x - *iptr; + } + + /* Set invalid operation for sNaN. */ + *iptr = x * 1.0f; + if ((e == 0x80) && (t & MANTISSA_MASK)) + return *iptr; + return asfloat (t & SIGN_MASK); +#endif } +#ifndef __modff libm_alias_float (__modf, modf) +#endif diff --git a/sysdeps/ieee754/flt-32/s_sinpif.c b/sysdeps/ieee754/flt-32/s_sinpif.c index 99a8bbb..c0d15e7 100644 --- a/sysdeps/ieee754/flt-32/s_sinpif.c +++ b/sysdeps/ieee754/flt-32/s_sinpif.c @@ -3,7 +3,7 @@ Copyright (c) 2022-2025 Alexei Sibidanov. The original version of this file was copied from the CORE-MATH -project (src/binary32/sinpi/sinpif.c, revision f786e13). +project (src/binary32/sinpi/sinpif.c, revision bbfabd99. Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal @@ -51,7 +51,7 @@ __sinpif (float x) { if (__glibc_unlikely (s < -6)) return copysignf (0.0f, x); - int32_t iq = m << (-s - 1); + int32_t iq = (uint32_t)m << (-s - 1); iq &= 127; if (iq == 0 || iq == 64) return copysignf (0.0f, x); @@ -63,10 +63,10 @@ __sinpif (float x) return z * (0x1.921fb54442d18p+1 + z2 * (-0x1.4abbce625be53p+2)); } int32_t si = 25 - s; - if (__glibc_unlikely (si >= 0 && (m << si) == 0)) + if (__glibc_unlikely (si >= 0 && ((uint32_t)m << si) == 0)) return copysignf (0.0f, x); - int32_t k = m << (31 - s); + int32_t k = (uint32_t)m << (31 - s); double z = k, z2 = z * z; double fs = SN[0] + z2 * (SN[1] + z2 * SN[2]); double fc = CN[0] + z2 * (CN[1] + z2 * CN[2]); diff --git a/sysdeps/ieee754/flt-32/w_ilogbf-impl.h b/sysdeps/ieee754/flt-32/w_ilogbf-impl.h new file mode 100644 index 0000000..5aa8bf0 --- /dev/null +++ b/sysdeps/ieee754/flt-32/w_ilogbf-impl.h @@ -0,0 +1,38 @@ +/* Get integer exponent of a floating-point value. + Copyright (C) 1999-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/>. */ + +static inline RET_TYPE +IMPL_NAME (float x) +{ + uint32_t ux = asuint (x); + int ex = (ux & ~SIGN_MASK) >> MANTISSA_WIDTH; + if (__glibc_unlikely (ex == 0)) + { + /* Zero or subnormal. + Clear sign and exponent. */ + ux <<= 1 + EXPONENT_WIDTH; + if (ux == 0) + return RET_INVALID (RET_LOGB0); + /* subnormal */ + return (RET_TYPE)-127 - stdc_leading_zeros (ux); + } + if (__glibc_unlikely (ex == EXPONENT_MASK >> MANTISSA_WIDTH)) + /* NaN or Inf */ + return RET_INVALID (ux << (1 + EXPONENT_WIDTH) ? RET_LOGBNAN : RET_LOGMAX); + return ex - 127; +} diff --git a/sysdeps/ieee754/flt-32/w_ilogbf.c b/sysdeps/ieee754/flt-32/w_ilogbf.c new file mode 100644 index 0000000..4e2a707 --- /dev/null +++ b/sysdeps/ieee754/flt-32/w_ilogbf.c @@ -0,0 +1,53 @@ +/* Get integer exponent of a floating-point value. + Copyright (C) 1999-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 <limits.h> +#include <math.h> +#include <stdbit.h> +#include <libm-alias-float.h> +#include <math-type-macros-float.h> +#include "math_config.h" + +#ifdef DEF_AS_LLOGBF +# define DECL_NAME __llogb +# define FUNC_NAME llogb +# define RET_TYPE long int +# define RET_LOGB0 FP_LLOGB0 +# define RET_LOGBNAN FP_LLOGBNAN +# define RET_LOGMAX LONG_MAX +# define RET_INVALID __math_invalidf_li +#else +# define DECL_NAME __ilogb +# define FUNC_NAME ilogb +# define RET_TYPE int +# define RET_LOGB0 FP_ILOGB0 +# define RET_LOGBNAN FP_ILOGBNAN +# define RET_LOGMAX INT_MAX +# define RET_INVALID __math_invalidf_i +#endif +#define __IMPL_NAME(x,y) x ## _ ## y +#define _IMPL_NAME(x,y) __IMPL_NAME(x,y) +#define IMPL_NAME _IMPL_NAME(FUNC_NAME, impl) +#include <w_ilogbf-impl.h> + +RET_TYPE +M_DECL_FUNC (DECL_NAME) (float x) +{ + return IMPL_NAME (x); +} +libm_alias_float (DECL_NAME, FUNC_NAME); diff --git a/sysdeps/ieee754/flt-32/w_llogbf.c b/sysdeps/ieee754/flt-32/w_llogbf.c new file mode 100644 index 0000000..8676434 --- /dev/null +++ b/sysdeps/ieee754/flt-32/w_llogbf.c @@ -0,0 +1,2 @@ +#define DEF_AS_LLOGBF +#include "w_ilogbf.c" diff --git a/sysdeps/ieee754/ldbl-128/Makefile b/sysdeps/ieee754/ldbl-128/Makefile index 5476a55..e666bdc 100644 --- a/sysdeps/ieee754/ldbl-128/Makefile +++ b/sysdeps/ieee754/ldbl-128/Makefile @@ -83,7 +83,7 @@ CFLAGS-w_j1l.c += -fno-builtin-j1f64x -fno-builtin-j1f128 CFLAGS-w_jnl.c += -fno-builtin-jnf64x -fno-builtin-jnf128 CFLAGS-s_ldexpl.c += -fno-builtin-ldexpf64x -fno-builtin-ldexpf128 CFLAGS-w_lgammal.c += -fno-builtin-lgammaf64x -fno-builtin-lgammaf128 -CFLAGS-w_lgammal_r.c += -fno-builtin-lgammaf64x_r +CFLAGS-w_lgammal_r.c += -fno-builtin-lgammaf64x_r -fno-builtin-lgammaf128_r CFLAGS-w_llogbl.c += -fno-builtin-llogbf64x -fno-builtin-llogbf128 CFLAGS-s_llrintl.c += -fno-builtin-llrintf64x -fno-builtin-llrintf128 CFLAGS-s_llroundl.c += -fno-builtin-llroundf64x -fno-builtin-llroundf128 diff --git a/sysdeps/ieee754/ldbl-128ibm-compat/Versions b/sysdeps/ieee754/ldbl-128ibm-compat/Versions index cd39b6a..ae4bd5b 100644 --- a/sysdeps/ieee754/ldbl-128ibm-compat/Versions +++ b/sysdeps/ieee754/ldbl-128ibm-compat/Versions @@ -154,8 +154,10 @@ libm { __tanpiieee128; } GLIBC_2.42 { + __compoundnieee128; __pownieee128; __powrieee128; + __rootnieee128; __rsqrtieee128; } } diff --git a/sysdeps/ieee754/ldbl-opt/Makefile b/sysdeps/ieee754/ldbl-opt/Makefile index beaed61..ef7da1f 100644 --- a/sysdeps/ieee754/ldbl-opt/Makefile +++ b/sysdeps/ieee754/ldbl-opt/Makefile @@ -42,6 +42,7 @@ libnldbl-calls = \ cimag \ clog \ clog10 \ + compoundn \ conj \ copysign \ cos \ @@ -180,6 +181,7 @@ libnldbl-calls = \ remainder \ remquo \ rint \ + rootn \ round \ roundeven \ rsqrt \ @@ -264,7 +266,7 @@ extra-objs += $(addsuffix .oS, $(libnldbl-routines)) CFLAGS-nldbl-acos.c = -fno-builtin-acosl CFLAGS-nldbl-acosh.c = -fno-builtin-acoshl -CFLAGS-nldbl-acospi.c = -fno-builtin-acospi +CFLAGS-nldbl-acospi.c = -fno-builtin-acospil CFLAGS-nldbl-asin.c = -fno-builtin-asinl CFLAGS-nldbl-asinh.c = -fno-builtin-asinhl CFLAGS-nldbl-asinpi.c = -fno-builtin-asinpil @@ -290,11 +292,12 @@ CFLAGS-nldbl-cexp.c = -fno-builtin-cexpl CFLAGS-nldbl-cimag.c = -fno-builtin-cimagl CFLAGS-nldbl-clog.c = -fno-builtin-clogl CFLAGS-nldbl-clog10.c = -fno-builtin-clog10l +CFLAGS-nldbl-compoundn.c = -fno-builtin-compoundnl CFLAGS-nldbl-conj.c = -fno-builtin-conjl CFLAGS-nldbl-copysign.c = -fno-builtin-copysignl CFLAGS-nldbl-cos.c = -fno-builtin-cosl CFLAGS-nldbl-cosh.c = -fno-builtin-coshl -CFLAGS-nldbl-cospi.c = -fno-builtin-cospi +CFLAGS-nldbl-cospi.c = -fno-builtin-cospil CFLAGS-nldbl-cpow.c = -fno-builtin-cpowl CFLAGS-nldbl-cproj.c = -fno-builtin-cprojl CFLAGS-nldbl-creal.c = -fno-builtin-creall @@ -382,6 +385,7 @@ CFLAGS-nldbl-powr.c = -fno-builtin-powrl CFLAGS-nldbl-remainder.c = -fno-builtin-remainderl -fno-builtin-dreml CFLAGS-nldbl-remquo.c = -fno-builtin-remquol CFLAGS-nldbl-rint.c = -fno-builtin-rintl +CFLAGS-nldbl-rootn.c = -fno-builtin-rootnl CFLAGS-nldbl-round.c = -fno-builtin-roundl CFLAGS-nldbl-roundeven.c = -fno-builtin-roundevenl CFLAGS-nldbl-rsqrt.c = -fno-builtin-rsqrtl @@ -394,11 +398,11 @@ CFLAGS-nldbl-significand.c = -fno-builtin-significandl CFLAGS-nldbl-sin.c = -fno-builtin-sinl CFLAGS-nldbl-sincos.c = -fno-builtin-sincosl CFLAGS-nldbl-sinh.c = -fno-builtin-sinhl -CFLAGS-nldbl-sinpi.c = -fno-builtin-sinpi +CFLAGS-nldbl-sinpi.c = -fno-builtin-sinpil CFLAGS-nldbl-sqrt.c = -fno-builtin-sqrtl CFLAGS-nldbl-tan.c = -fno-builtin-tanl CFLAGS-nldbl-tanh.c = -fno-builtin-tanhl -CFLAGS-nldbl-tanpi.c = -fno-builtin-tanpi +CFLAGS-nldbl-tanpi.c = -fno-builtin-tanpil CFLAGS-nldbl-tgamma.c = -fno-builtin-tgammal CFLAGS-nldbl-totalorder.c = -fno-builtin-totalorderl CFLAGS-nldbl-totalordermag.c = -fno-builtin-totalordermagl diff --git a/sysdeps/ieee754/ldbl-opt/nldbl-compoundn.c b/sysdeps/ieee754/ldbl-opt/nldbl-compoundn.c new file mode 100644 index 0000000..43da519 --- /dev/null +++ b/sysdeps/ieee754/ldbl-opt/nldbl-compoundn.c @@ -0,0 +1,8 @@ +#include "nldbl-compat.h" + +double +attribute_hidden +compoundnl (double x, long long int y) +{ + return compoundn (x, y); +} diff --git a/sysdeps/ieee754/ldbl-opt/nldbl-rootn.c b/sysdeps/ieee754/ldbl-opt/nldbl-rootn.c new file mode 100644 index 0000000..fb0d860 --- /dev/null +++ b/sysdeps/ieee754/ldbl-opt/nldbl-rootn.c @@ -0,0 +1,8 @@ +#include "nldbl-compat.h" + +double +attribute_hidden +rootnl (double x, long long int y) +{ + return rootn (x, y); +} |