math: fix tgamma to raise underflow for large negative values
[oweals/musl.git] / src / math / rintl.c
1 /* origin: FreeBSD /usr/src/lib/msun/src/s_rintl.c */
2 /*-
3  * Copyright (c) 2008 David Schultz <das@FreeBSD.ORG>
4  * All rights reserved.
5  *
6  * Redistribution and use in source and binary forms, with or without
7  * modification, are permitted provided that the following conditions
8  * are met:
9  * 1. Redistributions of source code must retain the above copyright
10  *    notice, this list of conditions and the following disclaimer.
11  * 2. Redistributions in binary form must reproduce the above copyright
12  *    notice, this list of conditions and the following disclaimer in the
13  *    documentation and/or other materials provided with the distribution.
14  *
15  * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
16  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
17  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
18  * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
19  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
20  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
21  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
22  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
23  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
24  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
25  * SUCH DAMAGE.
26  */
27
28 #include "libm.h"
29
30 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
31 long double rintl(long double x)
32 {
33         return rint(x);
34 }
35 #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
36
37 #define BIAS    (LDBL_MAX_EXP - 1)
38
39 static const float
40 shift[2] = {
41 #if LDBL_MANT_DIG == 64
42         0x1.0p63, -0x1.0p63
43 #elif LDBL_MANT_DIG == 113
44         0x1.0p112, -0x1.0p112
45 #endif
46 };
47 static const float zero[2] = { 0.0, -0.0 };
48
49 long double rintl(long double x)
50 {
51         union IEEEl2bits u;
52         uint32_t expsign;
53         int ex, sign;
54
55         u.e = x;
56         expsign = u.xbits.expsign;
57         ex = expsign & 0x7fff;
58
59         if (ex >= BIAS + LDBL_MANT_DIG - 1) {
60                 if (ex == BIAS + LDBL_MAX_EXP)
61                         return x + x; /* Inf, NaN, or unsupported format */
62                 return x;             /* finite and already an integer */
63         }
64         sign = expsign >> 15;
65
66         /*
67          * The following code assumes that intermediate results are
68          * evaluated in long double precision. If they are evaluated in
69          * greater precision, double rounding may occur, and if they are
70          * evaluated in less precision (as on i386), results will be
71          * wildly incorrect.
72          */
73         x += shift[sign];
74         x -= shift[sign];
75
76         /*
77          * If the result is +-0, then it must have the same sign as x, but
78          * the above calculation doesn't always give this.  Fix up the sign.
79          */
80         if (ex < BIAS && x == 0.0)
81                 return zero[sign];
82
83         return x;
84 }
85 #endif