diff options
Diffstat (limited to 'lib/msun')
| -rw-r--r-- | lib/msun/bsdsrc/b_tgamma.c | 2 | ||||
| -rw-r--r-- | lib/msun/ld128/e_rem_pio2l.h | 2 | ||||
| -rw-r--r-- | lib/msun/ld128/s_logl.c | 2 | ||||
| -rw-r--r-- | lib/msun/ld80/e_rem_pio2l.h | 2 | ||||
| -rw-r--r-- | lib/msun/ld80/s_logl.c | 2 | ||||
| -rw-r--r-- | lib/msun/src/e_fmod.c | 30 | ||||
| -rw-r--r-- | lib/msun/src/e_fmodf.c | 16 | ||||
| -rw-r--r-- | lib/msun/src/e_rem_pio2.c | 2 | ||||
| -rw-r--r-- | lib/msun/src/e_rem_pio2f.c | 2 | ||||
| -rw-r--r-- | lib/msun/src/e_remainder.c | 4 | ||||
| -rw-r--r-- | lib/msun/src/math_private.h | 35 | ||||
| -rw-r--r-- | lib/msun/src/s_ccosh.c | 23 | ||||
| -rw-r--r-- | lib/msun/src/s_ccoshf.c | 18 | ||||
| -rw-r--r-- | lib/msun/src/s_csinh.c | 23 | ||||
| -rw-r--r-- | lib/msun/src/s_csinhf.c | 18 | ||||
| -rw-r--r-- | lib/msun/src/s_ilogb.c | 13 | ||||
| -rw-r--r-- | lib/msun/src/s_ilogbf.c | 5 | ||||
| -rw-r--r-- | lib/msun/src/s_remquo.c | 30 | ||||
| -rw-r--r-- | lib/msun/src/s_remquof.c | 16 |
19 files changed, 148 insertions, 97 deletions
diff --git a/lib/msun/bsdsrc/b_tgamma.c b/lib/msun/bsdsrc/b_tgamma.c index cc9f8f70297e..5ed69a28156f 100644 --- a/lib/msun/bsdsrc/b_tgamma.c +++ b/lib/msun/bsdsrc/b_tgamma.c @@ -261,7 +261,7 @@ small_gam(double x) static double smaller_gam(double x) { - double d, rhi, rlo, t, xhi, xlo; + double d, t, xhi, xlo; struct Double r; if (x < x0 + left) { diff --git a/lib/msun/ld128/e_rem_pio2l.h b/lib/msun/ld128/e_rem_pio2l.h index 35ed0b865a7c..d38a027dc215 100644 --- a/lib/msun/ld128/e_rem_pio2l.h +++ b/lib/msun/ld128/e_rem_pio2l.h @@ -56,7 +56,7 @@ pio2_2t = 2.0670321098263988236496903051604844e-43L, /* 0x127044533e63a0105df5 pio2_3 = 2.0670321098263988236499468110329591e-43L, /* 0x127044533e63a0105e00000000000.0p-254 */ pio2_3t = -2.5650587247459238361625433492959285e-65L; /* -0x159c4ec64ddaeb5f78671cbfb2210.0p-327 */ -static inline __always_inline int +static __always_inline int __ieee754_rem_pio2l(long double x, long double *y) { union IEEEl2bits u,u1; diff --git a/lib/msun/ld128/s_logl.c b/lib/msun/ld128/s_logl.c index 8961dd0c96db..734e0a8c2d4d 100644 --- a/lib/msun/ld128/s_logl.c +++ b/lib/msun/ld128/s_logl.c @@ -445,7 +445,7 @@ struct ld { #endif #ifdef STRUCT_RETURN -static inline __always_inline void +static __always_inline void k_logl(long double x, struct ld *rp) #else long double diff --git a/lib/msun/ld80/e_rem_pio2l.h b/lib/msun/ld80/e_rem_pio2l.h index d9e4d3a0941f..ee35b28c83fd 100644 --- a/lib/msun/ld80/e_rem_pio2l.h +++ b/lib/msun/ld80/e_rem_pio2l.h @@ -68,7 +68,7 @@ pio2_2t = 6.36831716351095013979e-25L, /* 0xc51701b839a25205.0p-144 */ pio2_3t = -2.75299651904407171810e-37L; /* -0xbb5bf6c7ddd660ce.0p-185 */ #endif -static inline __always_inline int +static __always_inline int __ieee754_rem_pio2l(long double x, long double *y) { union IEEEl2bits u,u1; diff --git a/lib/msun/ld80/s_logl.c b/lib/msun/ld80/s_logl.c index 459374d7d164..2305a906bd4e 100644 --- a/lib/msun/ld80/s_logl.c +++ b/lib/msun/ld80/s_logl.c @@ -445,7 +445,7 @@ struct ld { #endif #ifdef STRUCT_RETURN -static inline __always_inline void +static __always_inline void k_logl(long double x, struct ld *rp) #else long double diff --git a/lib/msun/src/e_fmod.c b/lib/msun/src/e_fmod.c index fdfb56c2a475..973ab2ef30c4 100644 --- a/lib/msun/src/e_fmod.c +++ b/lib/msun/src/e_fmod.c @@ -28,14 +28,14 @@ static const double one = 1.0, Zero[] = {0.0, -0.0,}; double fmod(double x, double y) { - int32_t n,hx,hy,hz,ix,iy,sx,i; - u_int32_t lx,ly,lz; + int32_t hx, hy, hz, ix, iy, n, sx; + u_int32_t lx, ly, lz; EXTRACT_WORDS(hx,lx,x); EXTRACT_WORDS(hy,ly,y); sx = hx&0x80000000; /* sign of x */ - hx ^=sx; /* |x| */ - hy &= 0x7fffffff; /* |y| */ + hx ^= sx; /* |x| */ + hy &= 0x7fffffff; /* |y| */ /* purge off exception values */ if((hy|ly)==0||(hx>=0x7ff00000)|| /* y=0,or x not finite */ @@ -48,22 +48,16 @@ fmod(double x, double y) } /* determine ix = ilogb(x) */ - if(hx<0x00100000) { /* subnormal x */ - if(hx==0) { - for (ix = -1043, i=lx; i>0; i<<=1) ix -=1; - } else { - for (ix = -1022,i=(hx<<11); i>0; i<<=1) ix -=1; - } - } else ix = (hx>>20)-1023; + if(hx<0x00100000) + ix = subnormal_ilogb(hx, lx); + else + ix = (hx>>20)-1023; /* determine iy = ilogb(y) */ - if(hy<0x00100000) { /* subnormal y */ - if(hy==0) { - for (iy = -1043, i=ly; i>0; i<<=1) iy -=1; - } else { - for (iy = -1022,i=(hy<<11); i>0; i<<=1) iy -=1; - } - } else iy = (hy>>20)-1023; + if(hy<0x00100000) + iy = subnormal_ilogb(hy, ly); + else + iy = (hy>>20)-1023; /* set up {hx,lx}, {hy,ly} and align y to x */ if(ix >= -1022) diff --git a/lib/msun/src/e_fmodf.c b/lib/msun/src/e_fmodf.c index 0e6633fbe739..fc6fee879f52 100644 --- a/lib/msun/src/e_fmodf.c +++ b/lib/msun/src/e_fmodf.c @@ -28,7 +28,7 @@ static const float one = 1.0, Zero[] = {0.0, -0.0,}; float fmodf(float x, float y) { - int32_t n,hx,hy,hz,ix,iy,sx,i; + int32_t hx, hy, hz, ix, iy, n, sx; GET_FLOAT_WORD(hx,x); GET_FLOAT_WORD(hy,y); @@ -45,14 +45,16 @@ fmodf(float x, float y) return Zero[(u_int32_t)sx>>31]; /* |x|=|y| return x*0*/ /* determine ix = ilogb(x) */ - if(hx<0x00800000) { /* subnormal x */ - for (ix = -126,i=(hx<<8); i>0; i<<=1) ix -=1; - } else ix = (hx>>23)-127; + if(hx<0x00800000) + ix = subnormal_ilogbf(hx); + else + ix = (hx>>23)-127; /* determine iy = ilogb(y) */ - if(hy<0x00800000) { /* subnormal y */ - for (iy = -126,i=(hy<<8); i>=0; i<<=1) iy -=1; - } else iy = (hy>>23)-127; + if(hy<0x00800000) + iy = subnormal_ilogbf(hy); + else + iy = (hy>>23)-127; /* set up {hx,lx}, {hy,ly} and align y to x */ if(ix >= -126) diff --git a/lib/msun/src/e_rem_pio2.c b/lib/msun/src/e_rem_pio2.c index ef4107af94cb..d31cce4d6912 100644 --- a/lib/msun/src/e_rem_pio2.c +++ b/lib/msun/src/e_rem_pio2.c @@ -47,7 +47,7 @@ pio2_3 = 2.02226624871116645580e-21, /* 0x3BA3198A, 0x2E000000 */ pio2_3t = 8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */ #ifdef INLINE_REM_PIO2 -static __inline __always_inline +static __always_inline #endif int __ieee754_rem_pio2(double x, double *y) diff --git a/lib/msun/src/e_rem_pio2f.c b/lib/msun/src/e_rem_pio2f.c index 26f6bc85791b..1740b4f49a58 100644 --- a/lib/msun/src/e_rem_pio2f.c +++ b/lib/msun/src/e_rem_pio2f.c @@ -39,7 +39,7 @@ pio2_1 = 1.57079631090164184570e+00, /* 0x3FF921FB, 0x50000000 */ pio2_1t = 1.58932547735281966916e-08; /* 0x3E5110b4, 0x611A6263 */ #ifdef INLINE_REM_PIO2F -static __inline __always_inline +static __always_inline #endif int __ieee754_rem_pio2f(float x, double *y) diff --git a/lib/msun/src/e_remainder.c b/lib/msun/src/e_remainder.c index bd1ce8950619..5b71088911a2 100644 --- a/lib/msun/src/e_remainder.c +++ b/lib/msun/src/e_remainder.c @@ -66,8 +66,8 @@ remainder(double x, double p) if(x>=p_half) x -= p; } } - GET_HIGH_WORD(hx,x); - if ((hx&0x7fffffff)==0) hx = 0; + EXTRACT_WORDS(hx, lx, x); + if (((hx&0x7fffffff)|lx) == 0) hx = 0; SET_HIGH_WORD(x,hx^sx); return x; } diff --git a/lib/msun/src/math_private.h b/lib/msun/src/math_private.h index b71a8d03c0a1..6513c2ce14fc 100644 --- a/lib/msun/src/math_private.h +++ b/lib/msun/src/math_private.h @@ -740,6 +740,41 @@ irintl(long double x) (ar) = (x) - (ai); \ } while (0) +/* + * For a subnormal double entity split into high and low parts, compute ilogb. + */ +static inline int32_t +subnormal_ilogb(int32_t hi, int32_t lo) +{ + int32_t j; + uint32_t i; + + j = -1022; + if (hi == 0) { + j -= 21; + i = (uint32_t)lo; + } else + i = (uint32_t)hi << 11; + + for (; i < 0x7fffffff; i <<= 1) j -= 1; + + return (j); +} + +/* + * For a subnormal float entity represented as an int32_t, compute ilogb. + */ +static inline int32_t +subnormal_ilogbf(int32_t hx) +{ + int32_t j; + uint32_t i; + i = (uint32_t) hx << 8; + for (j = -126; i < 0x7fffffff; i <<= 1) j -=1; + + return (j); +} + #ifdef DEBUG #if defined(__amd64__) || defined(__i386__) #define breakpoint() asm("int $3") diff --git a/lib/msun/src/s_ccosh.c b/lib/msun/src/s_ccosh.c index 95ed3a32ddd7..145d9673f392 100644 --- a/lib/msun/src/s_ccosh.c +++ b/lib/msun/src/s_ccosh.c @@ -1,7 +1,7 @@ /*- * SPDX-License-Identifier: BSD-2-Clause * - * Copyright (c) 2005 Bruce D. Evans and Steven G. Kargl + * Copyright (c) 2005-2025 Bruce D. Evans and Steven G. Kargl * All rights reserved. * * Redistribution and use in source and binary forms, with or without @@ -49,7 +49,7 @@ static const double huge = 0x1p1023; double complex ccosh(double complex z) { - double x, y, h; + double c, h, s, x, y; int32_t hx, hy, ix, iy, lx, ly; x = creal(z); @@ -65,14 +65,16 @@ ccosh(double complex z) if (ix < 0x7ff00000 && iy < 0x7ff00000) { if ((iy | ly) == 0) return (CMPLX(cosh(x), x * y)); + + sincos(y, &s, &c); if (ix < 0x40360000) /* |x| < 22: normal case */ - return (CMPLX(cosh(x) * cos(y), sinh(x) * sin(y))); + return (CMPLX(cosh(x) * c, sinh(x) * s)); /* |x| >= 22, so cosh(x) ~= exp(|x|) */ if (ix < 0x40862e42) { /* x < 710: exp(|x|) won't overflow */ - h = exp(fabs(x)) * 0.5; - return (CMPLX(h * cos(y), copysign(h, x) * sin(y))); + h = exp(fabs(x)) / 2; + return (CMPLX(h * c, copysign(h, x) * s)); } else if (ix < 0x4096bbaa) { /* x < 1455: scale to avoid overflow */ z = __ldexp_cexp(CMPLX(fabs(x), y), -1); @@ -80,7 +82,7 @@ ccosh(double complex z) } else { /* x >= 1455: the result always overflows */ h = huge * x; - return (CMPLX(h * h * cos(y), h * sin(y))); + return (CMPLX(h * h * c, h * s)); } } @@ -130,7 +132,9 @@ ccosh(double complex z) if (ix == 0x7ff00000 && lx == 0) { if (iy >= 0x7ff00000) return (CMPLX(INFINITY, x * (y - y))); - return (CMPLX(INFINITY * cos(y), x * sin(y))); + + sincos(y, &s, &c); + return (CMPLX(INFINITY * c, x * s)); } /* @@ -155,3 +159,8 @@ ccos(double complex z) /* ccos(z) = ccosh(I * z) */ return (ccosh(CMPLX(-cimag(z), creal(z)))); } + +#if (LDBL_MANT_DIG == 53) +__weak_reference(ccosh, ccoshl); +__weak_reference(ccos, ccosl); +#endif diff --git a/lib/msun/src/s_ccoshf.c b/lib/msun/src/s_ccoshf.c index ba97a390c832..2851157ead42 100644 --- a/lib/msun/src/s_ccoshf.c +++ b/lib/msun/src/s_ccoshf.c @@ -1,7 +1,7 @@ /*- * SPDX-License-Identifier: BSD-2-Clause * - * Copyright (c) 2005 Bruce D. Evans and Steven G. Kargl + * Copyright (c) 2005-2025 Bruce D. Evans and Steven G. Kargl * All rights reserved. * * Redistribution and use in source and binary forms, with or without @@ -41,7 +41,7 @@ static const float huge = 0x1p127; float complex ccoshf(float complex z) { - float x, y, h; + float c, h, s, x, y; int32_t hx, hy, ix, iy; x = crealf(z); @@ -56,14 +56,16 @@ ccoshf(float complex z) if (ix < 0x7f800000 && iy < 0x7f800000) { if (iy == 0) return (CMPLXF(coshf(x), x * y)); + + sincosf(y, &s, &c); if (ix < 0x41100000) /* |x| < 9: normal case */ - return (CMPLXF(coshf(x) * cosf(y), sinhf(x) * sinf(y))); + return (CMPLXF(coshf(x) * c, sinhf(x) * s)); /* |x| >= 9, so cosh(x) ~= exp(|x|) */ if (ix < 0x42b17218) { /* x < 88.7: expf(|x|) won't overflow */ - h = expf(fabsf(x)) * 0.5F; - return (CMPLXF(h * cosf(y), copysignf(h, x) * sinf(y))); + h = expf(fabsf(x)) / 2; + return (CMPLXF(h * c, copysignf(h, x) * s)); } else if (ix < 0x4340b1e7) { /* x < 192.7: scale to avoid overflow */ z = __ldexp_cexpf(CMPLXF(fabsf(x), y), -1); @@ -71,7 +73,7 @@ ccoshf(float complex z) } else { /* x >= 192.7: the result always overflows */ h = huge * x; - return (CMPLXF(h * h * cosf(y), h * sinf(y))); + return (CMPLXF(h * h * c, h * s)); } } @@ -87,7 +89,9 @@ ccoshf(float complex z) if (ix == 0x7f800000) { if (iy >= 0x7f800000) return (CMPLXF(INFINITY, x * (y - y))); - return (CMPLXF(INFINITY * cosf(y), x * sinf(y))); + + sincosf(y, &s, &c); + return (CMPLXF(INFINITY * c, x * s)); } return (CMPLXF(((long double)x * x) * (y - y), diff --git a/lib/msun/src/s_csinh.c b/lib/msun/src/s_csinh.c index 1bd78b1e49bf..8ea97970af7f 100644 --- a/lib/msun/src/s_csinh.c +++ b/lib/msun/src/s_csinh.c @@ -1,7 +1,7 @@ /*- * SPDX-License-Identifier: BSD-2-Clause * - * Copyright (c) 2005 Bruce D. Evans and Steven G. Kargl + * Copyright (c) 2005-2025 Bruce D. Evans and Steven G. Kargl * All rights reserved. * * Redistribution and use in source and binary forms, with or without @@ -49,7 +49,7 @@ static const double huge = 0x1p1023; double complex csinh(double complex z) { - double x, y, h; + double c, h, s, x, y; int32_t hx, hy, ix, iy, lx, ly; x = creal(z); @@ -65,14 +65,16 @@ csinh(double complex z) if (ix < 0x7ff00000 && iy < 0x7ff00000) { if ((iy | ly) == 0) return (CMPLX(sinh(x), y)); + + sincos(y, &s, &c); if (ix < 0x40360000) /* |x| < 22: normal case */ - return (CMPLX(sinh(x) * cos(y), cosh(x) * sin(y))); + return (CMPLX(sinh(x) * c, cosh(x) * s)); /* |x| >= 22, so cosh(x) ~= exp(|x|) */ if (ix < 0x40862e42) { /* x < 710: exp(|x|) won't overflow */ - h = exp(fabs(x)) * 0.5; - return (CMPLX(copysign(h, x) * cos(y), h * sin(y))); + h = exp(fabs(x)) / 2; + return (CMPLX(copysign(h, x) * c, h * s)); } else if (ix < 0x4096bbaa) { /* x < 1455: scale to avoid overflow */ z = __ldexp_cexp(CMPLX(fabs(x), y), -1); @@ -80,7 +82,7 @@ csinh(double complex z) } else { /* x >= 1455: the result always overflows */ h = huge * x; - return (CMPLX(h * cos(y), h * h * sin(y))); + return (CMPLX(h * c, h * h * s)); } } @@ -129,7 +131,9 @@ csinh(double complex z) if (ix == 0x7ff00000 && lx == 0) { if (iy >= 0x7ff00000) return (CMPLX(x, y - y)); - return (CMPLX(x * cos(y), INFINITY * sin(y))); + + sincos(y, &s, &c); + return (CMPLX(x * c, INFINITY * s)); } /* @@ -155,3 +159,8 @@ csin(double complex z) z = csinh(CMPLX(cimag(z), creal(z))); return (CMPLX(cimag(z), creal(z))); } + +#if (LDBL_MANT_DIG == 53) +__weak_reference(csinh, csinhl); +__weak_reference(csin, csinl); +#endif diff --git a/lib/msun/src/s_csinhf.c b/lib/msun/src/s_csinhf.c index b1f333955e53..530ac63b6a27 100644 --- a/lib/msun/src/s_csinhf.c +++ b/lib/msun/src/s_csinhf.c @@ -1,7 +1,7 @@ /*- * SPDX-License-Identifier: BSD-2-Clause * - * Copyright (c) 2005 Bruce D. Evans and Steven G. Kargl + * Copyright (c) 2005-2025 Bruce D. Evans and Steven G. Kargl * All rights reserved. * * Redistribution and use in source and binary forms, with or without @@ -41,7 +41,7 @@ static const float huge = 0x1p127; float complex csinhf(float complex z) { - float x, y, h; + float c, h, s, x, y; int32_t hx, hy, ix, iy; x = crealf(z); @@ -56,14 +56,16 @@ csinhf(float complex z) if (ix < 0x7f800000 && iy < 0x7f800000) { if (iy == 0) return (CMPLXF(sinhf(x), y)); + + sincosf(y, &s, &c); if (ix < 0x41100000) /* |x| < 9: normal case */ - return (CMPLXF(sinhf(x) * cosf(y), coshf(x) * sinf(y))); + return (CMPLXF(sinhf(x) * c, coshf(x) * s)); /* |x| >= 9, so cosh(x) ~= exp(|x|) */ if (ix < 0x42b17218) { /* x < 88.7: expf(|x|) won't overflow */ - h = expf(fabsf(x)) * 0.5F; - return (CMPLXF(copysignf(h, x) * cosf(y), h * sinf(y))); + h = expf(fabsf(x)) / 2; + return (CMPLXF(copysignf(h, x) * c, h * s)); } else if (ix < 0x4340b1e7) { /* x < 192.7: scale to avoid overflow */ z = __ldexp_cexpf(CMPLXF(fabsf(x), y), -1); @@ -71,7 +73,7 @@ csinhf(float complex z) } else { /* x >= 192.7: the result always overflows */ h = huge * x; - return (CMPLXF(h * cosf(y), h * h * sinf(y))); + return (CMPLXF(h * c, h * h * s)); } } @@ -87,7 +89,9 @@ csinhf(float complex z) if (ix == 0x7f800000) { if (iy >= 0x7f800000) return (CMPLXF(x, y - y)); - return (CMPLXF(x * cosf(y), INFINITY * sinf(y))); + + sincosf(y, &s, &c); + return (CMPLXF(x * c, INFINITY * s)); } return (CMPLXF(((long double)x + x) * (y - y), diff --git a/lib/msun/src/s_ilogb.c b/lib/msun/src/s_ilogb.c index 0b076edbd9b5..a958182a2a35 100644 --- a/lib/msun/src/s_ilogb.c +++ b/lib/msun/src/s_ilogb.c @@ -23,21 +23,18 @@ #include "math.h" #include "math_private.h" - int ilogb(double x) +int +ilogb(double x) { - int32_t hx,lx,ix; + int32_t hx, ix, lx; EXTRACT_WORDS(hx,lx,x); hx &= 0x7fffffff; if(hx<0x00100000) { if((hx|lx)==0) return 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; - } + else + ix = subnormal_ilogb(hx, lx); return ix; } else if (hx<0x7ff00000) return (hx>>20)-1023; diff --git a/lib/msun/src/s_ilogbf.c b/lib/msun/src/s_ilogbf.c index ff3df1fc5b90..75508ab77120 100644 --- a/lib/msun/src/s_ilogbf.c +++ b/lib/msun/src/s_ilogbf.c @@ -19,7 +19,8 @@ #include "math.h" #include "math_private.h" - int ilogbf(float x) +int +ilogbf(float x) { int32_t hx,ix; @@ -29,7 +30,7 @@ if(hx==0) return FP_ILOGB0; else /* subnormal x */ - for (ix = -126,hx<<=8; hx>0; hx<<=1) ix -=1; + ix = subnormal_ilogbf(hx); return ix; } else if (hx<0x7f800000) return (hx>>23)-127; diff --git a/lib/msun/src/s_remquo.c b/lib/msun/src/s_remquo.c index e3aac25230e0..fba788ed5104 100644 --- a/lib/msun/src/s_remquo.c +++ b/lib/msun/src/s_remquo.c @@ -5,7 +5,7 @@ * * Developed at SunSoft, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice + * software is freely granted, provided that this notice * is preserved. * ==================================================== */ @@ -29,7 +29,7 @@ static const double Zero[] = {0.0, -0.0,}; double remquo(double x, double y, int *quo) { - int32_t n,hx,hy,hz,ix,iy,sx,i; + int32_t hx,hy,hz,ix,iy,n,sx; u_int32_t lx,ly,lz,q,sxy; EXTRACT_WORDS(hx,lx,x); @@ -55,25 +55,19 @@ remquo(double x, double y, int *quo) } /* determine ix = ilogb(x) */ - if(hx<0x00100000) { /* subnormal x */ - if(hx==0) { - for (ix = -1043, i=lx; i>0; i<<=1) ix -=1; - } else { - for (ix = -1022,i=(hx<<11); i>0; i<<=1) ix -=1; - } - } else ix = (hx>>20)-1023; + if(hx<0x00100000) + ix = subnormal_ilogb(hx, lx); + else + ix = (hx>>20)-1023; /* determine iy = ilogb(y) */ - if(hy<0x00100000) { /* subnormal y */ - if(hy==0) { - for (iy = -1043, i=ly; i>0; i<<=1) iy -=1; - } else { - for (iy = -1022,i=(hy<<11); i>0; i<<=1) iy -=1; - } - } else iy = (hy>>20)-1023; + if(hy<0x00100000) + iy = subnormal_ilogb(hy, ly); + else + iy = (hy>>20)-1023; /* set up {hx,lx}, {hy,ly} and align y to x */ - if(ix >= -1022) + if(ix >= -1022) hx = 0x00100000|(0x000fffff&hx); else { /* subnormal x, shift x to normal */ n = -1022-ix; @@ -85,7 +79,7 @@ remquo(double x, double y, int *quo) lx = 0; } } - if(iy >= -1022) + if(iy >= -1022) hy = 0x00100000|(0x000fffff&hy); else { /* subnormal y, shift y to normal */ n = -1022-iy; diff --git a/lib/msun/src/s_remquof.c b/lib/msun/src/s_remquof.c index c42bd8c4320d..39b87a91845f 100644 --- a/lib/msun/src/s_remquof.c +++ b/lib/msun/src/s_remquof.c @@ -27,7 +27,7 @@ static const float Zero[] = {0.0, -0.0,}; float remquof(float x, float y, int *quo) { - int32_t n,hx,hy,hz,ix,iy,sx,i; + int32_t hx, hy, hz, ix, iy, n, sx; u_int32_t q,sxy; GET_FLOAT_WORD(hx,x); @@ -49,14 +49,16 @@ remquof(float x, float y, int *quo) } /* determine ix = ilogb(x) */ - if(hx<0x00800000) { /* subnormal x */ - for (ix = -126,i=(hx<<8); i>0; i<<=1) ix -=1; - } else ix = (hx>>23)-127; + if(hx<0x00800000) + ix = subnormal_ilogbf(hx); + else + ix = (hx>>23)-127; /* determine iy = ilogb(y) */ - if(hy<0x00800000) { /* subnormal y */ - for (iy = -126,i=(hy<<8); i>0; i<<=1) iy -=1; - } else iy = (hy>>23)-127; + if(hy<0x00800000) + iy = subnormal_ilogbf(hy); + else + iy = (hy>>23)-127; /* set up {hx,lx}, {hy,ly} and align y to x */ if(ix >= -126) |
