1 /*
2  * SPDX-License-Identifier: BSD-3-Clause
3  *
4  * Copyright © 2023 Keith Packard
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  *
10  * 1. Redistributions of source code must retain the above copyright
11  *    notice, this list of conditions and the following disclaimer.
12  *
13  * 2. Redistributions in binary form must reproduce the above
14  *    copyright notice, this list of conditions and the following
15  *    disclaimer in the documentation and/or other materials provided
16  *    with the distribution.
17  *
18  * 3. Neither the name of the copyright holder nor the names of its
19  *    contributors may be used to endorse or promote products derived
20  *    from this software without specific prior written permission.
21  *
22  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
23  * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
24  * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
25  * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
26  * COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
27  * INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
28  * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
29  * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
31  * STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
32  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
33  * OF THE POSSIBILITY OF SUCH DAMAGE.
34  */
35 
36 /*
37  * Convert an 80-bit or 128-bit float to a string in hex ('a') format
38  *
39  * This chunk of code gets inserted into the vfprintf function in the
40  * long double case when long double is larger than 64 bits.
41  *
42  * This code uses the 'u128' type to hold the floating point value as
43  * an integer value.
44  *
45  * This code only works with long double type.
46  */
47 
48 #include "stdio_private.h"
49 
50 #if __SIZEOF_LONG_DOUBLE__ > 8
51 
52 #define _NEED_IO_FLOAT_LARGE
53 
54 #include "dtoa.h"
55 
56 #if __LDBL_MANT_DIG__ == 64
57 # define LEXP_BIAS      (__LDBL_MAX_EXP__ + 2)
58 # define LEXP_INF       (__LDBL_MAX_EXP__ - 3)
59 # define LSIG_BITS      (__LDBL_MANT_DIG__)
60 # ifdef __m68k__
61 #  define LDENORM_EXP_BIAS 0
62 # else
63 #  define LDENORM_EXP_BIAS 1
64 #  define LSIG_MSB_INF    _u128_lshift(to_u128(1), __LDBL_MANT_DIG__-1)
65 # endif
66 # if __BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__
67 #  define LEXP_SHIFT       __LDBL_MANT_DIG__
68 #  define LSIGN_SHIFT        79
69 # endif
70 # if __BYTE_ORDER__ == __ORDER_BIG_ENDIAN__
71 #  define LEXP_SHIFT       (__LDBL_MANT_DIG__ + 16)
72 #  define LSIGN_SHIFT        (79 + 16)
73 # endif
74 #else
75 # define LDENORM_EXP_BIAS 1
76 # define LSIGN_SHIFT        127
77 # define LEXP_BIAS       (__LDBL_MAX_EXP__ - 1)
78 # define LEXP_INF        (__LDBL_MAX_EXP__)
79 # define LSIG_MSB        _u128_lshift(to_u128(1), __LDBL_MANT_DIG__-1)
80 # define LSIG_BITS       (__LDBL_MANT_DIG__ - 1)
81 # define LEXP_SHIFT      (__LDBL_MANT_DIG__ - 1)
82 #endif
83 
84 #define LEXP_MASK        ((__LDBL_MAX_EXP__ - 1) + __LDBL_MAX_EXP__)
85 #define LSIG_SHIFT       0
86 #define LSIG_MASK        _u128_minus_64(_u128_lshift(to_u128(1), LSIG_BITS), 1)
87 
88 #define TOCASE(c)       ((c) - case_convert)
89 
90 #define LDTOX_NDIGS     (__LDBL_MANT_DIG__ + 3) / 4
91 
92 int
__ldtox_engine(long double x,struct dtoa * dtoa,int prec,unsigned char case_convert)93 __ldtox_engine(long double x, struct dtoa *dtoa, int prec, unsigned char case_convert)
94 {
95     _u128 fi, s;
96     int exp;
97 
98     dtoa->flags = 0;
99     fi = asuintld(x);
100     if (_u128_and_64(_u128_rshift(fi, LSIGN_SHIFT), 1))
101         dtoa->flags = DTOA_MINUS;
102 
103     exp = _u128_and_64(_u128_rshift(fi, LEXP_SHIFT), LEXP_MASK);
104     s = fi = _u128_lshift(_u128_and(fi, LSIG_MASK), LSIG_SHIFT);
105     if (!_u128_is_zero(s) || exp != 0) {
106         if (exp == 0)
107             exp = LDENORM_EXP_BIAS;
108         else
109         {
110 #ifdef LSIG_MSB
111             s = _u128_or(s, LSIG_MSB);
112 #endif
113         }
114         exp -= LEXP_BIAS;
115     }
116 
117     if (prec < 0)
118         prec = 0;
119     else if (prec >= (LDTOX_NDIGS - 1))
120         prec = LDTOX_NDIGS - 1;
121     else {
122         int     bits = ((LDTOX_NDIGS - 1) - prec) << 2;
123         _u128   half = _u128_lshift(to_u128(1), bits - 1);
124         _u128   mask = _u128_not(_u128_minus_64(_u128_lshift(half, 1), 1));
125 
126         /* round even */
127         if (_u128_gt(_u128_and(s, _u128_not(mask)), half) || _u128_and_64(_u128_rshift(s, bits), 1) != 0) {
128             s = _u128_plus(s, half);
129 
130 #if !defined(LSIG_MSB) && (__LDBL_MANT_DIG__ & 3) == 0
131             /*
132              * Formats without an implicit '1' for the MSB and which
133              * fill all four bits of the top digit might actually
134              * overflow when rounding. Check for that and shift right.
135              * We can do this without loss of accuracy because we just
136              * carried all the way from 'half' to the top digit of the
137              * value
138              */
139             if (_u128_ge(s, _u128_lshift(to_u128(1), __LDBL_MANT_DIG__))) {
140                 s = _u128_rshift(s, 4);
141                 exp += 4;
142             }
143 #endif
144         }
145 
146         s = _u128_and(s, mask);
147     }
148 
149     if (exp == LEXP_INF) {
150 #ifdef LSIG_MSB_INF
151         if (!_u128_eq(fi, LSIG_MSB_INF))
152 #else
153             if (!_u128_is_zero(fi))
154 #endif
155                 dtoa->flags |= DTOA_NAN;
156             else
157                 dtoa->flags |= DTOA_INF;
158     } else {
159         int8_t d;
160         for (d = LDTOX_NDIGS - 1; d >= 0; d--) {
161             int dig = _u128_and_64(s, 0xf);
162             s = _u128_rshift(s, 4);
163             if (dig == 0 && d > prec)
164                 continue;
165             if (dig <= 9)
166                 dig += '0';
167             else
168                 dig += TOCASE('a' - 10);
169             dtoa->digits[d] = dig;
170             if (prec < d)
171                 prec = d;
172         }
173     }
174     dtoa->exp = exp;
175     return prec;
176 }
177 
178 #endif
179