1 | /* Round double value to long long int. |
2 | Copyright (C) 1997-2024 Free Software Foundation, Inc. |
3 | This file is part of the GNU C Library. |
4 | |
5 | The GNU C Library is free software; you can redistribute it and/or |
6 | modify it under the terms of the GNU Lesser General Public |
7 | License as published by the Free Software Foundation; either |
8 | version 2.1 of the License, or (at your option) any later version. |
9 | |
10 | The GNU C Library is distributed in the hope that it will be useful, |
11 | but WITHOUT ANY WARRANTY; without even the implied warranty of |
12 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
13 | Lesser General Public License for more details. |
14 | |
15 | You should have received a copy of the GNU Lesser General Public |
16 | License along with the GNU C Library; if not, see |
17 | <https://www.gnu.org/licenses/>. */ |
18 | |
19 | #include <limits.h> |
20 | #include <math.h> |
21 | #include <math_private.h> |
22 | #include <stdint.h> |
23 | #include <libm-alias-double.h> |
24 | #include <math-barriers.h> |
25 | |
26 | /* Round to the nearest integer, with values exactly on a 0.5 boundary |
27 | rounded away from zero, regardless of the current rounding mode. |
28 | If (long long)x, when x is out of range of a long long, clips at |
29 | LLONG_MAX or LLONG_MIN, then this implementation also clips. */ |
30 | |
31 | long long int |
32 | __llround (double x) |
33 | { |
34 | #ifdef _ARCH_PWR5X |
35 | x = round (x); |
36 | /* The barrier prevents compiler from optimizing it to llround when |
37 | compiled with -fno-math-errno */ |
38 | math_opt_barrier (x); |
39 | return x; |
40 | #else |
41 | long long xr; |
42 | if (HAVE_PPC_FCTIDZ) |
43 | { |
44 | /* IEEE 1003.1 lround function. IEEE specifies "round to the nearest |
45 | integer value, rounding halfway cases away from zero, regardless of |
46 | the current rounding mode." However PowerPC Architecture defines |
47 | "round to Nearest" as "Choose the best approximation. In case of a |
48 | tie, choose the one that is even (least significant bit o).". |
49 | So we can't use the PowerPC "round to Nearest" mode. Instead we set |
50 | "round toward Zero" mode and round by adding +-0.5 before rounding |
51 | to the integer value. |
52 | |
53 | It is necessary to detect when x is (+-)0x1.fffffffffffffp-2 |
54 | because adding +-0.5 in this case will cause an erroneous shift, |
55 | carry and round. We simply return 0 if 0.5 > x > -0.5. Likewise |
56 | if x is and odd number between +-(2^52 and 2^53-1) a shift and |
57 | carry will erroneously round if biased with +-0.5. Therefore if x |
58 | is greater/less than +-2^52 we don't need to bias the number with |
59 | +-0.5. */ |
60 | double ax = fabs (x: x); |
61 | |
62 | if (ax < 0.5) |
63 | return 0; |
64 | |
65 | if (ax < 0x1p+52) |
66 | { |
67 | /* Test whether an integer to avoid spurious "inexact". */ |
68 | double t = ax + 0x1p+52; |
69 | t = t - 0x1p+52; |
70 | if (ax != t) |
71 | { |
72 | ax = ax + 0.5; |
73 | if (x < 0.0) |
74 | ax = -fabs (x: ax); |
75 | x = ax; |
76 | } |
77 | } |
78 | |
79 | return x; |
80 | } |
81 | else |
82 | { |
83 | /* Avoid incorrect exceptions from libgcc conversions (as of GCC |
84 | 5): <https://gcc.gnu.org/bugzilla/show_bug.cgi?id=59412>. */ |
85 | if (fabs (x: x) < 0x1p31) |
86 | xr = (long long int) (long int) x; |
87 | else |
88 | { |
89 | uint64_t i0; |
90 | EXTRACT_WORDS64 (i0, x); |
91 | int exponent = ((i0 >> 52) & 0x7ff) - 0x3ff; |
92 | if (exponent < 63) |
93 | { |
94 | unsigned long long int mant |
95 | = (i0 & ((1ULL << 52) - 1)) | (1ULL << 52); |
96 | if (exponent < 52) |
97 | /* llround is not required to raise "inexact". */ |
98 | mant >>= 52 - exponent; |
99 | else |
100 | mant <<= exponent - 52; |
101 | xr = (long long int) ((i0 & (1ULL << 63)) != 0 ? -mant : mant); |
102 | } |
103 | else if (x == (double) LLONG_MIN) |
104 | xr = LLONG_MIN; |
105 | else |
106 | xr = (long long int) (long int) x << 32; |
107 | } |
108 | } |
109 | /* Avoid spurious "inexact" converting LLONG_MAX to double, and from |
110 | subtraction when the result is out of range, by returning early |
111 | for arguments large enough that no rounding is needed. */ |
112 | if (!(fabs (x: x) < 0x1p52)) |
113 | return xr; |
114 | double xrf = (double) xr; |
115 | |
116 | if (x >= 0.0) |
117 | { |
118 | if (x - xrf >= 0.5) |
119 | xr += (long long) ((unsigned long long) xr + 1) > 0; |
120 | } |
121 | else |
122 | { |
123 | if (xrf - x >= 0.5) |
124 | xr -= (long long) ((unsigned long long) xr - 1) < 0; |
125 | } |
126 | return xr; |
127 | #endif |
128 | } |
129 | #ifndef __llround |
130 | libm_alias_double (__llround, llround) |
131 | #endif |
132 | |