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_synonyms.h" 44 #include "libm_protos.h" 45 #include <math.h> /* fabs() */ 46 #include <sys/isa_defs.h> 47 48 #if defined(_BIG_ENDIAN) 49 #define HIWORD 0 50 #define LOWORD 1 51 #else 52 #define HIWORD 1 53 #define LOWORD 0 54 #endif 55 #define __HI(x) ((int *) &x)[HIWORD] 56 #define __LO(x) ((int *) &x)[LOWORD] 57 58 static const double one = 1.0, Zero[] = {0.0, -0.0}; 59 60 static double 61 fmodquo(double x, double y, int *quo) { 62 int n, hx, hy, hz, ix, iy, sx, sq, i, m; 63 unsigned lx, ly, lz; 64 65 hx = __HI(x); /* high word of x */ 66 lx = __LO(x); /* low word of x */ 67 hy = __HI(y); /* high word of y */ 68 ly = __LO(y); /* low word of y */ 69 sx = hx & 0x80000000; /* sign of x */ 70 sq = (hx ^ hy) & 0x80000000; /* sign of x/y */ 71 hx ^= sx; /* |x| */ 72 hy &= 0x7fffffff; /* |y| */ 73 74 /* purge off exception values */ 75 *quo = 0; 76 if ((hy | ly) == 0 || hx >= 0x7ff00000 || /* y=0, or x !finite */ 77 (hy | ((ly | -ly) >> 31)) > 0x7ff00000) /* or y is NaN */ 78 return ((x * y) / (x * y)); 79 if (hx <= hy) { 80 if (hx < hy || lx < ly) 81 return (x); /* |x|<|y| return x */ 82 if (lx == ly) { 83 *quo = 1 + (sq >> 30); 84 /* |x|=|y| return x*0 */ 85 return (Zero[(unsigned) sx >> 31]); 86 } 87 } 88 89 /* determine ix = ilogb(x) */ 90 if (hx < 0x00100000) { /* subnormal x */ 91 if (hx == 0) { 92 for (ix = -1043, i = lx; i > 0; i <<= 1) 93 ix -= 1; 94 } else { 95 for (ix = -1022, i = (hx << 11); i > 0; i <<= 1) 96 ix -= 1; 97 } 98 } else 99 ix = (hx >> 20) - 1023; 100 101 /* determine iy = ilogb(y) */ 102 if (hy < 0x00100000) { /* subnormal y */ 103 if (hy == 0) { 104 for (iy = -1043, i = ly; i > 0; i <<= 1) 105 iy -= 1; 106 } else { 107 for (iy = -1022, i = (hy << 11); i > 0; i <<= 1) 108 iy -= 1; 109 } 110 } else 111 iy = (hy >> 20) - 1023; 112 113 /* set up {hx,lx}, {hy,ly} and align y to x */ 114 if (ix >= -1022) 115 hx = 0x00100000 | (0x000fffff & hx); 116 else { /* subnormal x, shift x to normal */ 117 n = -1022 - ix; 118 if (n <= 31) { 119 hx = (hx << n) | (lx >> (32 - n)); 120 lx <<= n; 121 } else { 122 hx = lx << (n - 32); 123 lx = 0; 124 } 125 } 126 if (iy >= -1022) 127 hy = 0x00100000 | (0x000fffff & hy); 128 else { /* subnormal y, shift y to normal */ 129 n = -1022 - iy; 130 if (n <= 31) { 131 hy = (hy << n) | (ly >> (32 - n)); 132 ly <<= n; 133 } else { 134 hy = ly << (n - 32); 135 ly = 0; 136 } 137 } 138 139 /* fix point fmod */ 140 n = ix - iy; 141 m = 0; 142 while (n--) { 143 hz = hx - hy; 144 lz = lx - ly; 145 if (lx < ly) 146 hz -= 1; 147 if (hz < 0) { 148 hx = hx + hx + (lx >> 31); 149 lx = lx + lx; 150 } else { 151 m += 1; 152 if ((hz | lz) == 0) { /* return sign(x)*0 */ 153 if (n < 31) 154 m <<= 1 + n; 155 else 156 m = 0; 157 m &= 0x7fffffff; 158 *quo = sq >= 0 ? m : -m; 159 return (Zero[(unsigned) sx >> 31]); 160 } 161 hx = hz + hz + (lz >> 31); 162 lx = lz + lz; 163 } 164 m += m; 165 } 166 hz = hx - hy; 167 lz = lx - ly; 168 if (lx < ly) 169 hz -= 1; 170 if (hz >= 0) { 171 hx = hz; 172 lx = lz; 173 m += 1; 174 } 175 m &= 0x7fffffff; 176 *quo = sq >= 0 ? m : -m; 177 178 /* convert back to floating value and restore the sign */ 179 if ((hx | lx) == 0) { /* return sign(x)*0 */ 180 return (Zero[(unsigned) sx >> 31]); 181 } 182 while (hx < 0x00100000) { /* normalize x */ 183 hx = hx + hx + (lx >> 31); 184 lx = lx + lx; 185 iy -= 1; 186 } 187 if (iy >= -1022) { /* normalize output */ 188 hx = (hx - 0x00100000) | ((iy + 1023) << 20); 189 __HI(x) = hx | sx; 190 __LO(x) = lx; 191 } else { /* subnormal output */ 192 n = -1022 - iy; 193 if (n <= 20) { 194 lx = (lx >> n) | ((unsigned) hx << (32 - n)); 195 hx >>= n; 196 } else if (n <= 31) { 197 lx = (hx << (32 - n)) | (lx >> n); 198 hx = sx; 199 } else { 200 lx = hx >> (n - 32); 201 hx = sx; 202 } 203 __HI(x) = hx | sx; 204 __LO(x) = lx; 205 x *= one; /* create necessary signal */ 206 } 207 return (x); /* exact output */ 208 } 209 210 #define zero Zero[0] 211 212 double 213 remquo(double x, double y, int *quo) { 214 int hx, hy, sx, sq; 215 double v; 216 unsigned ly; 217 218 hx = __HI(x); /* high word of x */ 219 hy = __HI(y); /* high word of y */ 220 ly = __LO(y); /* low word of y */ 221 sx = hx & 0x80000000; /* sign of x */ 222 sq = (hx ^ hy) & 0x80000000; /* sign of x/y */ 223 hx ^= sx; /* |x| */ 224 hy &= 0x7fffffff; /* |y| */ 225 226 /* purge off exception values */ 227 *quo = 0; 228 if ((hy | ly) == 0 || hx >= 0x7ff00000 || /* y=0, or x !finite */ 229 (hy | ((ly | -ly) >> 31)) > 0x7ff00000) /* or y is NaN */ 230 return ((x * y) / (x * y)); 231 232 y = fabs(y); 233 x = fabs(x); 234 if (hy <= 0x7fdfffff) { 235 x = fmodquo(x, y + y, quo); 236 *quo = ((*quo) & 0x3fffffff) << 1; 237 } 238 if (hy < 0x00200000) { 239 if (x + x > y) { 240 *quo += 1; 241 if (x == y) 242 x = zero; 243 else 244 x -= y; 245 if (x + x >= y) { 246 x -= y; 247 *quo += 1; 248 } 249 } 250 } else { 251 v = 0.5 * y; 252 if (x > v) { 253 *quo += 1; 254 if (x == y) 255 x = zero; 256 else 257 x -= y; 258 if (x >= v) { 259 x -= y; 260 *quo += 1; 261 } 262 } 263 } 264 if (sq != 0) 265 *quo = -(*quo); 266 return (sx == 0 ? x : -x); 267 } 268