1 /* 2 * CDDL HEADER START 3 * 4 * The contents of this file are subject to the terms of the 5 * Common Development and Distribution License (the "License"). 6 * You may not use this file except in compliance with the License. 7 * 8 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE 9 * or http://www.opensolaris.org/os/licensing. 10 * See the License for the specific language governing permissions 11 * and limitations under the License. 12 * 13 * When distributing Covered Code, include this CDDL HEADER in each 14 * file and include the License file at usr/src/OPENSOLARIS.LICENSE. 15 * If applicable, add the following below this CDDL HEADER, with the 16 * fields enclosed by brackets "[]" replaced with your own identifying 17 * information: Portions Copyright [yyyy] [name of copyright owner] 18 * 19 * CDDL HEADER END 20 */ 21 22 /* 23 * Copyright 2011 Nexenta Systems, Inc. All rights reserved. 24 */ 25 /* 26 * Copyright 2006 Sun Microsystems, Inc. All rights reserved. 27 * Use is subject to license terms. 28 */ 29 30 #pragma weak __remquo = remquo 31 32 /* INDENT OFF */ 33 /* 34 * double remquo(double x, double y, int *quo) return remainder(x,y) and an 35 * integer pointer quo such that *quo = N mod {2**31}, where N is the 36 * exact integral part of x/y rounded to nearest even. 37 * 38 * remquo call internal fmodquo 39 */ 40 /* INDENT ON */ 41 42 #include "libm.h" 43 #include "libm_protos.h" 44 #include <math.h> /* fabs() */ 45 #include <sys/isa_defs.h> 46 47 #if defined(_BIG_ENDIAN) 48 #define HIWORD 0 49 #define LOWORD 1 50 #else 51 #define HIWORD 1 52 #define LOWORD 0 53 #endif 54 #define __HI(x) ((int *) &x)[HIWORD] 55 #define __LO(x) ((int *) &x)[LOWORD] 56 57 static const double one = 1.0, Zero[] = {0.0, -0.0}; 58 59 static double 60 fmodquo(double x, double y, int *quo) { 61 int n, hx, hy, hz, ix, iy, sx, sq, i, m; 62 unsigned lx, ly, lz; 63 64 hx = __HI(x); /* high word of x */ 65 lx = __LO(x); /* low word of x */ 66 hy = __HI(y); /* high word of y */ 67 ly = __LO(y); /* low word of y */ 68 sx = hx & 0x80000000; /* sign of x */ 69 sq = (hx ^ hy) & 0x80000000; /* sign of x/y */ 70 hx ^= sx; /* |x| */ 71 hy &= 0x7fffffff; /* |y| */ 72 73 /* purge off exception values */ 74 *quo = 0; 75 if ((hy | ly) == 0 || hx >= 0x7ff00000 || /* y=0, or x !finite */ 76 (hy | ((ly | -ly) >> 31)) > 0x7ff00000) /* or y is NaN */ 77 return ((x * y) / (x * y)); 78 if (hx <= hy) { 79 if (hx < hy || lx < ly) 80 return (x); /* |x|<|y| return x */ 81 if (lx == ly) { 82 *quo = 1 + (sq >> 30); 83 /* |x|=|y| return x*0 */ 84 return (Zero[(unsigned) sx >> 31]); 85 } 86 } 87 88 /* determine ix = ilogb(x) */ 89 if (hx < 0x00100000) { /* subnormal x */ 90 if (hx == 0) { 91 for (ix = -1043, i = lx; i > 0; i <<= 1) 92 ix -= 1; 93 } else { 94 for (ix = -1022, i = (hx << 11); i > 0; i <<= 1) 95 ix -= 1; 96 } 97 } else 98 ix = (hx >> 20) - 1023; 99 100 /* determine iy = ilogb(y) */ 101 if (hy < 0x00100000) { /* subnormal y */ 102 if (hy == 0) { 103 for (iy = -1043, i = ly; i > 0; i <<= 1) 104 iy -= 1; 105 } else { 106 for (iy = -1022, i = (hy << 11); i > 0; i <<= 1) 107 iy -= 1; 108 } 109 } else 110 iy = (hy >> 20) - 1023; 111 112 /* set up {hx,lx}, {hy,ly} and align y to x */ 113 if (ix >= -1022) 114 hx = 0x00100000 | (0x000fffff & hx); 115 else { /* subnormal x, shift x to normal */ 116 n = -1022 - ix; 117 if (n <= 31) { 118 hx = (hx << n) | (lx >> (32 - n)); 119 lx <<= n; 120 } else { 121 hx = lx << (n - 32); 122 lx = 0; 123 } 124 } 125 if (iy >= -1022) 126 hy = 0x00100000 | (0x000fffff & hy); 127 else { /* subnormal y, shift y to normal */ 128 n = -1022 - iy; 129 if (n <= 31) { 130 hy = (hy << n) | (ly >> (32 - n)); 131 ly <<= n; 132 } else { 133 hy = ly << (n - 32); 134 ly = 0; 135 } 136 } 137 138 /* fix point fmod */ 139 n = ix - iy; 140 m = 0; 141 while (n--) { 142 hz = hx - hy; 143 lz = lx - ly; 144 if (lx < ly) 145 hz -= 1; 146 if (hz < 0) { 147 hx = hx + hx + (lx >> 31); 148 lx = lx + lx; 149 } else { 150 m += 1; 151 if ((hz | lz) == 0) { /* return sign(x)*0 */ 152 if (n < 31) 153 m <<= 1 + n; 154 else 155 m = 0; 156 m &= 0x7fffffff; 157 *quo = sq >= 0 ? m : -m; 158 return (Zero[(unsigned) sx >> 31]); 159 } 160 hx = hz + hz + (lz >> 31); 161 lx = lz + lz; 162 } 163 m += m; 164 } 165 hz = hx - hy; 166 lz = lx - ly; 167 if (lx < ly) 168 hz -= 1; 169 if (hz >= 0) { 170 hx = hz; 171 lx = lz; 172 m += 1; 173 } 174 m &= 0x7fffffff; 175 *quo = sq >= 0 ? m : -m; 176 177 /* convert back to floating value and restore the sign */ 178 if ((hx | lx) == 0) { /* return sign(x)*0 */ 179 return (Zero[(unsigned) sx >> 31]); 180 } 181 while (hx < 0x00100000) { /* normalize x */ 182 hx = hx + hx + (lx >> 31); 183 lx = lx + lx; 184 iy -= 1; 185 } 186 if (iy >= -1022) { /* normalize output */ 187 hx = (hx - 0x00100000) | ((iy + 1023) << 20); 188 __HI(x) = hx | sx; 189 __LO(x) = lx; 190 } else { /* subnormal output */ 191 n = -1022 - iy; 192 if (n <= 20) { 193 lx = (lx >> n) | ((unsigned) hx << (32 - n)); 194 hx >>= n; 195 } else if (n <= 31) { 196 lx = (hx << (32 - n)) | (lx >> n); 197 hx = sx; 198 } else { 199 lx = hx >> (n - 32); 200 hx = sx; 201 } 202 __HI(x) = hx | sx; 203 __LO(x) = lx; 204 x *= one; /* create necessary signal */ 205 } 206 return (x); /* exact output */ 207 } 208 209 #define zero Zero[0] 210 211 double 212 remquo(double x, double y, int *quo) { 213 int hx, hy, sx, sq; 214 double v; 215 unsigned ly; 216 217 hx = __HI(x); /* high word of x */ 218 hy = __HI(y); /* high word of y */ 219 ly = __LO(y); /* low word of y */ 220 sx = hx & 0x80000000; /* sign of x */ 221 sq = (hx ^ hy) & 0x80000000; /* sign of x/y */ 222 hx ^= sx; /* |x| */ 223 hy &= 0x7fffffff; /* |y| */ 224 225 /* purge off exception values */ 226 *quo = 0; 227 if ((hy | ly) == 0 || hx >= 0x7ff00000 || /* y=0, or x !finite */ 228 (hy | ((ly | -ly) >> 31)) > 0x7ff00000) /* or y is NaN */ 229 return ((x * y) / (x * y)); 230 231 y = fabs(y); 232 x = fabs(x); 233 if (hy <= 0x7fdfffff) { 234 x = fmodquo(x, y + y, quo); 235 *quo = ((*quo) & 0x3fffffff) << 1; 236 } 237 if (hy < 0x00200000) { 238 if (x + x > y) { 239 *quo += 1; 240 if (x == y) 241 x = zero; 242 else 243 x -= y; 244 if (x + x >= y) { 245 x -= y; 246 *quo += 1; 247 } 248 } 249 } else { 250 v = 0.5 * y; 251 if (x > v) { 252 *quo += 1; 253 if (x == y) 254 x = zero; 255 else 256 x -= y; 257 if (x >= v) { 258 x -= y; 259 *quo += 1; 260 } 261 } 262 } 263 if (sq != 0) 264 *quo = -(*quo); 265 return (sx == 0 ? x : -x); 266 } 267