1 /* lrint adapted to be llrint for Newlib, 2009 by Craig Howland. */
2 /* @(#)s_lrint.c 5.1 93/09/24 */
3 /*
4 * ====================================================
5 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
6 *
7 * Developed at SunPro, a Sun Microsystems, Inc. business.
8 * Permission to use, copy, modify, and distribute this
9 * software is freely granted, provided that this notice
10 * is preserved.
11 * ====================================================
12 */
13
14 /*
15 * llrint(x)
16 * Return x rounded to integral value according to the prevailing
17 * rounding mode.
18 * Method:
19 * Using floating addition.
20 * Exception:
21 * Inexact flag raised if x not equal to llrint(x).
22 */
23
24 #include "fdlibm.h"
25 #include <limits.h>
26
27 #ifdef _NEED_FLOAT64
28
29 static const double
30 /* Adding a double, x, to 2^52 will cause the result to be rounded based on
31 the fractional part of x, according to the implementation's current rounding
32 mode. 2^52 is the smallest double that can be represented using all 52 significant
33 digits. */
34 TWO52[2]={
35 _F_64(4.50359962737049600000e+15), /* 0x43300000, 0x00000000 */
36 _F_64(-4.50359962737049600000e+15), /* 0xC3300000, 0x00000000 */
37 };
38
39 long long int
llrint64(__float64 x)40 llrint64(__float64 x)
41 {
42 __int32_t i0,j0,sx;
43 __uint32_t i1;
44 __float64 t;
45 volatile __float64 w;
46 long long int result;
47
48 EXTRACT_WORDS(i0,i1,x);
49
50 /* Extract sign bit. */
51 sx = (i0>>31)&1;
52
53 /* Extract exponent field. */
54 j0 = ((i0 & 0x7ff00000) >> 20) - 1023;
55 /* j0 in [-1023,1024] */
56
57 if(j0 < 20)
58 {
59 /* j0 in [-1023,19] */
60 w = TWO52[sx] + x;
61 t = w - TWO52[sx];
62 GET_HIGH_WORD(i0, t);
63 /* Detect the all-zeros representation of plus and
64 minus zero, which fails the calculation below. */
65 if ((i0 & ~((__int32_t)1 << 31)) == 0)
66 return 0;
67 j0 = ((i0 & 0x7ff00000) >> 20) - 1023;
68 i0 &= 0x000fffff;
69 i0 |= 0x00100000;
70 result = (j0 < 0 ? 0 : i0 >> (20 - j0));
71 }
72 else if (j0 < (int)(8 * sizeof (long long int)) - 1)
73 {
74 /* 64bit return: j0 in [20,62] */
75 if (j0 >= 52)
76 /* 64bit return: j0 in [52,62] */
77 /* 64bit return: left shift amt in [32,42] */
78 result = ((long long int) ((i0 & 0x000fffff) | 0x00100000) << (j0 - 20)) |
79 /* 64bit return: right shift amt in [0,10] */
80 ((long long int) i1 << (j0 - 52));
81 else
82 {
83 /* 64bit return: j0 in [20,51] */
84 w = TWO52[sx] + x;
85 t = w - TWO52[sx];
86 EXTRACT_WORDS (i0, i1, t);
87 j0 = ((i0 & 0x7ff00000) >> 20) - 1023;
88 i0 &= 0x000fffff;
89 i0 |= 0x00100000;
90 /* After round:
91 * 64bit return: j0 in [20,52] */
92 /* 64bit return: left shift amt in [0,32] */
93 /* ***64bit return: right shift amt in [32,0] */
94 result = ((long long int) i0 << (j0 - 20))
95 | SAFE_RIGHT_SHIFT (i1, (52 - j0));
96 }
97 }
98 else
99 {
100 if (sizeof (long long) == 4 && (__float64) LLONG_MIN - _F_64(1.0) < x && x < (__float64) LLONG_MIN) {
101 if (nearbyint(x) == LLONG_MIN)
102 __math_set_inexact64();
103 else
104 __math_set_invalid();
105 return LLONG_MIN;
106 }
107 else if (x != LLONG_MIN)
108 {
109 __math_set_invalid();
110 return sx ? LLONG_MIN : LLONG_MAX;
111 }
112 return (long long) x;
113 }
114
115
116 return sx ? -result : result;
117 }
118
119 _MATH_ALIAS_k_d(llrint)
120
121 #endif /* _NEED_FLOAT64 */
122
