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