Branch data Line data Source code
1 : : // Copyright 2019-2026 David Robillard <d@drobilla.net>
2 : : // SPDX-License-Identifier: ISC
3 : :
4 : : #include "digits.h"
5 : :
6 : : #include "attributes.h"
7 : : #include "bigint.h"
8 : : #include "ieee_float.h"
9 : : #include "int_math.h"
10 : : #include "soft_float.h"
11 : : #include "warnings.h"
12 : :
13 : : #include <assert.h>
14 : : #include <math.h>
15 : : #include <stdbool.h>
16 : : #include <stddef.h>
17 : : #include <stdint.h>
18 : :
19 : : /*
20 : : This is more or less just an implementation of the classic rational number
21 : : based floating point print routine ("Dragon4"). See "How to Print
22 : : Floating-Point Numbers Accurately" by Guy L. Steele Jr. and Jon L White for
23 : : the canonical source. The basic idea is to find a big rational between 1 and
24 : : 10 where value = (numer / denom) * 10^e, then continuously divide it to
25 : : generate decimal digits.
26 : :
27 : : Unfortunately, this algorithm requires pretty massive bigints to work
28 : : correctly for all doubles, and isn't particularly fast. Something like
29 : : Grisu3 could be added to improve performance, but that has the annoying
30 : : property of needing a more precise fallback in some cases, meaning it would
31 : : only add more code, not replace any. Since this is already a pretty
32 : : ridiculous amount of code, I'll hold off on this until it becomes a problem,
33 : : or somebody comes up with a better algorithm.
34 : : */
35 : :
36 : : /// Return true if the number is within the lower boundary
37 : : static bool
38 : 406901 : within_lower(const ExessBigint* const numer,
39 : : const ExessBigint* const d_lower,
40 : : const bool is_even)
41 : : {
42 : 271887 : return is_even ? exess_bigint_compare(numer, d_lower) <= 0
43 [ + + ]: 406901 : : exess_bigint_compare(numer, d_lower) < 0;
44 : : }
45 : :
46 : : /// Return true if the number is within the upper boundary
47 : : static bool
48 : 437568 : within_upper(const ExessBigint* const numer,
49 : : const ExessBigint* const denom,
50 : : const ExessBigint* const d_upper,
51 : : const bool is_even)
52 : : {
53 : 294306 : return is_even ? exess_bigint_plus_compare(numer, d_upper, denom) >= 0
54 [ + + ]: 437568 : : exess_bigint_plus_compare(numer, d_upper, denom) > 0;
55 : : }
56 : :
57 : : /**
58 : : Find values so that 0.1 <= numer/denom < 1 or 1 <= numer/denom < 10.
59 : :
60 : : @param significand Double significand.
61 : : @param exponent Double exponent (base 2).
62 : : @param decimal_power Decimal exponent (log10 of the double).
63 : : @param[out] numer Numerator of rational number.
64 : : @param[out] denom Denominator of rational number.
65 : : @param[out] delta Distance to the lower and upper boundaries.
66 : : */
67 : : static EXESS_I_NONBLOCKING void
68 : 30667 : calculate_initial_values(const uint64_t significand,
69 : : const int exponent,
70 : : const int decimal_power,
71 : : const bool lower_is_closer,
72 : : ExessBigint* const numer,
73 : : ExessBigint* const denom,
74 : : ExessBigint* const delta)
75 : : {
76 : : /* Use a common denominator of 2^1 so that boundary distance is an integer.
77 : : If the lower boundary is closer, we need to scale everything but the
78 : : lower boundary to compensate, so add another factor of two here (this is
79 : : faster than shifting them again later as in the paper). */
80 : 30667 : const unsigned lg_denom = 1U + lower_is_closer;
81 : :
82 [ + + ]: 30667 : if (exponent >= 0) {
83 : : // delta = 2^e
84 : 13012 : exess_bigint_set_u32(delta, 1);
85 : 13012 : exess_bigint_shift_left(delta, (unsigned)exponent);
86 : :
87 : : // numer = f * 2^e
88 : 13012 : exess_bigint_set_u64(numer, significand);
89 : 13012 : exess_bigint_shift_left(numer, (unsigned)exponent + lg_denom);
90 : :
91 : : // denom = 10^d
92 : 13012 : exess_bigint_set_pow10(denom, (unsigned)decimal_power);
93 : 13012 : exess_bigint_shift_left(denom, lg_denom);
94 [ + + ]: 17655 : } else if (decimal_power >= 0) {
95 : : // delta = 2^e, which is just 1 here since 2^-e is in the denominator
96 : 3750 : exess_bigint_set_u32(delta, 1);
97 : :
98 : : // numer = f
99 : 3750 : exess_bigint_set_u64(numer, significand);
100 : 3750 : exess_bigint_shift_left(numer, lg_denom);
101 : :
102 : : // denom = 10^d * 2^-e
103 : 3750 : exess_bigint_set_pow10(denom, (unsigned)decimal_power);
104 : 3750 : exess_bigint_shift_left(denom, (unsigned)-exponent + lg_denom);
105 : : } else {
106 : : // delta = 10^d
107 : 13905 : exess_bigint_set_pow10(delta, (unsigned)-decimal_power);
108 : :
109 : : // numer = f * 10^-d
110 : 13905 : exess_bigint_set(numer, delta);
111 : 13905 : exess_bigint_multiply_u64(numer, significand);
112 : 13905 : exess_bigint_shift_left(numer, lg_denom);
113 : :
114 : : // denom = 2^-exponent
115 : 13905 : exess_bigint_set_u32(denom, 1);
116 : 13905 : exess_bigint_shift_left(denom, (unsigned)-exponent + lg_denom);
117 : : }
118 : 30667 : }
119 : :
120 : : #ifndef NDEBUG
121 : : static EXESS_I_NONBLOCKING bool
122 : 30667 : check_initial_values(const ExessBigint* const numer,
123 : : const ExessBigint* const denom,
124 : : const ExessBigint* const d_upper)
125 : : {
126 : 30667 : ExessBigint upper = *numer;
127 : 30667 : exess_bigint_add(&upper, d_upper);
128 [ - + ]: 30667 : assert(exess_bigint_compare(&upper, denom) >= 0);
129 : :
130 : 30667 : const uint32_t div = exess_bigint_divmod(&upper, denom);
131 [ + - - + ]: 30667 : assert(div >= 1 && div < 10);
132 : 30667 : return true;
133 : : }
134 : : #endif
135 : :
136 : : static unsigned
137 : 30667 : emit_digits(ExessBigint* const numer,
138 : : const ExessBigint* const denom,
139 : : ExessBigint* const d_lower,
140 : : ExessBigint* const d_upper,
141 : : const bool is_even,
142 : : const size_t digits_size,
143 : : char* const digits)
144 : : {
145 : 30667 : unsigned length = 0;
146 [ + + ]: 418136 : for (size_t i = 0; i < digits_size; ++i) {
147 : : // Emit the next digit
148 : 406901 : const uint32_t digit = exess_bigint_divmod(numer, denom);
149 [ - + ]: 406901 : assert(digit <= 9);
150 : 406901 : digits[length++] = (char)('0' + digit);
151 : :
152 : : // Check for termination
153 : 406901 : const bool within_low = within_lower(numer, d_lower, is_even);
154 : 406901 : const bool within_high = within_upper(numer, denom, d_upper, is_even);
155 [ + + + + ]: 406901 : if (!within_low && !within_high) {
156 : 387469 : exess_bigint_multiply_u32(numer, 10);
157 : 387469 : exess_bigint_multiply_u32(d_lower, 10);
158 [ + + ]: 387469 : if (d_lower != d_upper) {
159 : 31715 : exess_bigint_multiply_u32(d_upper, 10);
160 : : }
161 : : } else {
162 [ + + + + : 19432 : if (!within_low || (within_high && exess_bigint_plus_compare(
+ + ]
163 : : numer, numer, denom) >= 0)) {
164 : : // In high only, or halfway and the next digit is > 5, round up
165 [ - + ]: 9214 : assert(digits[length - 1] != '9');
166 : 9214 : digits[length - 1]++;
167 : : }
168 : :
169 : 19432 : break;
170 : : }
171 : : }
172 : :
173 : 30667 : return length;
174 : : }
175 : :
176 : : /// Return true if the lower boundary is closer than the upper boundary
177 : : static bool
178 : 30667 : double_lower_boundary_is_closer(const double d)
179 : : {
180 : 30667 : const uint64_t rep = double_to_rep(d);
181 : 30667 : const uint64_t mant = rep & dbl_mant_mask;
182 : 30667 : const uint64_t expt = rep & dbl_expt_mask;
183 : 30667 : const bool is_subnormal = expt == 0;
184 : :
185 : : // True when f = 2^(p-1) (except for the smallest normal)
186 [ + + + + ]: 30667 : return !is_subnormal && mant == 0;
187 : : }
188 : :
189 : : /// Return an estimate for the decimal power of value, undershot by at most 1
190 : : static EXESS_I_NONBLOCKING int
191 : 30667 : approximate_power(const ExessSoftFloat value)
192 : : {
193 : : /* See "Printing Floating-Point Numbers Quickly and Accurately" by Robert
194 : : G. Burger and R. Kent Dybvig. The trick of undershooting by 0.69 comes
195 : : from the implementation by Ryan Juckett, see
196 : : https://ryanjuckett.com/printing-floating-point-numbers-part-2-dragon4/
197 : : */
198 : :
199 : 30667 : EXESS_I_CONSTEXPR double log10_2 = 0.30102999566398119521373889472449;
200 : :
201 : 30667 : const double f_msb_index = 64U - exess_clz64(value.f);
202 : 30667 : const double power = ceil(((f_msb_index + value.e) * log10_2) - 0.69);
203 : :
204 : 30667 : return (int)power;
205 : : }
206 : :
207 : : EXESS_I_NONBLOCKING ExessDigitCount
208 : 30667 : generate_digits(const double d, const size_t digits_size, char* const digits)
209 : : {
210 : : EXESS_DISABLE_CONVERSION_WARNINGS
211 [ + - + - : 30667 : assert(isfinite(d) && fpclassify(d) != FP_ZERO);
+ - + + -
+ - + ]
212 : : EXESS_RESTORE_WARNINGS
213 : :
214 : 30667 : const ExessSoftFloat value = soft_float_from_double(d);
215 : 30667 : const int power = approximate_power(value);
216 : 30667 : const bool is_even = !(value.f & 1U);
217 : 30667 : const bool lower_is_closer = double_lower_boundary_is_closer(d);
218 : :
219 : : // Calculate initial values so that v = (numer / denom) * 10^power
220 : : ExessBigint numer;
221 : : ExessBigint denom;
222 : : ExessBigint d_lower;
223 : 30667 : calculate_initial_values(
224 : 30667 : value.f, value.e, power, lower_is_closer, &numer, &denom, &d_lower);
225 : :
226 : : ExessBigint d_upper_storage;
227 : 30667 : ExessBigint* d_upper = NULL;
228 [ + + ]: 30667 : if (lower_is_closer) {
229 : : // Scale upper boundary to account for the closer lower boundary
230 : : // (the numerator and denominator were already scaled above)
231 : 2166 : d_upper_storage = d_lower;
232 : 2166 : d_upper = &d_upper_storage;
233 : 2166 : exess_bigint_shift_left(d_upper, 1);
234 : : } else {
235 : 28501 : d_upper = &d_lower; // Boundaries are the same, reuse the lower
236 : : }
237 : :
238 : : // Scale if necessary to make 1 <= (numer + delta) / denom < 10
239 : 30667 : ExessDigitCount count = {0, 0};
240 [ + + ]: 30667 : if (within_upper(&numer, &denom, d_upper, is_even)) {
241 : 17012 : count.expt = power;
242 : : } else {
243 : 13655 : count.expt = power - 1;
244 : 13655 : exess_bigint_multiply_u32(&numer, 10);
245 : 13655 : exess_bigint_multiply_u32(&d_lower, 10);
246 [ + + ]: 13655 : if (d_upper != &d_lower) {
247 : 1278 : exess_bigint_multiply_u32(d_upper, 10);
248 : : }
249 : : }
250 : :
251 : : // Write digits to output
252 [ - + ]: 30667 : assert(check_initial_values(&numer, &denom, d_upper));
253 : 30667 : count.count = emit_digits(
254 : : &numer, &denom, &d_lower, d_upper, is_even, digits_size, digits);
255 : :
256 : : // Trim trailing zeros
257 [ + + + + ]: 31818 : while (count.count > 1 && digits[count.count - 1] == '0') {
258 : 1151 : digits[--count.count] = 0;
259 : : }
260 : :
261 : 30667 : digits[count.count] = '\0';
262 : 30667 : return count;
263 : : }
|
