1 /* Round to integer type. ldbl-128ibm version.
2 Copyright (C) 2016-2023 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 <errno.h>
20 #include <fenv.h>
21 #include <math.h>
22 #include <math_private.h>
23 #include <stdbool.h>
24 #include <stdint.h>
25
26 #define BIAS 0x3ff
27 #define MANT_DIG 53
28
29 #if UNSIGNED
30 # define RET_TYPE uintmax_t
31 #else
32 # define RET_TYPE intmax_t
33 #endif
34
35 #include <fromfp.h>
36
37 RET_TYPE
38 FUNC (long double x, int round, unsigned int width)
39 {
40 double hi, lo;
41 if (width > INTMAX_WIDTH)
42 width = INTMAX_WIDTH;
43 uint64_t hx, lx;
44 ldbl_unpack (x, &hi, &lo);
45 EXTRACT_WORDS64 (hx, hi);
46 EXTRACT_WORDS64 (lx, lo);
47 bool negative = (hx & 0x8000000000000000ULL) != 0;
48 bool lo_negative = (lx & 0x8000000000000000ULL) != 0;
49 if (width == 0)
50 return fromfp_domain_error (negative, width);
51 hx &= 0x7fffffffffffffffULL;
52 lx &= 0x7fffffffffffffffULL;
53 if ((hx | lx) == 0)
54 return 0;
55 int hi_exponent = hx >> (MANT_DIG - 1);
56 hi_exponent -= BIAS;
57 int exponent = hi_exponent;
58 hx &= ((1ULL << (MANT_DIG - 1)) - 1);
59 if (hx == 0 && lx != 0 && lo_negative != negative)
60 exponent--;
61 int max_exponent = fromfp_max_exponent (negative, width);
62 if (exponent > max_exponent)
63 return fromfp_domain_error (negative, width);
64 int lo_exponent = lx >> (MANT_DIG - 1);
65 lo_exponent -= BIAS;
66
67 /* Convert the high part to integer. */
68 hx |= 1ULL << (MANT_DIG - 1);
69 uintmax_t uret;
70 bool half_bit, more_bits;
71 if (hi_exponent >= MANT_DIG - 1)
72 {
73 uret = hx;
74 uret <<= hi_exponent - (MANT_DIG - 1);
75 half_bit = false;
76 more_bits = false;
77 }
78 else if (hi_exponent >= -1)
79 {
80 uint64_t h = 1ULL << (MANT_DIG - 2 - hi_exponent);
81 half_bit = (hx & h) != 0;
82 more_bits = (hx & (h - 1)) != 0;
83 uret = hx >> (MANT_DIG - 1 - hi_exponent);
84 }
85 else
86 {
87 uret = 0;
88 half_bit = false;
89 more_bits = true;
90 }
91
92 /* Likewise, the low part. */
93 if (lx != 0)
94 {
95 uintmax_t lo_uret;
96 bool lo_half_bit, lo_more_bits;
97 lx &= ((1ULL << (MANT_DIG - 1)) - 1);
98 lx |= 1ULL << (MANT_DIG - 1);
99 /* The high part exponent is at most 64, so the low part
100 exponent is at most 11. */
101 if (lo_exponent >= -1)
102 {
103 uint64_t h = 1ULL << (MANT_DIG - 2 - lo_exponent);
104 lo_half_bit = (lx & h) != 0;
105 lo_more_bits = (lx & (h - 1)) != 0;
106 lo_uret = lx >> (MANT_DIG - 1 - lo_exponent);
107 }
108 else
109 {
110 lo_uret = 0;
111 lo_half_bit = false;
112 lo_more_bits = true;
113 }
114 if (lo_negative == negative)
115 {
116 uret += lo_uret;
117 half_bit |= lo_half_bit;
118 more_bits |= lo_more_bits;
119 }
120 else
121 {
122 uret -= lo_uret;
123 if (lo_half_bit)
124 {
125 uret--;
126 half_bit = true;
127 }
128 if (lo_more_bits && !more_bits)
129 {
130 if (half_bit)
131 {
132 half_bit = false;
133 more_bits = true;
134 }
135 else
136 {
137 uret--;
138 half_bit = true;
139 more_bits = true;
140 }
141 }
142 }
143 }
144
145 return fromfp_round_and_return (negative, uret, half_bit, more_bits, round,
146 exponent, max_exponent, width);
147 }