xref: /illumos-gate/usr/src/lib/libm/common/m9x/nearbyintf.c (revision c124a83e09115de88ecccd4f689983f42a1d53bd)
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 nearbyintf = __nearbyintf
31 
32 #include "libm.h"
33 #include <fenv.h>
34 
35 float
36 __nearbyintf(float x) {
37 	union {
38 		unsigned i;
39 		float f;
40 	} xx;
41 	unsigned hx, sx, i, frac;
42 	int rm;
43 
44 	xx.f = x;
45 	sx = xx.i & 0x80000000;
46 	hx = xx.i & ~0x80000000;
47 
48 	/* handle trivial cases */
49 	if (hx >= 0x4b000000) {	/* x is nan, inf, or already integral */
50 		if (hx > 0x7f800000)	/* x is nan */
51 			return (x * x);		/* + -> * for Cheetah */
52 		return (x);
53 	} else if (hx == 0)		/* x is zero */
54 		return (x);
55 
56 	/* get the rounding mode */
57 	rm = fegetround();
58 
59 	/* flip the sense of directed roundings if x is negative */
60 	if (sx && (rm == FE_UPWARD || rm == FE_DOWNWARD))
61 		rm = (FE_UPWARD + FE_DOWNWARD) - rm;
62 
63 	/* handle |x| < 1 */
64 	if (hx < 0x3f800000) {
65 		if (rm == FE_UPWARD || (rm == FE_TONEAREST && hx > 0x3f000000))
66 			xx.i = sx | 0x3f800000;
67 		else
68 			xx.i = sx;
69 		return (xx.f);
70 	}
71 
72 	/* round x at the integer bit */
73 	i = 1 << (0x96 - (hx >> 23));
74 	frac = hx & (i - 1);
75 	if (!frac)
76 		return (x);
77 
78 	hx &= ~(i - 1);
79 	if (rm == FE_UPWARD || (rm == FE_TONEAREST && (frac > (i >> 1) ||
80 		((frac == (i >> 1)) && (hx & i)))))
81 		xx.i = sx | (hx + i);
82 	else
83 		xx.i = sx | hx;
84 	return (xx.f);
85 }
86 
87 #if 0
88 
89 /*
90  * Alternate implementations for SPARC, x86, using fp ops.  These may
91  * be faster depending on how expensive saving and restoring the fp
92  * modes and status flags is.
93  */
94 
95 #include "libm.h"
96 #include "fma.h"
97 
98 #if defined(__sparc)
99 
100 float
101 __nearbyintf(float x) {
102 	union {
103 		unsigned i;
104 		float f;
105 	} xx, yy;
106 	float z;
107 	unsigned hx, sx, fsr, oldfsr;
108 	int rm;
109 
110 	xx.f = x;
111 	sx = xx.i & 0x80000000;
112 	hx = xx.i & ~0x80000000;
113 
114 	/* handle trivial cases */
115 	if (hx >= 0x4b000000)	/* x is nan, inf, or already integral */
116 		return (x + 0.0f);
117 	else if (hx == 0)	/* x is zero */
118 		return (x);
119 
120 	/* save the fsr */
121 	__fenv_getfsr(&oldfsr);
122 
123 	/* handle |x| < 1 */
124 	if (hx < 0x3f800000) {
125 		/* flip the sense of directed roundings if x is negative */
126 		rm = oldfsr >> 30;
127 		if (sx)
128 			rm ^= rm >> 1;
129 		if (rm == FSR_RP || (rm == FSR_RN && hx > 0x3f000000))
130 			xx.i = sx | 0x3f800000;
131 		else
132 			xx.i = sx;
133 		return (xx.f);
134 	}
135 
136 	/* clear the inexact trap */
137 	fsr = oldfsr & ~FSR_NXM;
138 	__fenv_setfsr(&fsr);
139 
140 	/* round x at the integer bit */
141 	yy.i = sx | 0x4b000000;
142 	z = (x + yy.f) - yy.f;
143 
144 	/* restore the old fsr */
145 	__fenv_setfsr(&oldfsr);
146 
147 	return (z);
148 }
149 
150 #elif defined(__x86)
151 
152 /* inline template */
153 extern long double frndint(long double);
154 
155 float
156 __nearbyintf(float x) {
157 	long double z;
158 	unsigned oldcwsw, cwsw;
159 
160 	/* save the control and status words, mask the inexact exception */
161 	__fenv_getcwsw(&oldcwsw);
162 	cwsw = oldcwsw | 0x00200000;
163 	__fenv_setcwsw(&cwsw);
164 
165 	z = frndint((long double) x);
166 
167 	/*
168 	 * restore the control and status words, preserving all but the
169 	 * inexact flag
170 	 */
171 	__fenv_getcwsw(&cwsw);
172 	oldcwsw |= (cwsw & 0x1f);
173 	__fenv_setcwsw(&oldcwsw);
174 
175 	/* note: the value of z is representable in single precision */
176 	return (z);
177 }
178 
179 #else
180 #error Unknown architecture
181 #endif
182 
183 #endif
184