1 /* From: @(#)s_floor.c 5.1 93/09/24 */
2 /*
3 * ====================================================
4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5 *
6 * Developed at SunPro, a Sun Microsystems, Inc. business.
7 * Permission to use, copy, modify, and distribute this
8 * software is freely granted, provided that this notice
9 * is preserved.
10 * ====================================================
11 */
12
13 /*
14 * truncl(x)
15 * Return x rounded toward 0 to integral value
16 * Method:
17 * Bit twiddling.
18 * Exception:
19 * Inexact flag raised if x not equal to truncl(x).
20 */
21
22
23
24
25 #ifdef LDBL_IMPLICIT_NBIT
26 #define MANH_SIZE (EXT_FRACHBITS + 1)
27 #else
28 #define MANH_SIZE EXT_FRACHBITS
29 #endif
30
31 static const long double huge = 1.0e300L;
32 static const float zero[] = { 0.0, -0.0 };
33
34 long double
truncl(long double x)35 truncl(long double x)
36 {
37 int e, es;
38 uint32_t ix0, ix1;
39
40 GET_LDOUBLE_WORDS(es,ix0,ix1,x);
41 e = (es&0x7fff) - LDBL_MAX_EXP + 1;
42
43 if (e < MANH_SIZE - 1) {
44 if (e < 0) { /* raise inexact if x != 0 */
45 if (huge + x > 0.0L)
46 return ((long double)zero[(es&0x8000)!=0]);
47 } else {
48 uint64_t m = ((1llu << MANH_SIZE) - 1) >> (e + 1);
49 if (((ix0 & m) | ix1) == 0)
50 return (x); /* x is integral */
51 if (huge + x > 0.0L) { /* raise inexact flag */
52 ix0 &= ~m;
53 ix1 = 0;
54 }
55 }
56 } else if (e < LDBL_MANT_DIG - 1) {
57 uint64_t m = (uint64_t)-1 >> (64 - LDBL_MANT_DIG + e + 1);
58 if ((ix1 & m) == 0)
59 return (x); /* x is integral */
60 if (huge + x > 0.0L) /* raise inexact flag */
61 ix1 &= ~m;
62 } else if (e == 0x7fff - LDBL_MAX_EXP + 1) {
63 return x + x;
64 }
65 SET_LDOUBLE_WORDS(x,es,ix0,ix1);
66 return (x);
67 }
68
