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