Index: lib/libc/gen/modf.c =================================================================== --- lib/libc/gen/modf.c +++ lib/libc/gen/modf.c @@ -1,138 +1,2 @@ -/* @(#)s_modf.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. - * ==================================================== - */ - -#include -__FBSDID("$FreeBSD$"); - -/* - * modf(double x, double *iptr) - * return fraction part of x, and return x's integral part in *iptr. - * Method: - * Bit twiddling. - * - * Exception: - * No exception. - */ - -#include -#include -#include - -/* Bit fiddling routines copied from msun/src/math_private.h,v 1.15 */ - -#if BYTE_ORDER == BIG_ENDIAN - -typedef union -{ - double value; - struct - { - u_int32_t msw; - u_int32_t lsw; - } parts; -} ieee_double_shape_type; - -#endif - -#if BYTE_ORDER == LITTLE_ENDIAN - -typedef union -{ - double value; - struct - { - u_int32_t lsw; - u_int32_t msw; - } parts; -} ieee_double_shape_type; - -#endif - -/* Get two 32 bit ints from a double. */ - -#define EXTRACT_WORDS(ix0,ix1,d) \ -do { \ - ieee_double_shape_type ew_u; \ - ew_u.value = (d); \ - (ix0) = ew_u.parts.msw; \ - (ix1) = ew_u.parts.lsw; \ -} while (0) - -/* Get the more significant 32 bit int from a double. */ - -#define GET_HIGH_WORD(i,d) \ -do { \ - ieee_double_shape_type gh_u; \ - gh_u.value = (d); \ - (i) = gh_u.parts.msw; \ -} while (0) - -/* Set a double from two 32 bit ints. */ - -#define INSERT_WORDS(d,ix0,ix1) \ -do { \ - ieee_double_shape_type iw_u; \ - iw_u.parts.msw = (ix0); \ - iw_u.parts.lsw = (ix1); \ - (d) = iw_u.value; \ -} while (0) - -static const double one = 1.0; - -double -modf(double x, double *iptr) -{ - int32_t i0,i1,j0; - u_int32_t i; - EXTRACT_WORDS(i0,i1,x); - j0 = ((i0>>20)&0x7ff)-0x3ff; /* exponent of x */ - if(j0<20) { /* integer part in high x */ - if(j0<0) { /* |x|<1 */ - INSERT_WORDS(*iptr,i0&0x80000000,0); /* *iptr = +-0 */ - return x; - } else { - i = (0x000fffff)>>j0; - if(((i0&i)|i1)==0) { /* x is integral */ - u_int32_t high; - *iptr = x; - GET_HIGH_WORD(high,x); - INSERT_WORDS(x,high&0x80000000,0); /* return +-0 */ - return x; - } else { - INSERT_WORDS(*iptr,i0&(~i),0); - return x - *iptr; - } - } - } else if (j0>51) { /* no fraction part */ - u_int32_t high; - if (j0 == 0x400) { /* inf/NaN */ - *iptr = x; - return 0.0 / x; - } - *iptr = x*one; - GET_HIGH_WORD(high,x); - INSERT_WORDS(x,high&0x80000000,0); /* return +-0 */ - return x; - } else { /* fraction part in low x */ - i = ((u_int32_t)(0xffffffff))>>(j0-20); - if((i1&i)==0) { /* x is integral */ - u_int32_t high; - *iptr = x; - GET_HIGH_WORD(high,x); - INSERT_WORDS(x,high&0x80000000,0); /* return +-0 */ - return x; - } else { - INSERT_WORDS(*iptr,i0,i1&(~i)); - return x - *iptr; - } - } -} +/* Libc also provides modf() for backwards-compatibility. */ +#include "../../msun/src/s_modf.c" Index: lib/msun/src/e_fmod.c =================================================================== --- lib/msun/src/e_fmod.c +++ lib/msun/src/e_fmod.c @@ -1,136 +1,71 @@ - -/* @(#)e_fmod.c 1.3 95/01/18 */ -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunSoft, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - -#include -__FBSDID("$FreeBSD$"); - -/* - * __ieee754_fmod(x,y) - * Return x mod y in exact arithmetic - * Method: shift and subtract - */ - #include +#include +#include -#include "math.h" -#include "math_private.h" - -static const double one = 1.0, Zero[] = {0.0, -0.0,}; - -double -__ieee754_fmod(double x, double y) +double fmod(double x, double y) { - int32_t n,hx,hy,hz,ix,iy,sx,i; - u_int32_t lx,ly,lz; + union {double f; uint64_t i;} ux = {x}, uy = {y}; + int ex = ux.i>>52 & 0x7ff; + int ey = uy.i>>52 & 0x7ff; + int sx = ux.i>>63; + uint64_t i; - EXTRACT_WORDS(hx,lx,x); - EXTRACT_WORDS(hy,ly,y); - sx = hx&0x80000000; /* sign of x */ - hx ^=sx; /* |x| */ - hy &= 0x7fffffff; /* |y| */ + /* in the followings uxi should be ux.i, but then gcc wrongly adds */ + /* float load/store to inner loops ruining performance and code size */ + uint64_t uxi = ux.i; - /* purge off exception values */ - if((hy|ly)==0||(hx>=0x7ff00000)|| /* y=0,or x not finite */ - ((hy|((ly|-ly)>>31))>0x7ff00000)) /* or y is NaN */ - return nan_mix_op(x, y, *)/nan_mix_op(x, y, *); - if(hx<=hy) { - if((hx>31]; /* |x|=|y| return x*0*/ + if (uy.i<<1 == 0 || isnan(y) || ex == 0x7ff) + return (x*y)/(x*y); + if (uxi<<1 <= uy.i<<1) { + if (uxi<<1 == uy.i<<1) + return 0*x; + return x; } - /* 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; - - /* 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; - - /* set up {hx,lx}, {hy,ly} and align y to x */ - if(ix >= -1022) - hx = 0x00100000|(0x000fffff&hx); - else { /* subnormal x, shift x to normal */ - n = -1022-ix; - if(n<=31) { - hx = (hx<>(32-n)); - lx <<= n; - } else { - hx = lx<<(n-32); - lx = 0; - } + /* normalize x and y */ + if (!ex) { + for (i = uxi<<12; i>>63 == 0; ex--, i <<= 1); + uxi <<= -ex + 1; + } else { + uxi &= -1ULL >> 12; + uxi |= 1ULL << 52; } - if(iy >= -1022) - hy = 0x00100000|(0x000fffff&hy); - else { /* subnormal y, shift y to normal */ - n = -1022-iy; - if(n<=31) { - hy = (hy<>(32-n)); - ly <<= n; - } else { - hy = ly<<(n-32); - ly = 0; - } + if (!ey) { + for (i = uy.i<<12; i>>63 == 0; ey--, i <<= 1); + uy.i <<= -ey + 1; + } else { + uy.i &= -1ULL >> 12; + uy.i |= 1ULL << 52; } - /* fix point fmod */ - n = ix - iy; - while(n--) { - hz=hx-hy;lz=lx-ly; if(lx>31); lx = lx+lx;} - else { - if((hz|lz)==0) /* return sign(x)*0 */ - return Zero[(u_int32_t)sx>>31]; - hx = hz+hz+(lz>>31); lx = lz+lz; - } + /* x mod y */ + for (; ex > ey; ex--) { + i = uxi - uy.i; + if (i >> 63 == 0) { + if (i == 0) + return 0*x; + uxi = i; + } + uxi <<= 1; } - hz=hx-hy;lz=lx-ly; if(lx=0) {hx=hz;lx=lz;} - - /* convert back to floating value and restore the sign */ - if((hx|lx)==0) /* return sign(x)*0 */ - return Zero[(u_int32_t)sx>>31]; - while(hx<0x00100000) { /* normalize x */ - hx = hx+hx+(lx>>31); lx = lx+lx; - iy -= 1; + i = uxi - uy.i; + if (i >> 63 == 0) { + if (i == 0) + return 0*x; + uxi = i; } - if(iy>= -1022) { /* normalize output */ - hx = ((hx-0x00100000)|((iy+1023)<<20)); - INSERT_WORDS(x,hx|sx,lx); - } else { /* subnormal output */ - n = -1022 - iy; - if(n<=20) { - lx = (lx>>n)|((u_int32_t)hx<<(32-n)); - hx >>= n; - } else if (n<=31) { - lx = (hx<<(32-n))|(lx>>n); hx = sx; - } else { - lx = hx>>(n-32); hx = sx; - } - INSERT_WORDS(x,hx|sx,lx); - x *= one; /* create necessary signal */ + for (; uxi>>52 == 0; uxi <<= 1, ex--); + + /* scale result */ + if (ex > 0) { + uxi -= 1ULL << 52; + uxi |= (uint64_t)ex << 52; + } else { + uxi >>= -ex + 1; } - return x; /* exact output */ + uxi |= (uint64_t)sx << 63; + ux.i = uxi; + return ux.f; } #if (LDBL_MANT_DIG == 53) Index: lib/msun/src/e_fmodf.c =================================================================== --- lib/msun/src/e_fmodf.c +++ lib/msun/src/e_fmodf.c @@ -1,104 +1,65 @@ -/* e_fmodf.c -- float version of e_fmod.c. - * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. - */ +#include +#include -/* - * ==================================================== - * 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. - * ==================================================== - */ - -#include -__FBSDID("$FreeBSD$"); - -/* - * __ieee754_fmodf(x,y) - * Return x mod y in exact arithmetic - * Method: shift and subtract - */ - -#include "math.h" -#include "math_private.h" - -static const float one = 1.0, Zero[] = {0.0, -0.0,}; - -float -__ieee754_fmodf(float x, float y) +float fmodf(float x, float y) { - int32_t n,hx,hy,hz,ix,iy,sx,i; - - GET_FLOAT_WORD(hx,x); - GET_FLOAT_WORD(hy,y); - sx = hx&0x80000000; /* sign of x */ - hx ^=sx; /* |x| */ - hy &= 0x7fffffff; /* |y| */ + union {float f; uint32_t i;} ux = {x}, uy = {y}; + int ex = ux.i>>23 & 0xff; + int ey = uy.i>>23 & 0xff; + uint32_t sx = ux.i & 0x80000000; + uint32_t i; + uint32_t uxi = ux.i; - /* purge off exception values */ - if(hy==0||(hx>=0x7f800000)|| /* y=0,or x not finite */ - (hy>0x7f800000)) /* or y is NaN */ - return nan_mix_op(x, y, *)/nan_mix_op(x, y, *); - if(hx>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; - - /* 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 (uy.i<<1 == 0 || isnan(y) || ex == 0xff) + return (x*y)/(x*y); + if (uxi<<1 <= uy.i<<1) { + if (uxi<<1 == uy.i<<1) + return 0*x; + return x; + } - /* set up {hx,lx}, {hy,ly} and align y to x */ - if(ix >= -126) - hx = 0x00800000|(0x007fffff&hx); - else { /* subnormal x, shift x to normal */ - n = -126-ix; - hx = hx<>31 == 0; ex--, i <<= 1); + uxi <<= -ex + 1; + } else { + uxi &= -1U >> 9; + uxi |= 1U << 23; } - if(iy >= -126) - hy = 0x00800000|(0x007fffff&hy); - else { /* subnormal y, shift y to normal */ - n = -126-iy; - hy = hy<>31 == 0; ey--, i <<= 1); + uy.i <<= -ey + 1; + } else { + uy.i &= -1U >> 9; + uy.i |= 1U << 23; } - /* fix point fmod */ - n = ix - iy; - while(n--) { - hz=hx-hy; - if(hz<0){hx = hx+hx;} - else { - if(hz==0) /* return sign(x)*0 */ - return Zero[(u_int32_t)sx>>31]; - hx = hz+hz; - } + /* x mod y */ + for (; ex > ey; ex--) { + i = uxi - uy.i; + if (i >> 31 == 0) { + if (i == 0) + return 0*x; + uxi = i; + } + uxi <<= 1; } - hz=hx-hy; - if(hz>=0) {hx=hz;} - - /* convert back to floating value and restore the sign */ - if(hx==0) /* return sign(x)*0 */ - return Zero[(u_int32_t)sx>>31]; - while(hx<0x00800000) { /* normalize x */ - hx = hx+hx; - iy -= 1; + i = uxi - uy.i; + if (i >> 31 == 0) { + if (i == 0) + return 0*x; + uxi = i; } - if(iy>= -126) { /* normalize output */ - hx = ((hx-0x00800000)|((iy+127)<<23)); - SET_FLOAT_WORD(x,hx|sx); - } else { /* subnormal output */ - n = -126 - iy; - hx >>= n; - SET_FLOAT_WORD(x,hx|sx); - x *= one; /* create necessary signal */ + for (; uxi>>23 == 0; uxi <<= 1, ex--); + + /* scale result up */ + if (ex > 0) { + uxi -= 1U << 23; + uxi |= (uint32_t)ex << 23; + } else { + uxi >>= -ex + 1; } - return x; /* exact output */ + uxi |= sx; + ux.i = uxi; + return ux.f; } Index: lib/msun/src/e_fmodl.c =================================================================== --- lib/msun/src/e_fmodl.c +++ lib/msun/src/e_fmodl.c @@ -1,149 +1,100 @@ -/* @(#)e_fmod.c 1.3 95/01/18 */ -/*- - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunSoft, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ +#include "libm.h" -#include -__FBSDID("$FreeBSD$"); - -#include -#include - -#include "fpmath.h" -#include "math.h" -#include "math_private.h" - -#define BIAS (LDBL_MAX_EXP - 1) - -#if LDBL_MANL_SIZE > 32 -typedef uint64_t manl_t; -#else -typedef uint32_t manl_t; -#endif - -#if LDBL_MANH_SIZE > 32 -typedef uint64_t manh_t; -#else -typedef uint32_t manh_t; -#endif - -/* - * These macros add and remove an explicit integer bit in front of the - * fractional mantissa, if the architecture doesn't have such a bit by - * default already. - */ -#ifdef LDBL_IMPLICIT_NBIT -#define SET_NBIT(hx) ((hx) | (1ULL << LDBL_MANH_SIZE)) -#define HFRAC_BITS LDBL_MANH_SIZE -#else -#define SET_NBIT(hx) (hx) -#define HFRAC_BITS (LDBL_MANH_SIZE - 1) -#endif - -#define MANL_SHIFT (LDBL_MANL_SIZE - 1) - -static const long double one = 1.0, Zero[] = {0.0, -0.0,}; - -/* - * fmodl(x,y) - * Return x mod y in exact arithmetic - * Method: shift and subtract - * - * Assumptions: - * - The low part of the mantissa fits in a manl_t exactly. - * - The high part of the mantissa fits in an int64_t with enough room - * for an explicit integer bit in front of the fractional bits. - */ -long double -fmodl(long double x, long double y) +#if (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 +long double fmodl(long double x, long double y) { - union IEEEl2bits ux, uy; - int64_t hx,hz; /* We need a carry bit even if LDBL_MANH_SIZE is 32. */ - manh_t hy; - manl_t lx,ly,lz; - int ix,iy,n,sx; - - ux.e = x; - uy.e = y; - sx = ux.bits.sign; - - /* purge off exception values */ - if((uy.bits.exp|uy.bits.manh|uy.bits.manl)==0 || /* y=0 */ - (ux.bits.exp == BIAS + LDBL_MAX_EXP) || /* or x not finite */ - (uy.bits.exp == BIAS + LDBL_MAX_EXP && - ((uy.bits.manh&~LDBL_NBIT)|uy.bits.manl)!=0)) /* or y is NaN */ - return nan_mix_op(x, y, *)/nan_mix_op(x, y, *); - if(ux.bits.exp<=uy.bits.exp) { - if((ux.bits.exp>MANL_SHIFT); lx = lx+lx;} - else { - if ((hz|lz)==0) /* return sign(x)*0 */ - return Zero[sx]; - hx = hz+hz+(lz>>MANL_SHIFT); lx = lz+lz; - } + /* x mod y */ +#if LDBL_MANT_DIG == 64 + uint64_t i, mx, my; + mx = ux.i.m; + my = uy.i.m; + for (; ex > ey; ex--) { + i = mx - my; + if (mx >= my) { + if (i == 0) + return 0*x; + mx = 2*i; + } else if (2*mx < mx) { + mx = 2*mx - my; + } else { + mx = 2*mx; + } } - hz=hx-hy;lz=lx-ly; if(lx=0) {hx=hz;lx=lz;} - - /* convert back to floating value and restore the sign */ - if((hx|lx)==0) /* return sign(x)*0 */ - return Zero[sx]; - while(hx<(1ULL<>MANL_SHIFT); lx = lx+lx; - iy -= 1; + i = mx - my; + if (mx >= my) { + if (i == 0) + return 0*x; + mx = i; + } + for (; mx >> 63 == 0; mx *= 2, ex--); + ux.i.m = mx; +#elif LDBL_MANT_DIG == 113 + uint64_t hi, lo, xhi, xlo, yhi, ylo; + xhi = (ux.i2.hi & -1ULL>>16) | 1ULL<<48; + yhi = (uy.i2.hi & -1ULL>>16) | 1ULL<<48; + xlo = ux.i2.lo; + ylo = uy.i2.lo; + for (; ex > ey; ex--) { + hi = xhi - yhi; + lo = xlo - ylo; + if (xlo < ylo) + hi -= 1; + if (hi >> 63 == 0) { + if ((hi|lo) == 0) + return 0*x; + xhi = 2*hi + (lo>>63); + xlo = 2*lo; + } else { + xhi = 2*xhi + (xlo>>63); + xlo = 2*xlo; + } } - ux.bits.manh = hx; /* The mantissa is truncated here if needed. */ - ux.bits.manl = lx; - if (iy < LDBL_MIN_EXP) { - ux.bits.exp = iy + (BIAS + 512); - ux.e *= 0x1p-512; - } else { - ux.bits.exp = iy + BIAS; + hi = xhi - yhi; + lo = xlo - ylo; + if (xlo < ylo) + hi -= 1; + if (hi >> 63 == 0) { + if ((hi|lo) == 0) + return 0*x; + xhi = hi; + xlo = lo; } - x = ux.e * one; /* create necessary signal */ - return x; /* exact output */ + for (; xhi >> 48 == 0; xhi = 2*xhi + (xlo>>63), xlo = 2*xlo, ex--); + ux.i2.hi = xhi; + ux.i2.lo = xlo; +#endif + + /* scale result */ + if (ex <= 0) { + ux.i.se = (ex+120)|sx; + ux.f *= 0x1p-120f; + } else + ux.i.se = ex|sx; + return ux.f; } +#endif Index: lib/msun/src/e_remainder.c =================================================================== --- lib/msun/src/e_remainder.c +++ lib/msun/src/e_remainder.c @@ -1,79 +1,12 @@ +#include -/* @(#)e_remainder.c 1.3 95/01/18 */ -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunSoft, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - -#include -__FBSDID("$FreeBSD$"); - -/* __ieee754_remainder(x,p) - * Return : - * returns x REM p = x - [x/p]*p as if in infinite - * precise arithmetic, where [x/p] is the (infinite bit) - * integer nearest x/p (in half way case choose the even one). - * Method : - * Based on fmod() return x-[x/p]chopped*p exactlp. - */ - -#include - -#include "math.h" -#include "math_private.h" - -static const double zero = 0.0; - - -double -__ieee754_remainder(double x, double p) +double remainder(double x, double y) { - int32_t hx,hp; - u_int32_t sx,lx,lp; - double p_half; - - EXTRACT_WORDS(hx,lx,x); - EXTRACT_WORDS(hp,lp,p); - sx = hx&0x80000000; - hp &= 0x7fffffff; - hx &= 0x7fffffff; - - /* purge off exception values */ - if(((hp|lp)==0)|| /* p = 0 */ - (hx>=0x7ff00000)|| /* x not finite */ - ((hp>=0x7ff00000)&& /* p is NaN */ - (((hp-0x7ff00000)|lp)!=0))) - return nan_mix_op(x, p, *)/nan_mix_op(x, p, *); - - - if (hp<=0x7fdfffff) x = __ieee754_fmod(x,p+p); /* now x < 2p */ - if (((hx-hp)|(lx-lp))==0) return zero*x; - x = fabs(x); - p = fabs(p); - if (hp<0x00200000) { - if(x+x>p) { - x-=p; - if(x+x>=p) x -= p; - } - } else { - p_half = 0.5*p; - if(x>p_half) { - x-=p; - if(x>=p_half) x -= p; - } - } - GET_HIGH_WORD(hx,x); - if ((hx&0x7fffffff)==0) hx = 0; - SET_HIGH_WORD(x,hx^sx); - return x; + int q; + return remquo(x, y, &q); } #if LDBL_MANT_DIG == 53 __weak_reference(remainder, remainderl); +__weak_reference(remainder, drem); #endif Index: lib/msun/src/e_remainderf.c =================================================================== --- lib/msun/src/e_remainderf.c +++ lib/msun/src/e_remainderf.c @@ -1,65 +1,9 @@ -/* e_remainderf.c -- float version of e_remainder.c. - * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. - */ +#include -/* - * ==================================================== - * 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. - * ==================================================== - */ - -#include -__FBSDID("$FreeBSD$"); - -#include "math.h" -#include "math_private.h" - -static const float zero = 0.0; - - -float -__ieee754_remainderf(float x, float p) +float remainderf(float x, float y) { - int32_t hx,hp; - u_int32_t sx; - float p_half; - - GET_FLOAT_WORD(hx,x); - GET_FLOAT_WORD(hp,p); - sx = hx&0x80000000; - hp &= 0x7fffffff; - hx &= 0x7fffffff; - - /* purge off exception values */ - if((hp==0)|| /* p = 0 */ - (hx>=0x7f800000)|| /* x not finite */ - ((hp>0x7f800000))) /* p is NaN */ - return nan_mix_op(x, p, *)/nan_mix_op(x, p, *); - - - if (hp<=0x7effffff) x = __ieee754_fmodf(x,p+p); /* now x < 2p */ - if ((hx-hp)==0) return zero*x; - x = fabsf(x); - p = fabsf(p); - if (hp<0x01000000) { - if(x+x>p) { - x-=p; - if(x+x>=p) x -= p; - } - } else { - p_half = (float)0.5*p; - if(x>p_half) { - x-=p; - if(x>=p_half) x -= p; - } - } - GET_FLOAT_WORD(hx,x); - if ((hx&0x7fffffff)==0) hx = 0; - SET_FLOAT_WORD(x,hx^sx); - return x; + int q; + return remquof(x, y, &q); } + +__weak_reference(remainderf, dremf); Index: lib/msun/src/e_remainderl.c =================================================================== --- lib/msun/src/e_remainderl.c +++ lib/msun/src/e_remainderl.c @@ -1,40 +1,15 @@ -/*- - * SPDX-License-Identifier: BSD-2-Clause-FreeBSD - * - * Copyright (c) 2008 David Schultz - * All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * - * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND - * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE - * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE - * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS - * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) - * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT - * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY - * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF - * SUCH DAMAGE. - */ - -#include -__FBSDID("$FreeBSD$"); - #include +#include -long double -remainderl(long double x, long double y) +#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 +long double remainderl(long double x, long double y) { - int quo; - - return (remquol(x, y, &quo)); + return remainder(x, y); +} +#else +long double remainderl(long double x, long double y) +{ + int q; + return remquol(x, y, &q); } +#endif Index: lib/msun/src/libm.h =================================================================== --- /dev/null +++ lib/msun/src/libm.h @@ -0,0 +1,61 @@ +/* Compat header to allow using musl sources without local modifications. */ +#include +#include + +#include +#include +#include + +#include "math_private.h" + +#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 +#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 && __BYTE_ORDER == __LITTLE_ENDIAN +union ldshape { + long double f; + struct { + uint64_t m; + uint16_t se; + } i; +}; +#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 && __BYTE_ORDER == __BIG_ENDIAN +/* This is the m68k variant of 80-bit long double, and this definition only works + * on archs where the alignment requirement of uint64_t is <= 4. */ +union ldshape { + long double f; + struct { + uint16_t se; + uint16_t pad; + uint64_t m; + } i; +}; +#elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384 && __BYTE_ORDER == __LITTLE_ENDIAN +union ldshape { + long double f; + struct { + uint64_t lo; + uint32_t mid; + uint16_t top; + uint16_t se; + } i; + struct { + uint64_t lo; + uint64_t hi; + } i2; +}; +#elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384 && __BYTE_ORDER == __BIG_ENDIAN +union ldshape { + long double f; + struct { + uint16_t se; + uint16_t top; + uint32_t mid; + uint64_t lo; + } i; + struct { + uint64_t hi; + uint64_t lo; + } i2; +}; +#else +#error Unsupported long double representation +#endif Index: lib/msun/src/s_modf.c =================================================================== --- lib/msun/src/s_modf.c +++ lib/msun/src/s_modf.c @@ -1,79 +1,34 @@ -/* @(#)s_modf.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. - * ==================================================== - */ +#include "libm.h" -#ifndef lint -static char rcsid[] = "$FreeBSD$"; -#endif - -/* - * modf(double x, double *iptr) - * return fraction part of x, and return x's integral part in *iptr. - * Method: - * Bit twiddling. - * - * Exception: - * No exception. - */ - -#include "math.h" -#include "math_private.h" +double modf(double x, double *iptr) +{ + union {double f; uint64_t i;} u = {x}; + uint64_t mask; + int e = (int)(u.i>>52 & 0x7ff) - 0x3ff; -static const double one = 1.0; + /* no fractional part */ + if (e >= 52) { + *iptr = x; + if (e == 0x400 && u.i<<12 != 0) /* nan */ + return x; + u.i &= 1ULL<<63; + return u.f; + } -double -modf(double x, double *iptr) -{ - int32_t i0,i1,j0; - u_int32_t i; - EXTRACT_WORDS(i0,i1,x); - j0 = ((i0>>20)&0x7ff)-0x3ff; /* exponent of x */ - if(j0<20) { /* integer part in high x */ - if(j0<0) { /* |x|<1 */ - INSERT_WORDS(*iptr,i0&0x80000000,0); /* *iptr = +-0 */ + /* no integral part*/ + if (e < 0) { + u.i &= 1ULL<<63; + *iptr = u.f; return x; - } else { - i = (0x000fffff)>>j0; - if(((i0&i)|i1)==0) { /* x is integral */ - u_int32_t high; - *iptr = x; - GET_HIGH_WORD(high,x); - INSERT_WORDS(x,high&0x80000000,0); /* return +-0 */ - return x; - } else { - INSERT_WORDS(*iptr,i0&(~i),0); - return x - *iptr; - } - } - } else if (j0>51) { /* no fraction part */ - u_int32_t high; - if (j0 == 0x400) { /* inf/NaN */ - *iptr = x; - return 0.0 / x; - } - *iptr = x*one; - GET_HIGH_WORD(high,x); - INSERT_WORDS(x,high&0x80000000,0); /* return +-0 */ - return x; - } else { /* fraction part in low x */ - i = ((u_int32_t)(0xffffffff))>>(j0-20); - if((i1&i)==0) { /* x is integral */ - u_int32_t high; + } + + mask = -1ULL>>12>>e; + if ((u.i & mask) == 0) { *iptr = x; - GET_HIGH_WORD(high,x); - INSERT_WORDS(x,high&0x80000000,0); /* return +-0 */ - return x; - } else { - INSERT_WORDS(*iptr,i0,i1&(~i)); - return x - *iptr; - } + u.i &= 1ULL<<63; + return u.f; } + u.i &= ~mask; + *iptr = u.f; + return x - u.f; } Index: lib/msun/src/s_modff.c =================================================================== --- lib/msun/src/s_modff.c +++ lib/msun/src/s_modff.c @@ -1,57 +1,34 @@ -/* s_modff.c -- float version of s_modf.c. - * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. - */ +#include "libm.h" -/* - * ==================================================== - * 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. - * ==================================================== - */ - -#include -__FBSDID("$FreeBSD$"); - -#include "math.h" -#include "math_private.h" - -static const float one = 1.0; - -float -modff(float x, float *iptr) +float modff(float x, float *iptr) { - int32_t i0,j0; - u_int32_t i; - GET_FLOAT_WORD(i0,x); - j0 = ((i0>>23)&0xff)-0x7f; /* exponent of x */ - if(j0<23) { /* 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 */ - u_int32_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; + union {float f; uint32_t i;} u = {x}; + uint32_t mask; + int e = (int)(u.i>>23 & 0xff) - 0x7f; + + /* no fractional part */ + if (e >= 23) { + *iptr = x; + if (e == 0x80 && u.i<<9 != 0) { /* nan */ + return x; } - } - } else { /* no fraction part */ - u_int32_t ix; - *iptr = x*one; - if (x != x) /* NaN */ + u.i &= 0x80000000; + return u.f; + } + /* no integral part */ + if (e < 0) { + u.i &= 0x80000000; + *iptr = u.f; return x; - GET_FLOAT_WORD(ix,x); - SET_FLOAT_WORD(x,ix&0x80000000); /* return +-0 */ - return x; } + + mask = 0x007fffff>>e; + if ((u.i & mask) == 0) { + *iptr = x; + u.i &= 0x80000000; + return u.f; + } + u.i &= ~mask; + *iptr = u.f; + return x - u.f; } Index: lib/msun/src/s_remquo.c =================================================================== --- lib/msun/src/s_remquo.c +++ lib/msun/src/s_remquo.c @@ -1,159 +1,82 @@ -/* @(#)e_fmod.c 1.3 95/01/18 */ -/*- - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunSoft, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ +#include +#include -#include -__FBSDID("$FreeBSD$"); - -#include - -#include "math.h" -#include "math_private.h" - -static const double Zero[] = {0.0, -0.0,}; - -/* - * Return the IEEE remainder and set *quo to the last n bits of the - * quotient, rounded to the nearest integer. We choose n=31 because - * we wind up computing all the integer bits of the quotient anyway as - * a side-effect of computing the remainder by the shift and subtract - * method. In practice, this is far more bits than are needed to use - * remquo in reduction algorithms. - */ -double -remquo(double x, double y, int *quo) +double remquo(double x, double y, int *quo) { - int32_t n,hx,hy,hz,ix,iy,sx,i; - u_int32_t lx,ly,lz,q,sxy; - - EXTRACT_WORDS(hx,lx,x); - EXTRACT_WORDS(hy,ly,y); - sxy = (hx ^ hy) & 0x80000000; - sx = hx&0x80000000; /* sign of x */ - hx ^=sx; /* |x| */ - hy &= 0x7fffffff; /* |y| */ + union {double f; uint64_t i;} ux = {x}, uy = {y}; + int ex = ux.i>>52 & 0x7ff; + int ey = uy.i>>52 & 0x7ff; + int sx = ux.i>>63; + int sy = uy.i>>63; + uint32_t q; + uint64_t i; + uint64_t uxi = ux.i; - /* purge off exception values */ - if((hy|ly)==0||(hx>=0x7ff00000)|| /* y=0,or x not finite */ - ((hy|((ly|-ly)>>31))>0x7ff00000)) /* or y is NaN */ - return nan_mix_op(x, y, *)/nan_mix_op(x, y, *); - if(hx<=hy) { - if((hx>31]; /* |x|=|y| return x*0*/ - } - } - - /* 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; - - /* 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; + *quo = 0; + if (uy.i<<1 == 0 || isnan(y) || ex == 0x7ff) + return (x*y)/(x*y); + if (ux.i<<1 == 0) + return x; - /* set up {hx,lx}, {hy,ly} and align y to x */ - if(ix >= -1022) - hx = 0x00100000|(0x000fffff&hx); - else { /* subnormal x, shift x to normal */ - n = -1022-ix; - if(n<=31) { - hx = (hx<>(32-n)); - lx <<= n; - } else { - hx = lx<<(n-32); - lx = 0; - } + /* normalize x and y */ + if (!ex) { + for (i = uxi<<12; i>>63 == 0; ex--, i <<= 1); + uxi <<= -ex + 1; + } else { + uxi &= -1ULL >> 12; + uxi |= 1ULL << 52; } - if(iy >= -1022) - hy = 0x00100000|(0x000fffff&hy); - else { /* subnormal y, shift y to normal */ - n = -1022-iy; - if(n<=31) { - hy = (hy<>(32-n)); - ly <<= n; - } else { - hy = ly<<(n-32); - ly = 0; - } + if (!ey) { + for (i = uy.i<<12; i>>63 == 0; ey--, i <<= 1); + uy.i <<= -ey + 1; + } else { + uy.i &= -1ULL >> 12; + uy.i |= 1ULL << 52; } - /* fix point fmod */ - n = ix - iy; q = 0; - while(n--) { - hz=hx-hy;lz=lx-ly; if(lx>31); lx = lx+lx;} - else {hx = hz+hz+(lz>>31); lx = lz+lz; q++;} - q <<= 1; + if (ex < ey) { + if (ex+1 == ey) + goto end; + return x; } - hz=hx-hy;lz=lx-ly; if(lx=0) {hx=hz;lx=lz;q++;} - /* convert back to floating value and restore the sign */ - if((hx|lx)==0) { /* return sign(x)*0 */ - q &= 0x7fffffff; - *quo = (sxy ? -q : q); - return Zero[(u_int32_t)sx>>31]; + /* x mod y */ + for (; ex > ey; ex--) { + i = uxi - uy.i; + if (i >> 63 == 0) { + uxi = i; + q++; + } + uxi <<= 1; + q <<= 1; } - while(hx<0x00100000) { /* normalize x */ - hx = hx+hx+(lx>>31); lx = lx+lx; - iy -= 1; + i = uxi - uy.i; + if (i >> 63 == 0) { + uxi = i; + q++; } - if(iy>= -1022) { /* normalize output */ - hx = ((hx-0x00100000)|((iy+1023)<<20)); - } else { /* subnormal output */ - n = -1022 - iy; - if(n<=20) { - lx = (lx>>n)|((u_int32_t)hx<<(32-n)); - hx >>= n; - } else if (n<=31) { - lx = (hx<<(32-n))|(lx>>n); hx = 0; - } else { - lx = hx>>(n-32); hx = 0; - } + if (uxi == 0) + ex = -60; + else + for (; uxi>>52 == 0; uxi <<= 1, ex--); + end: + /* scale result and decide between |x| and |x|-|y| */ + if (ex > 0) { + uxi -= 1ULL << 52; + uxi |= (uint64_t)ex << 52; + } else { + uxi >>= -ex + 1; } -fixup: - INSERT_WORDS(x,hx,lx); - y = fabs(y); - if (y < 0x1p-1021) { - if (x+x>y || (x+x==y && (q & 1))) { + ux.i = uxi; + x = ux.f; + if (sy) + y = -y; + if (ex == ey || (ex+1 == ey && (2*x > y || (2*x == y && q%2)))) { + x -= y; q++; - x-=y; - } - } else if (x>0.5*y || (x==0.5*y && (q & 1))) { - q++; - x-=y; } - GET_HIGH_WORD(hx,x); - SET_HIGH_WORD(x,hx^sx); q &= 0x7fffffff; - *quo = (sxy ? -q : q); - return x; + *quo = sx^sy ? -(int)q : (int)q; + return sx ? -x : x; } - -#if LDBL_MANT_DIG == 53 -__weak_reference(remquo, remquol); -#endif Index: lib/msun/src/s_remquof.c =================================================================== --- lib/msun/src/s_remquof.c +++ lib/msun/src/s_remquof.c @@ -1,122 +1,82 @@ -/* @(#)e_fmod.c 1.3 95/01/18 */ -/*- - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunSoft, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ +#include +#include -#include -__FBSDID("$FreeBSD$"); - -#include "math.h" -#include "math_private.h" - -static const float Zero[] = {0.0, -0.0,}; - -/* - * Return the IEEE remainder and set *quo to the last n bits of the - * quotient, rounded to the nearest integer. We choose n=31 because - * we wind up computing all the integer bits of the quotient anyway as - * a side-effect of computing the remainder by the shift and subtract - * method. In practice, this is far more bits than are needed to use - * remquo in reduction algorithms. - */ -float -remquof(float x, float y, int *quo) +float remquof(float x, float y, int *quo) { - int32_t n,hx,hy,hz,ix,iy,sx,i; - u_int32_t q,sxy; - - GET_FLOAT_WORD(hx,x); - GET_FLOAT_WORD(hy,y); - sxy = (hx ^ hy) & 0x80000000; - sx = hx&0x80000000; /* sign of x */ - hx ^=sx; /* |x| */ - hy &= 0x7fffffff; /* |y| */ + union {float f; uint32_t i;} ux = {x}, uy = {y}; + int ex = ux.i>>23 & 0xff; + int ey = uy.i>>23 & 0xff; + int sx = ux.i>>31; + int sy = uy.i>>31; + uint32_t q; + uint32_t i; + uint32_t uxi = ux.i; - /* purge off exception values */ - if(hy==0||hx>=0x7f800000||hy>0x7f800000) /* y=0,NaN;or x not finite */ - return nan_mix_op(x, y, *)/nan_mix_op(x, y, *); - if(hx>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; + *quo = 0; + if (uy.i<<1 == 0 || isnan(y) || ex == 0xff) + return (x*y)/(x*y); + if (ux.i<<1 == 0) + return x; - /* 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; - - /* set up {hx,lx}, {hy,ly} and align y to x */ - if(ix >= -126) - hx = 0x00800000|(0x007fffff&hx); - else { /* subnormal x, shift x to normal */ - n = -126-ix; - hx <<= n; + /* normalize x and y */ + if (!ex) { + for (i = uxi<<9; i>>31 == 0; ex--, i <<= 1); + uxi <<= -ex + 1; + } else { + uxi &= -1U >> 9; + uxi |= 1U << 23; } - if(iy >= -126) - hy = 0x00800000|(0x007fffff&hy); - else { /* subnormal y, shift y to normal */ - n = -126-iy; - hy <<= n; + if (!ey) { + for (i = uy.i<<9; i>>31 == 0; ey--, i <<= 1); + uy.i <<= -ey + 1; + } else { + uy.i &= -1U >> 9; + uy.i |= 1U << 23; } - /* fix point fmod */ - n = ix - iy; q = 0; - while(n--) { - hz=hx-hy; - if(hz<0) hx = hx << 1; - else {hx = hz << 1; q++;} - q <<= 1; + if (ex < ey) { + if (ex+1 == ey) + goto end; + return x; } - hz=hx-hy; - if(hz>=0) {hx=hz;q++;} - /* convert back to floating value and restore the sign */ - if(hx==0) { /* return sign(x)*0 */ - q &= 0x7fffffff; - *quo = (sxy ? -q : q); - return Zero[(u_int32_t)sx>>31]; + /* x mod y */ + for (; ex > ey; ex--) { + i = uxi - uy.i; + if (i >> 31 == 0) { + uxi = i; + q++; + } + uxi <<= 1; + q <<= 1; } - while(hx<0x00800000) { /* normalize x */ - hx <<= 1; - iy -= 1; + i = uxi - uy.i; + if (i >> 31 == 0) { + uxi = i; + q++; } - if(iy>= -126) { /* normalize output */ - hx = ((hx-0x00800000)|((iy+127)<<23)); - } else { /* subnormal output */ - n = -126 - iy; - hx >>= n; + if (uxi == 0) + ex = -30; + else + for (; uxi>>23 == 0; uxi <<= 1, ex--); + end: + /* scale result and decide between |x| and |x|-|y| */ + if (ex > 0) { + uxi -= 1U << 23; + uxi |= (uint32_t)ex << 23; + } else { + uxi >>= -ex + 1; } -fixup: - SET_FLOAT_WORD(x,hx); - y = fabsf(y); - if (y < 0x1p-125f) { - if (x+x>y || (x+x==y && (q & 1))) { + ux.i = uxi; + x = ux.f; + if (sy) + y = -y; + if (ex == ey || (ex+1 == ey && (2*x > y || (2*x == y && q%2)))) { + x -= y; q++; - x-=y; - } - } else if (x>0.5f*y || (x==0.5f*y && (q & 1))) { - q++; - x-=y; } - GET_FLOAT_WORD(hx,x); - SET_FLOAT_WORD(x,hx^sx); q &= 0x7fffffff; - *quo = (sxy ? -q : q); - return x; + *quo = sx^sy ? -(int)q : (int)q; + return sx ? -x : x; }