1 /* @(#)e_fmod.c 1.3 95/01/18 */ 2 /*- 3 * ==================================================== 4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 5 * 6 * Developed at SunSoft, a Sun Microsystems, Inc. business. 7 * Permission to use, copy, modify, and distribute this 8 * software is freely granted, provided that this notice 9 * is preserved. 10 * ==================================================== 11 */ 12 13 #include <sys/types.h> 14 #include <machine/ieee.h> 15 16 #include <float.h> 17 #include <math.h> 18 #include <stdint.h> 19 20 #include "math_private.h" 21 22 #define BIAS (LDBL_MAX_EXP - 1) 23 24 /* 25 * These macros add and remove an explicit integer bit in front of the 26 * fractional mantissa, if the architecture doesn't have such a bit by 27 * default already. 28 */ 29 #ifdef LDBL_IMPLICIT_NBIT 30 #define LDBL_NBIT 0 31 #define SET_NBIT(hx) ((hx) | (1ULL << LDBL_MANH_SIZE)) 32 #define HFRAC_BITS (EXT_FRACHBITS + EXT_FRACHMBITS) 33 #else 34 #define LDBL_NBIT 0x80000000 35 #define SET_NBIT(hx) (hx) 36 #define HFRAC_BITS (EXT_FRACHBITS + EXT_FRACHMBITS - 1) 37 #endif 38 39 #define MANL_SHIFT (EXT_FRACLMBITS + EXT_FRACLBITS - 1) 40 41 static const long double Zero[] = {0.0L, -0.0L}; 42 43 /* 44 * Return the IEEE remainder and set *quo to the last n bits of the 45 * quotient, rounded to the nearest integer. We choose n=31 because 46 * we wind up computing all the integer bits of the quotient anyway as 47 * a side-effect of computing the remainder by the shift and subtract 48 * method. In practice, this is far more bits than are needed to use 49 * remquo in reduction algorithms. 50 * 51 * Assumptions: 52 * - The low part of the mantissa fits in a manl_t exactly. 53 * - The high part of the mantissa fits in an int64_t with enough room 54 * for an explicit integer bit in front of the fractional bits. 55 */ 56 long double 57 remquol(long double x, long double y, int *quo) 58 { 59 int64_t hx,hz,hy,_hx; 60 uint64_t lx,ly,lz; 61 uint64_t sx,sxy; 62 int ix,iy,n,q; 63 64 GET_LDOUBLE_WORDS64(hx,lx,x); 65 GET_LDOUBLE_WORDS64(hy,ly,y); 66 sx = (hx>>48)&0x8000; 67 sxy = sx ^ ((hy>>48)&0x8000); 68 hx &= 0x7fffffffffffffffLL; /* |x| */ 69 hy &= 0x7fffffffffffffffLL; /* |y| */ 70 SET_LDOUBLE_WORDS64(x,hx,lx); 71 SET_LDOUBLE_WORDS64(y,hy,ly); 72 73 /* purge off exception values */ 74 if((hy|ly)==0 || /* y=0 */ 75 ((hx>>48) == BIAS + LDBL_MAX_EXP) || /* or x not finite */ 76 ((hy>>48) == BIAS + LDBL_MAX_EXP && 77 (((hy&0x0000ffffffffffffLL)&~LDBL_NBIT)|ly)!=0)) /* or y is NaN */ 78 return (x*y)/(x*y); 79 if((hx>>48)<=(hy>>48)) { 80 if(((hx>>48)<(hy>>48)) || 81 ((hx&0x0000ffffffffffffLL)<=(hy&0x0000ffffffffffffLL) && 82 ((hx&0x0000ffffffffffffLL)<(hy&0x0000ffffffffffffLL) || 83 lx<ly))) { 84 q = 0; 85 goto fixup; /* |x|<|y| return x or x-y */ 86 } 87 if((hx&0x0000ffffffffffffLL)==(hy&0x0000ffffffffffffLL) && 88 lx==ly) { 89 *quo = 1; 90 return Zero[sx!=0]; /* |x|=|y| return x*0*/ 91 } 92 } 93 94 /* determine ix = ilogb(x) */ 95 if((hx>>48) == 0) { /* subnormal x */ 96 x *= 0x1.0p512; 97 GET_LDOUBLE_WORDS64(hx,lx,x); 98 ix = (hx>>48) - (BIAS + 512); 99 } else { 100 ix = (hx>>48) - BIAS; 101 } 102 103 /* determine iy = ilogb(y) */ 104 if((hy>>48) == 0) { /* subnormal y */ 105 y *= 0x1.0p512; 106 GET_LDOUBLE_WORDS64(hy,ly,y); 107 iy = (hy>>48) - (BIAS + 512); 108 } else { 109 iy = (hy>>48) - BIAS; 110 } 111 112 /* set up {hx,lx}, {hy,ly} and align y to x */ 113 _hx = SET_NBIT(hx) & 0x0000ffffffffffffLL; 114 hy = SET_NBIT(hy); 115 116 /* fix point fmod */ 117 n = ix - iy; 118 q = 0; 119 120 while(n--) { 121 hz=_hx-hy;lz=lx-ly; if(lx<ly) hz -= 1; 122 if(hz<0){_hx = _hx+_hx+(lx>>MANL_SHIFT); lx = lx+lx;} 123 else {_hx = hz+hz+(lz>>MANL_SHIFT); lx = lz+lz; q++;} 124 q <<= 1; 125 } 126 hz=_hx-hy;lz=lx-ly; if(lx<ly) hz -= 1; 127 if(hz>=0) {_hx=hz;lx=lz;q++;} 128 129 /* convert back to floating value and restore the sign */ 130 if((_hx|lx)==0) { /* return sign(x)*0 */ 131 *quo = (sxy ? -q : q); 132 return Zero[sx!=0]; 133 } 134 while(_hx<(1ULL<<HFRAC_BITS)) { /* normalize x */ 135 _hx = _hx+_hx+(lx>>MANL_SHIFT); lx = lx+lx; 136 iy -= 1; 137 } 138 hx = (hx&0xffff000000000000LL) | (_hx&0x0000ffffffffffffLL); 139 if (iy < LDBL_MIN_EXP) { 140 hx = (hx&0x0000ffffffffffffLL) | (uint64_t)(iy + BIAS + 512)<<48; 141 SET_LDOUBLE_WORDS64(x,hx,lx); 142 x *= 0x1p-512; 143 GET_LDOUBLE_WORDS64(hx,lx,x); 144 } else { 145 hx = (hx&0x0000ffffffffffffLL) | (uint64_t)(iy + BIAS)<<48; 146 } 147 hx &= 0x7fffffffffffffffLL; 148 SET_LDOUBLE_WORDS64(x,hx,lx); 149 fixup: 150 y = fabsl(y); 151 if (y < LDBL_MIN * 2) { 152 if (x+x>y || (x+x==y && (q & 1))) { 153 q++; 154 x-=y; 155 } 156 } else if (x>0.5*y || (x==0.5*y && (q & 1))) { 157 q++; 158 x-=y; 159 } 160 161 GET_LDOUBLE_MSW64(hx,x); 162 hx ^= sx; 163 SET_LDOUBLE_MSW64(x,hx); 164 165 q &= 0x7fffffff; 166 *quo = (sxy ? -q : q); 167 return x; 168 } 169