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
26 #ifdef _NEED_FLOAT64
27
28 static const double
29 /* Adding a double, x, to 2^52 will cause the result to be rounded based on
30 the fractional part of x, according to the implementation's current rounding
31 mode. 2^52 is the smallest double that can be represented using all 52 significant
32 digits. */
33 TWO52[2]={
34 _F_64(4.50359962737049600000e+15), /* 0x43300000, 0x00000000 */
35 _F_64(-4.50359962737049600000e+15), /* 0xC3300000, 0x00000000 */
36 };
37
38 long long int
llrint64(__float64 x)39 llrint64(__float64 x)
40 {
41 __int32_t i0,j0,sx;
42 __uint32_t i1;
43 __float64 t;
44 volatile __float64 w;
45 long long int result;
46
47 EXTRACT_WORDS(i0,i1,x);
48
49 /* Extract sign bit. */
50 sx = (i0>>31)&1;
51
52 /* Extract exponent field. */
53 j0 = ((i0 & 0x7ff00000) >> 20) - 1023;
54 /* j0 in [-1023,1024] */
55
56 if(j0 < 20)
57 {
58 /* j0 in [-1023,19] */
59 if(j0 < -1)
60 return 0;
61 else
62 {
63 /* j0 in [0,19] */
64 /* shift amt in [0,19] */
65 w = TWO52[sx] + x;
66 t = w - TWO52[sx];
67 GET_HIGH_WORD(i0, t);
68 /* Detect the all-zeros representation of plus and
69 minus zero, which fails the calculation below. */
70 if ((i0 & ~((__int32_t)1 << 31)) == 0)
71 return 0;
72 /* After round: j0 in [0,20] */
73 j0 = ((i0 & 0x7ff00000) >> 20) - 1023;
74 i0 &= 0x000fffff;
75 i0 |= 0x00100000;
76 /* shift amt in [20,0] */
77 result = i0 >> (20 - j0);
78 }
79 }
80 else if (j0 < (int)(8 * sizeof (long long int)) - 1)
81 {
82 /* 64bit return: j0 in [20,62] */
83 if (j0 >= 52)
84 /* 64bit return: j0 in [52,62] */
85 /* 64bit return: left shift amt in [32,42] */
86 result = ((long long int) ((i0 & 0x000fffff) | 0x00100000) << (j0 - 20)) |
87 /* 64bit return: right shift amt in [0,10] */
88 ((long long int) i1 << (j0 - 52));
89 else
90 {
91 /* 64bit return: j0 in [20,51] */
92 w = TWO52[sx] + x;
93 t = w - TWO52[sx];
94 EXTRACT_WORDS (i0, i1, t);
95 j0 = ((i0 & 0x7ff00000) >> 20) - 1023;
96 i0 &= 0x000fffff;
97 i0 |= 0x00100000;
98 /* After round:
99 * 64bit return: j0 in [20,52] */
100 /* 64bit return: left shift amt in [0,32] */
101 /* ***64bit return: right shift amt in [32,0] */
102 result = ((long long int) i0 << (j0 - 20))
103 | SAFE_RIGHT_SHIFT (i1, (52 - j0));
104 }
105 }
106 else
107 {
108 return (long long int) x;
109 }
110
111 return sx ? -result : result;
112 }
113
114 _MATH_ALIAS_k_d(llrint)
115
116 #endif /* _NEED_FLOAT64 */
117
