src/libm/s_floor.c
 author Sam Lantinga Sat, 19 Sep 2009 13:29:40 +0000 changeset 3280 00cace2d9080 parent 3162 dc1eb82ffdaa child 3337 9ac6f0782dd6 permissions -rw-r--r--
Merged a cleaned up version of Jiang's code changes from Google Summer of Code 2009
1 /* @(#)s_floor.c 5.1 93/09/24 */
2 /*
3  * ====================================================
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  */
13 #if defined(LIBM_SCCS) && !defined(lint)
14 static const char rcsid[] =
15     "\$NetBSD: s_floor.c,v 1.8 1995/05/10 20:47:20 jtc Exp \$";
16 #endif
18 /*
19  * floor(x)
20  * Return x rounded toward -inf to integral value
21  * Method:
22  *	Bit twiddling.
23  * Exception:
24  *	Inexact flag raised if x not equal to floor(x).
25  */
27 #include "math.h"
28 #include "math_private.h"
30 #ifdef __STDC__
31 static const double huge = 1.0e300;
32 #else
33 static double huge = 1.0e300;
34 #endif
36 libm_hidden_proto(floor)
37 #ifdef __STDC__
38      double floor(double x)
39 #else
40      double floor(x)
41      double x;
42 #endif
43 {
44     int32_t i0, i1, j0;
45     u_int32_t i, j;
46     EXTRACT_WORDS(i0, i1, x);
47     j0 = ((i0 >> 20) & 0x7ff) - 0x3ff;
48     if (j0 < 20) {
49         if (j0 < 0) {           /* raise inexact if x != 0 */
50             if (huge + x > 0.0) {       /* return 0*sign(x) if |x|<1 */
51                 if (i0 >= 0) {
52                     i0 = i1 = 0;
53                 } else if (((i0 & 0x7fffffff) | i1) != 0) {
54                     i0 = 0xbff00000;
55                     i1 = 0;
56                 }
57             }
58         } else {
59             i = (0x000fffff) >> j0;
60             if (((i0 & i) | i1) == 0)
61                 return x;       /* x is integral */
62             if (huge + x > 0.0) {       /* raise inexact flag */
63                 if (i0 < 0)
64                     i0 += (0x00100000) >> j0;
65                 i0 &= (~i);
66                 i1 = 0;
67             }
68         }
69     } else if (j0 > 51) {
70         if (j0 == 0x400)
71             return x + x;       /* inf or NaN */
72         else
73             return x;           /* x is integral */
74     } else {
75         i = ((u_int32_t) (0xffffffff)) >> (j0 - 20);
76         if ((i1 & i) == 0)
77             return x;           /* x is integral */
78         if (huge + x > 0.0) {   /* raise inexact flag */
79             if (i0 < 0) {
80                 if (j0 == 20)
81                     i0 += 1;
82                 else {
83                     j = i1 + (1 << (52 - j0));
84                     if (j < (u_int32_t) i1)
85                         i0 += 1;        /* got a carry */
86                     i1 = j;
87                 }
88             }
89             i1 &= (~i);
90         }
91     }
92     INSERT_WORDS(x, i0, i1);
93     return x;
94 }
96 libm_hidden_def(floor)