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