diff options
Diffstat (limited to 'sysdeps/ieee754/ldbl-128ibm/e_fmodl.c')
| -rw-r--r-- | sysdeps/ieee754/ldbl-128ibm/e_fmodl.c | 149 |
1 files changed, 0 insertions, 149 deletions
diff --git a/sysdeps/ieee754/ldbl-128ibm/e_fmodl.c b/sysdeps/ieee754/ldbl-128ibm/e_fmodl.c deleted file mode 100644 index 5284fd0fd5..0000000000 --- a/sysdeps/ieee754/ldbl-128ibm/e_fmodl.c +++ /dev/null @@ -1,149 +0,0 @@ -/* e_fmodl.c -- long double version of e_fmod.c. - * Conversion to IEEE quad long double by Jakub Jelinek, jj@ultra.linux.cz. - */ -/* - * ==================================================== - * 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. - * ==================================================== - */ - -/* - * __ieee754_fmodl(x,y) - * Return x mod y in exact arithmetic - * Method: shift and subtract - */ - -#include <math.h> -#include <math_private.h> -#include <ieee754.h> - -static const long double one = 1.0, Zero[] = {0.0, -0.0,}; - -long double -__ieee754_fmodl (long double x, long double y) -{ - int64_t hx, hy, hz, sx, sy; - uint64_t lx, ly, lz; - int n, ix, iy; - double xhi, xlo, yhi, ylo; - - ldbl_unpack (x, &xhi, &xlo); - EXTRACT_WORDS64 (hx, xhi); - EXTRACT_WORDS64 (lx, xlo); - ldbl_unpack (y, &yhi, &ylo); - EXTRACT_WORDS64 (hy, yhi); - EXTRACT_WORDS64 (ly, ylo); - sx = hx&0x8000000000000000ULL; /* sign of x */ - hx ^= sx; /* |x| */ - sy = hy&0x8000000000000000ULL; /* sign of y */ - hy ^= sy; /* |y| */ - - /* purge off exception values */ - if(__builtin_expect(hy==0 || - (hx>=0x7ff0000000000000LL)|| /* y=0,or x not finite */ - (hy>0x7ff0000000000000LL),0)) /* or y is NaN */ - return (x*y)/(x*y); - if (__glibc_unlikely (hx <= hy)) - { - /* If |x| < |y| return x. */ - if (hx < hy) - return x; - /* At this point the absolute value of the high doubles of - x and y must be equal. */ - if ((lx & 0x7fffffffffffffffLL) == 0 - && (ly & 0x7fffffffffffffffLL) == 0) - /* Both low parts are zero. The result should be an - appropriately signed zero, but the subsequent logic - could treat them as unequal, depending on the signs - of the low parts. */ - return Zero[(uint64_t) sx >> 63]; - /* If the low double of y is the same sign as the high - double of y (ie. the low double increases |y|)... */ - if (((ly ^ sy) & 0x8000000000000000LL) == 0 - /* ... then a different sign low double to high double - for x or same sign but lower magnitude... */ - && (int64_t) (lx ^ sx) < (int64_t) (ly ^ sy)) - /* ... means |x| < |y|. */ - return x; - /* If the low double of x differs in sign to the high - double of x (ie. the low double decreases |x|)... */ - if (((lx ^ sx) & 0x8000000000000000LL) != 0 - /* ... then a different sign low double to high double - for y with lower magnitude (we've already caught - the same sign for y case above)... */ - && (int64_t) (lx ^ sx) > (int64_t) (ly ^ sy)) - /* ... means |x| < |y|. */ - return x; - /* If |x| == |y| return x*0. */ - if ((lx ^ sx) == (ly ^ sy)) - return Zero[(uint64_t) sx >> 63]; - } - - /* Make the IBM extended format 105 bit mantissa look like the ieee854 112 - bit mantissa so the following operations will give the correct - result. */ - ldbl_extract_mantissa(&hx, &lx, &ix, x); - ldbl_extract_mantissa(&hy, &ly, &iy, y); - - if (__glibc_unlikely (ix == -IEEE754_DOUBLE_BIAS)) - { - /* subnormal x, shift x to normal. */ - while ((hx & (1LL << 48)) == 0) - { - hx = (hx << 1) | (lx >> 63); - lx = lx << 1; - ix -= 1; - } - } - - if (__glibc_unlikely (iy == -IEEE754_DOUBLE_BIAS)) - { - /* subnormal y, shift y to normal. */ - while ((hy & (1LL << 48)) == 0) - { - hy = (hy << 1) | (ly >> 63); - ly = ly << 1; - iy -= 1; - } - } - - /* fix point fmod */ - n = ix - iy; - while(n--) { - hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1; - if(hz<0){hx = hx+hx+(lx>>63); lx = lx+lx;} - else { - if((hz|lz)==0) /* return sign(x)*0 */ - return Zero[(u_int64_t)sx>>63]; - hx = hz+hz+(lz>>63); lx = lz+lz; - } - } - hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1; - if(hz>=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_int64_t)sx>>63]; - while(hx<0x0001000000000000LL) { /* normalize x */ - hx = hx+hx+(lx>>63); lx = lx+lx; - iy -= 1; - } - if(__builtin_expect(iy>= -1022,0)) { /* normalize output */ - x = ldbl_insert_mantissa((sx>>63), iy, hx, lx); - } else { /* subnormal output */ - n = -1022 - iy; - /* We know 1 <= N <= 52, and that there are no nonzero - bits in places below 2^-1074. */ - lx = (lx >> n) | ((u_int64_t) hx << (64 - n)); - hx >>= n; - x = ldbl_insert_mantissa((sx>>63), -1023, hx, lx); - x *= one; /* create necessary signal */ - } - return x; /* exact output */ -} -strong_alias (__ieee754_fmodl, __fmodl_finite) |