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 __remquol = remquol 31 32 #include "libm.h" 33 #if defined(__SUNPRO_C) 34 #include <sunmath.h> /* fabsl */ 35 #endif 36 /* INDENT OFF */ 37 static const int 38 is = -0x7fffffff - 1, 39 im = 0x0000ffff, 40 iu = 0x00010000; 41 42 static const long double zero = 0.0L, one = 1.0L; 43 /* INDENT ON */ 44 45 #if defined(__sparc) 46 #define __H0(x) ((int *) &x)[0] 47 #define __H1(x) ((int *) &x)[1] 48 #define __H2(x) ((int *) &x)[2] 49 #define __H3(x) ((int *) &x)[3] 50 #else 51 #error Unsupported architecture 52 #endif 53 54 /* 55 * On entrance: *quo is initialized to 0, x finite and y non-zero & ordered 56 */ 57 static long double 58 fmodquol(long double x, long double y, int *quo) { 59 long double a, b; 60 int n, ix, iy, k, sx, sq, m; 61 int hx; 62 int x0, y0, z0, carry; 63 unsigned x1, x2, x3, y1, y2, y3, z1, z2, z3; 64 65 hx = __H0(x); 66 x1 = __H1(x); 67 x2 = __H2(x); 68 x3 = __H3(x); 69 y0 = __H0(y); 70 y1 = __H1(y); 71 y2 = __H2(y); 72 y3 = __H3(y); 73 74 sx = hx & is; 75 sq = (hx ^ y0) & is; 76 x0 = hx ^ sx; 77 y0 &= ~0x80000000; 78 79 a = fabsl(x); 80 b = fabsl(y); 81 if (a <= b) { 82 if (a < b) 83 return (x); 84 else { 85 *quo = 1 + (sq >> 30); 86 return (zero * x); 87 } 88 } 89 /* determine ix = ilogbl(x) */ 90 if (x0 < iu) { /* subnormal x */ 91 ix = 0; 92 ix = -16382; 93 while (x0 == 0) { 94 ix -= 16; 95 x0 = x1 >> 16; 96 x1 = (x1 << 16) | (x2 >> 16); 97 x2 = (x2 << 16) | (x3 >> 16); 98 x3 = (x3 << 16); 99 } 100 while (x0 < iu) { 101 ix -= 1; 102 x0 = (x0 << 1) | (x1 >> 31); 103 x1 = (x1 << 1) | (x2 >> 31); 104 x2 = (x2 << 1) | (x3 >> 31); 105 x3 <<= 1; 106 } 107 } else { 108 ix = (x0 >> 16) - 16383; 109 x0 = iu | (x0 & im); 110 } 111 112 /* determine iy = ilogbl(y) */ 113 if (y0 < iu) { /* subnormal y */ 114 iy = -16382; 115 while (y0 == 0) { 116 iy -= 16; 117 y0 = y1 >> 16; 118 y1 = (y1 << 16) | (y2 >> 16); 119 y2 = (y2 << 16) | (y3 >> 16); 120 y3 = (y3 << 16); 121 } 122 while (y0 < iu) { 123 iy -= 1; 124 y0 = (y0 << 1) | (y1 >> 31); 125 y1 = (y1 << 1) | (y2 >> 31); 126 y2 = (y2 << 1) | (y3 >> 31); 127 y3 <<= 1; 128 } 129 } else { 130 iy = (y0 >> 16) - 16383; 131 y0 = iu | (y0 & im); 132 } 133 134 135 /* fix point fmod */ 136 n = ix - iy; 137 m = 0; 138 while (n--) { 139 while (x0 == 0 && n >= 16) { 140 m <<= 16; 141 n -= 16; 142 x0 = x1 >> 16; 143 x1 = (x1 << 16) | (x2 >> 16); 144 x2 = (x2 << 16) | (x3 >> 16); 145 x3 = (x3 << 16); 146 } 147 while (x0 < iu && n >= 1) { 148 m += m; 149 n -= 1; 150 x0 = (x0 << 1) | (x1 >> 31); 151 x1 = (x1 << 1) | (x2 >> 31); 152 x2 = (x2 << 1) | (x3 >> 31); 153 x3 = (x3 << 1); 154 } 155 carry = 0; 156 z3 = x3 - y3; 157 carry = z3 > x3; 158 if (carry == 0) { 159 z2 = x2 - y2; 160 carry = z2 > x2; 161 } else { 162 z2 = x2 - y2 - 1; 163 carry = z2 >= x2; 164 } 165 if (carry == 0) { 166 z1 = x1 - y1; 167 carry = z1 > x1; 168 } else { 169 z1 = x1 - y1 - 1; 170 carry = z1 >= x1; 171 } 172 z0 = x0 - y0 - carry; 173 if (z0 < 0) { /* double x */ 174 x0 = x0 + x0 + ((x1 & is) != 0); 175 x1 = x1 + x1 + ((x2 & is) != 0); 176 x2 = x2 + x2 + ((x3 & is) != 0); 177 x3 = x3 + x3; 178 m += m; 179 } else { 180 m += 1; 181 if (z0 == 0) { 182 if ((z1 | z2 | z3) == 0) { 183 /* 0: we are done */ 184 if (n < 31) 185 m <<= (1 + n); 186 else 187 m = 0; 188 m &= ~0x80000000; 189 *quo = sq >= 0 ? m : -m; 190 __H0(a) = hx & is; 191 __H1(a) = __H2(a) = __H3(a) = 0; 192 return (a); 193 } 194 } 195 /* x = z << 1 */ 196 z0 = z0 + z0 + ((z1 & is) != 0); 197 z1 = z1 + z1 + ((z2 & is) != 0); 198 z2 = z2 + z2 + ((z3 & is) != 0); 199 z3 = z3 + z3; 200 x0 = z0; 201 x1 = z1; 202 x2 = z2; 203 x3 = z3; 204 m += m; 205 } 206 } 207 carry = 0; 208 z3 = x3 - y3; 209 carry = z3 > x3; 210 if (carry == 0) { 211 z2 = x2 - y2; 212 carry = z2 > x2; 213 } else { 214 z2 = x2 - y2 - 1; 215 carry = z2 >= x2; 216 } 217 if (carry == 0) { 218 z1 = x1 - y1; 219 carry = z1 > x1; 220 } else { 221 z1 = x1 - y1 - 1; 222 carry = z1 >= x1; 223 } 224 z0 = x0 - y0 - carry; 225 if (z0 >= 0) { 226 x0 = z0; 227 x1 = z1; 228 x2 = z2; 229 x3 = z3; 230 m += 1; 231 } 232 m &= ~0x80000000; 233 *quo = sq >= 0 ? m : -m; 234 235 /* convert back to floating value and restore the sign */ 236 if ((x0 | x1 | x2 | x3) == 0) { 237 __H0(a) = hx & is; 238 __H1(a) = __H2(a) = __H3(a) = 0; 239 return (a); 240 } 241 while (x0 < iu) { 242 if (x0 == 0) { 243 iy -= 16; 244 x0 = x1 >> 16; 245 x1 = (x1 << 16) | (x2 >> 16); 246 x2 = (x2 << 16) | (x3 >> 16); 247 x3 = (x3 << 16); 248 } else { 249 x0 = x0 + x0 + ((x1 & is) != 0); 250 x1 = x1 + x1 + ((x2 & is) != 0); 251 x2 = x2 + x2 + ((x3 & is) != 0); 252 x3 = x3 + x3; 253 iy -= 1; 254 } 255 } 256 257 /* normalize output */ 258 if (iy >= -16382) { 259 __H0(a) = sx | (x0 - iu) | ((iy + 16383) << 16); 260 __H1(a) = x1; 261 __H2(a) = x2; 262 __H3(a) = x3; 263 } else { /* subnormal output */ 264 n = -16382 - iy; 265 k = n & 31; 266 if (k <= 16) { 267 x3 = (x2 << (32 - k)) | (x3 >> k); 268 x2 = (x1 << (32 - k)) | (x2 >> k); 269 x1 = (x0 << (32 - k)) | (x1 >> k); 270 x0 >>= k; 271 } else { 272 x3 = (x2 << (32 - k)) | (x3 >> k); 273 x2 = (x1 << (32 - k)) | (x2 >> k); 274 x1 = (x0 << (32 - k)) | (x1 >> k); 275 x0 = 0; 276 } 277 while (n >= 32) { 278 n -= 32; 279 x3 = x2; 280 x2 = x1; 281 x1 = x0; 282 x0 = 0; 283 } 284 __H0(a) = x0 | sx; 285 __H1(a) = x1; 286 __H2(a) = x2; 287 __H3(a) = x3; 288 a *= one; 289 } 290 return (a); 291 } 292 293 long double 294 remquol(long double x, long double y, int *quo) { 295 int hx, hy, sx, sq; 296 long double v; 297 298 hx = __H0(x); /* high word of x */ 299 hy = __H0(y); /* high word of y */ 300 sx = hx & is; /* sign of x */ 301 sq = (hx ^ hy) & is; /* sign of x/y */ 302 hx ^= sx; /* |x| */ 303 hy &= ~0x80000000; 304 305 /* purge off exception values */ 306 *quo = 0; 307 /* y=0, y is NaN, x is NaN or inf */ 308 if (y == 0.0L || y != y || hx >= 0x7fff0000) 309 return ((x * y) / (x * y)); 310 311 y = fabsl(y); 312 x = fabsl(x); 313 if (hy <= 0x7ffdffff) { 314 x = fmodquol(x, y + y, quo); 315 *quo = ((*quo) & 0x3fffffff) << 1; 316 } 317 if (hy < 0x00020000) { 318 if (x + x > y) { 319 *quo += 1; 320 if (x == y) 321 x = zero; 322 else 323 x -= y; 324 if (x + x >= y) { 325 x -= y; 326 *quo += 1; 327 } 328 } 329 } else { 330 v = 0.5L * y; 331 if (x > v) { 332 *quo += 1; 333 if (x == y) 334 x = zero; 335 else 336 x -= y; 337 if (x >= v) { 338 x -= y; 339 *quo += 1; 340 } 341 } 342 } 343 if (sq != 0) 344 *quo = -(*quo); 345 return (sx == 0 ? x : -x); 346 } 347