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 2006 Sun Microsystems, Inc. All rights reserved.
24 * Use is subject to license terms.
25 */
26
27 /*
28 * Copyright 2011, Richard Lowe
29 */
30
31 #ifndef _LIBM_INLINES_H
32 #define _LIBM_INLINES_H
33
34 #ifdef __GNUC__
35
36 #ifdef __cplusplus
37 extern "C" {
38 #endif
39
40 #include <sys/types.h>
41 #include <sys/ieeefp.h>
42
43 #define _LO_WORD(x) ((uint32_t *)&x)[0]
44 #define _HI_WORD(x) ((uint32_t *)&x)[1]
45 #define _HIER_WORD(x) ((uint32_t *)&x)[2]
46
47 extern __GNU_INLINE double
__inline_sqrt(double a)48 __inline_sqrt(double a)
49 {
50 double ret;
51
52 __asm__ __volatile__("fsqrt\n\t" : "=t" (ret) : "0" (a) : "cc");
53 return (ret);
54 }
55
56 extern __GNU_INLINE double
__ieee754_sqrt(double a)57 __ieee754_sqrt(double a)
58 {
59 return (__inline_sqrt(a));
60 }
61
62 extern __GNU_INLINE float
__inline_sqrtf(float a)63 __inline_sqrtf(float a)
64 {
65 float ret;
66
67 __asm__ __volatile__("fsqrt\n\t" : "=t" (ret) : "0" (a) : "cc");
68 return (ret);
69 }
70
71 extern __GNU_INLINE double
__inline_rint(double a)72 __inline_rint(double a)
73 {
74 __asm__ __volatile__(
75 "andl $0x7fffffff,%1\n\t"
76 "cmpl $0x43300000,%1\n\t"
77 "jae 1f\n\t"
78 "frndint\n\t"
79 "1: fwait\n\t"
80 : "+t" (a), "+&r" (_HI_WORD(a))
81 :
82 : "cc");
83
84 return (a);
85 }
86
87 /*
88 * 00 - 24 bits
89 * 01 - reserved
90 * 10 - 53 bits
91 * 11 - 64 bits
92 */
93 extern __GNU_INLINE int
__swapRP(int i)94 __swapRP(int i)
95 {
96 int ret;
97 uint16_t cw;
98
99 __asm__ __volatile__("fstcw %0\n\t" : "=m" (cw));
100
101 ret = (cw >> 8) & 0x3;
102 cw = (cw & 0xfcff) | ((i & 0x3) << 8);
103
104 __asm__ __volatile__("fldcw %0\n\t" : : "m" (cw));
105
106 return (ret);
107 }
108
109 /*
110 * 00 - Round to nearest, with even preferred
111 * 01 - Round down
112 * 10 - Round up
113 * 11 - Chop
114 */
115 extern __GNU_INLINE enum fp_direction_type
__swap87RD(enum fp_direction_type i)116 __swap87RD(enum fp_direction_type i)
117 {
118 int ret;
119 uint16_t cw;
120
121 __asm__ __volatile__("fstcw %0\n\t" : "=m" (cw));
122
123 ret = (cw >> 10) & 0x3;
124 cw = (cw & 0xf3ff) | ((i & 0x3) << 10);
125
126 __asm__ __volatile__("fldcw %0\n\t" : : "m" (cw));
127
128 return (ret);
129 }
130
131 extern __GNU_INLINE double
ceil(double d)132 ceil(double d)
133 {
134 /*
135 * Let's set a Rounding Control (RC) bits from x87 FPU Control Word
136 * to fp_positive and save old bits in rd.
137 */
138 short rd = __swap87RD(fp_positive);
139
140 /*
141 * The FRNDINT instruction returns a floating-point value that is the
142 * integral value closest to the source value in the direction of the
143 * rounding mode specified in the RC field of the x87 FPU control word.
144 *
145 * Rounds the source value in the ST(0) register to the nearest
146 * integral value, depending on the current rounding mode
147 * (setting of the RC field of the FPU control word),
148 * and stores the result in ST(0).
149 */
150 __asm__ __volatile__("frndint" : "+t" (d) : : "cc");
151
152 /* restore old RC bits */
153 __swap87RD(rd);
154
155 return (d);
156 }
157
158 extern __GNU_INLINE double
copysign(double d1,double d2)159 copysign(double d1, double d2)
160 {
161 __asm__ __volatile__(
162 "andl $0x7fffffff,%0\n\t" /* %0 <-- hi_32(abs(d)) */
163 "andl $0x80000000,%1\n\t" /* %1[31] <-- sign_bit(d2) */
164 "orl %1,%0\n\t" /* %0 <-- hi_32(copysign(x,y)) */
165 : "+&r" (_HI_WORD(d1)), "+r" (_HI_WORD(d2))
166 :
167 : "cc");
168
169 return (d1);
170 }
171
172 extern __GNU_INLINE double
fabs(double d)173 fabs(double d)
174 {
175 __asm__ __volatile__("fabs\n\t" : "+t" (d) : : "cc");
176 return (d);
177 }
178
179 extern __GNU_INLINE float
fabsf(float d)180 fabsf(float d)
181 {
182 __asm__ __volatile__("fabs\n\t" : "+t" (d) : : "cc");
183 return (d);
184 }
185
186 extern __GNU_INLINE long double
fabsl(long double d)187 fabsl(long double d)
188 {
189 __asm__ __volatile__("fabs\n\t" : "+t" (d) : : "cc");
190 return (d);
191 }
192
193 extern __GNU_INLINE int
finite(double d)194 finite(double d)
195 {
196 int ret = _HI_WORD(d);
197
198 __asm__ __volatile__(
199 "notl %0\n\t"
200 "andl $0x7ff00000,%0\n\t"
201 "negl %0\n\t"
202 "shrl $31,%0\n\t"
203 : "+r" (ret)
204 :
205 : "cc");
206 return (ret);
207 }
208
209 extern __GNU_INLINE double
floor(double d)210 floor(double d)
211 {
212 short rd = __swap87RD(fp_negative);
213
214 __asm__ __volatile__("frndint" : "+t" (d), "+r" (rd) : : "cc");
215 __swap87RD(rd);
216
217 return (d);
218 }
219
220 /*
221 * branchless __isnan
222 * ((0x7ff00000-[((lx|-lx)>>31)&1]|ahx)>>31)&1 = 1 iff x is NaN
223 */
224 extern __GNU_INLINE int
isnan(double d)225 isnan(double d)
226 {
227 int ret;
228
229 __asm__ __volatile__(
230 "movl %1,%%ecx\n\t"
231 "negl %%ecx\n\t" /* ecx <-- -lo_32(x) */
232 "orl %%ecx,%1\n\t"
233 "shrl $31,%1\n\t" /* 1 iff lx != 0 */
234 "andl $0x7fffffff,%2\n\t" /* ecx <-- hi_32(abs(x)) */
235 "orl %2,%1\n\t"
236 "subl $0x7ff00000,%1\n\t"
237 "negl %1\n\t"
238 "shrl $31,%1\n\t"
239 : "=r" (ret)
240 : "0" (_HI_WORD(d)), "r" (_LO_WORD(d))
241 : "ecx");
242
243 return (ret);
244 }
245
246 extern __GNU_INLINE int
isnanf(float f)247 isnanf(float f)
248 {
249 __asm__ __volatile__(
250 "andl $0x7fffffff,%0\n\t"
251 "negl %0\n\t"
252 "addl $0x7f800000,%0\n\t"
253 "shrl $31,%0\n\t"
254 : "+r" (f)
255 :
256 : "cc");
257
258 return (f);
259 }
260
261 extern __GNU_INLINE double
rint(double a)262 rint(double a)
263 {
264 return (__inline_rint(a));
265 }
266
267 extern __GNU_INLINE double
scalbn(double d,int n)268 scalbn(double d, int n)
269 {
270 double dummy;
271
272 __asm__ __volatile__(
273 "fildl %2\n\t" /* Convert N to extended */
274 "fxch\n\t"
275 "fscale\n\t"
276 : "+t" (d), "=u" (dummy)
277 : "m" (n)
278 : "cc");
279
280 return (d);
281 }
282
283 extern __GNU_INLINE int
signbit(double d)284 signbit(double d)
285 {
286 return (_HI_WORD(d) >> 31);
287 }
288
289 extern __GNU_INLINE int
signbitf(float f)290 signbitf(float f)
291 {
292 return ((*(uint32_t *)&f) >> 31);
293 }
294
295 extern __GNU_INLINE double
sqrt(double d)296 sqrt(double d)
297 {
298 return (__inline_sqrt(d));
299 }
300
301 extern __GNU_INLINE float
sqrtf(float f)302 sqrtf(float f)
303 {
304 return (__inline_sqrtf(f));
305 }
306
307 extern __GNU_INLINE long double
sqrtl(long double ld)308 sqrtl(long double ld)
309 {
310 __asm__ __volatile__("fsqrt" : "+t" (ld) : : "cc");
311 return (ld);
312 }
313
314 extern __GNU_INLINE int
isnanl(long double ld)315 isnanl(long double ld)
316 {
317 int ret = _HIER_WORD(ld);
318
319 __asm__ __volatile__(
320 "andl $0x00007fff,%0\n\t"
321 "jz 1f\n\t" /* jump if exp is all 0 */
322 "xorl $0x00007fff,%0\n\t"
323 "jz 2f\n\t" /* jump if exp is all 1 */
324 "testl $0x80000000,%1\n\t"
325 "jz 3f\n\t" /* jump if leading bit is 0 */
326 "xorl %0,%0\n\t"
327 "jmp 1f\n\t"
328 "2:\n\t" /* note that %0 = 0 from before */
329 "cmpl $0x80000000,%1\n\t" /* what is first half of significand? */
330 "jnz 3f\n\t" /* jump if not equal to 0x80000000 */
331 "testl $0xffffffff,%2\n\t" /* is second half of significand 0? */
332 "jnz 3f\n\t" /* jump if not equal to 0 */
333 "jmp 1f\n\t"
334 "3:\n\t"
335 "movl $1,%0\n\t"
336 "1:\n\t"
337 : "+&r" (ret)
338 : "r" (_HI_WORD(ld)), "r" (_LO_WORD(ld))
339 : "cc");
340
341 return (ret);
342 }
343
344 #ifdef __cplusplus
345 }
346 #endif
347
348 #endif /* __GNUC__ */
349
350 #endif /* _LIBM_INLINES_H */
351