math: fix tgamma to raise underflow for large negative values
[oweals/musl.git] / src / math / trunc.c
1 /* origin: FreeBSD /usr/src/lib/msun/src/s_trunc.c */
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  * trunc(x)
14  * Return x rounded toward 0 to integral value
15  * Method:
16  *      Bit twiddling.
17  * Exception:
18  *      Inexact flag raised if x not equal to trunc(x).
19  */
20
21 #include "libm.h"
22
23 static const double huge = 1.0e300;
24
25 double trunc(double x)
26 {
27         int32_t i0,i1,j0;
28         uint32_t i;
29
30         EXTRACT_WORDS(i0, i1, x);
31         j0 = ((i0>>20)&0x7ff) - 0x3ff;
32         if (j0 < 20) {
33                 if (j0 < 0) { /* |x|<1, return 0*sign(x) */
34                         /* raise inexact if x != 0 */
35                         if (huge+x > 0.0) {
36                                 i0 &= 0x80000000U;
37                                 i1 = 0;
38                         }
39                 } else {
40                         i = 0x000fffff>>j0;
41                         if (((i0&i)|i1) == 0)
42                                 return x; /* x is integral */
43                         /* raise inexact */
44                         if (huge+x > 0.0) {
45                                 i0 &= ~i;
46                                 i1 = 0;
47                         }
48                 }
49         } else if (j0 > 51) {
50                 if (j0 == 0x400)
51                         return x + x;  /* inf or NaN */
52                 return x;              /* x is integral */
53         } else {
54                 i = (uint32_t)0xffffffff>>(j0-20);
55                 if ((i1&i) == 0)
56                         return x;      /* x is integral */
57                 /* raise inexact */
58                 if (huge+x > 0.0)
59                         i1 &= ~i;
60         }
61         INSERT_WORDS(x, i0, i1);
62         return x;
63 }