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 #if defined(ELFOBJ) 31 #pragma weak llrintl = __llrintl 32 #if defined(__sparcv9) || defined(__amd64) 33 #pragma weak lrintl = __llrintl 34 #pragma weak __lrintl = __llrintl 35 #endif 36 #endif 37 38 #include "libm.h" 39 40 #if defined(__sparc) 41 42 #include "fma.h" 43 #include "fenv_inlines.h" 44 45 long long 46 llrintl(long double x) { 47 union { 48 unsigned i[4]; 49 long double q; 50 } xx; 51 union { 52 unsigned i[2]; 53 long long l; 54 } zz; 55 union { 56 unsigned i; 57 float f; 58 } tt; 59 unsigned int hx, sx, frac, fsr; 60 int rm, j; 61 volatile float dummy; 62 63 xx.q = x; 64 sx = xx.i[0] & 0x80000000; 65 hx = xx.i[0] & ~0x80000000; 66 67 /* handle trivial cases */ 68 if (hx > 0x403e0000) { /* |x| > 2^63 + ... or x is nan */ 69 /* convert an out-of-range float */ 70 tt.i = sx | 0x7f000000; 71 return ((long long) tt.f); 72 } else if ((hx | xx.i[1] | xx.i[2] | xx.i[3]) == 0) /* x is zero */ 73 return (0LL); 74 75 /* get the rounding mode */ 76 __fenv_getfsr32(&fsr); 77 rm = fsr >> 30; 78 79 /* flip the sense of directed roundings if x is negative */ 80 if (sx) 81 rm ^= rm >> 1; 82 83 /* handle |x| < 1 */ 84 if (hx < 0x3fff0000) { 85 dummy = 1.0e30f; /* x is nonzero, so raise inexact */ 86 dummy += 1.0e-30f; 87 if (rm == FSR_RP || (rm == FSR_RN && (hx >= 0x3ffe0000 && 88 ((hx & 0xffff) | xx.i[1] | xx.i[2] | xx.i[3])))) 89 return (sx ? -1LL : 1LL); 90 return (0LL); 91 } 92 93 /* extract the integer and fractional parts of x */ 94 j = 0x406f - (hx >> 16); 95 xx.i[0] = 0x10000 | (xx.i[0] & 0xffff); 96 if (j >= 96) { 97 zz.i[0] = 0; 98 zz.i[1] = xx.i[0] >> (j - 96); 99 frac = ((xx.i[0] << 1) << (127 - j)) | (xx.i[1] >> (j - 96)); 100 if (((xx.i[1] << 1) << (127 - j)) | xx.i[2] | xx.i[3]) 101 frac |= 1; 102 } else if (j >= 64) { 103 zz.i[0] = xx.i[0] >> (j - 64); 104 zz.i[1] = ((xx.i[0] << 1) << (95 - j)) | (xx.i[1] >> (j - 64)); 105 frac = ((xx.i[1] << 1) << (95 - j)) | (xx.i[2] >> (j - 64)); 106 if (((xx.i[2] << 1) << (95 - j)) | xx.i[3]) 107 frac |= 1; 108 } else { 109 zz.i[0] = ((xx.i[0] << 1) << (63 - j)) | (xx.i[1] >> (j - 32)); 110 zz.i[1] = ((xx.i[1] << 1) << (63 - j)) | (xx.i[2] >> (j - 32)); 111 frac = ((xx.i[2] << 1) << (63 - j)) | (xx.i[3] >> (j - 32)); 112 if ((xx.i[3] << 1) << (63 - j)) 113 frac |= 1; 114 } 115 116 /* round */ 117 if (frac && (rm == FSR_RP || (rm == FSR_RN && (frac > 0x80000000u || 118 (frac == 0x80000000 && (zz.i[1] & 1)))))) { 119 if (++zz.i[1] == 0) 120 zz.i[0]++; 121 } 122 123 /* check for result out of range (note that z is |x| at this point) */ 124 if (zz.i[0] > 0x80000000u || (zz.i[0] == 0x80000000 && (zz.i[1] || 125 !sx))) { 126 tt.i = sx | 0x7f000000; 127 return ((long long) tt.f); 128 } 129 130 /* raise inexact if need be */ 131 if (frac) { 132 dummy = 1.0e30F; 133 dummy += 1.0e-30F; 134 } 135 136 /* negate result if need be */ 137 if (sx) { 138 zz.i[0] = ~zz.i[0]; 139 zz.i[1] = -zz.i[1]; 140 if (zz.i[1] == 0) 141 zz.i[0]++; 142 } 143 return (zz.l); 144 } 145 #elif defined(__x86) 146 long long 147 llrintl(long double x) { 148 /* 149 * Note: The following code works on x86 (in the default rounding 150 * precision mode), but one ought to just use the fistpll instruction 151 * instead. 152 */ 153 union { 154 unsigned i[3]; 155 long double e; 156 } xx, yy; 157 int ex; 158 159 xx.e = x; 160 ex = xx.i[2] & 0x7fff; 161 162 if (ex < 0x403e) { /* |x| < 2^63 */ 163 /* add and subtract a power of two to round x to an integer */ 164 yy.i[2] = (xx.i[2] & 0x8000) | 0x403e; 165 yy.i[1] = 0x80000000; 166 yy.i[0] = 0; 167 x = (x + yy.e) - yy.e; 168 } 169 170 /* now x is nan, inf, or integral */ 171 return ((long long) x); 172 } 173 #else 174 #error Unknown architecture 175 #endif 176