Branch data Line data Source code
1 : : // Copyright 2019-2026 David Robillard <d@drobilla.net>
2 : : // SPDX-License-Identifier: ISC
3 : :
4 : : #include <exess/exess.h>
5 : :
6 : : #include "big_decimal.h"
7 : : #include "bigint.h"
8 : : #include "decimal_to_double.h"
9 : : #include "ieee_float.h"
10 : : #include "int_math.h"
11 : : #include "macros.h"
12 : : #include "soft_float.h"
13 : :
14 : : #include <assert.h>
15 : : #include <float.h>
16 : : #include <math.h>
17 : : #include <stdbool.h>
18 : : #include <stdint.h>
19 : : #include <string.h>
20 : :
21 : : /// Number of digits always represented exactly by an unsigned 64-bit integer
22 : : static const int uint64_digits10 = 19;
23 : :
24 : : static EXESS_NONBLOCKING uint64_t
25 : 47260 : normalize(ExessSoftFloat* value, const uint64_t error)
26 : : {
27 : 47260 : const int original_e = value->e;
28 : :
29 : 47260 : *value = soft_float_normalize(*value);
30 : :
31 [ - + ]: 47260 : assert(value->e <= original_e);
32 : 47260 : return error << (unsigned)(original_e - value->e);
33 : : }
34 : :
35 : : /**
36 : : Return the error added by floating point multiplication.
37 : :
38 : : Should be l + r + l*r/(2^64) + 0.5, but we short the denominator to 63 due
39 : : to lack of precision, which effectively rounds up.
40 : : */
41 : : static inline uint64_t
42 : 23634 : product_error(const uint64_t lerror,
43 : : const uint64_t rerror,
44 : : const uint64_t half_ulp)
45 : : {
46 : 23634 : return lerror + rerror + ((lerror * rerror) >> 63U) + half_ulp;
47 : : }
48 : :
49 : : /**
50 : : Guess the binary floating point value for decimal input.
51 : :
52 : : @param significand Significand from the input.
53 : : @param expt10 Decimal exponent from the input.
54 : : @param n_digits Number of decimal digits in the significand.
55 : : @param[out] guess Either the exact number, or its predecessor.
56 : : @return True if `guess` is correct.
57 : : */
58 : : static EXESS_I_NONBLOCKING bool
59 : 23630 : sftod(const uint64_t significand,
60 : : const int expt10,
61 : : const int n_digits,
62 : : ExessSoftFloat* const guess)
63 : : {
64 [ - + ]: 23630 : assert(expt10 <= max_dec_expt);
65 [ - + ]: 23630 : assert(expt10 >= min_dec_expt);
66 : :
67 : : /* The general idea here is to try and find a power of 10 that we can
68 : : multiply by the significand to get the number. We get one from the
69 : : cache which is possibly too small, then multiply by another power of 10
70 : : to make up the difference if necessary. For example, with a target
71 : : power of 10^70, if we get 10^68 from the cache, then we multiply again
72 : : by 10^2. This, as well as normalization, accumulates error, which is
73 : : tracked throughout to know if we got the precise number. */
74 : :
75 : : // Use a common denominator of 2^3 to avoid fractions
76 : 23630 : EXESS_CONSTEXPR unsigned lg_denom = 3;
77 : 23630 : EXESS_CONSTEXPR uint64_t denom = 1U << 3U;
78 : 23630 : EXESS_CONSTEXPR uint64_t half_ulp = 4U;
79 : :
80 : : // Start out with just the significand, and no error
81 : 23630 : ExessSoftFloat input = {significand, 0};
82 : 23630 : uint64_t error = normalize(&input, 0);
83 : :
84 : : // Get a power of 10 that takes us most of the way without overshooting
85 : 23630 : int cached_expt10 = 0;
86 : 23630 : ExessSoftFloat pow10 = soft_float_pow10_under(expt10, &cached_expt10);
87 : :
88 : : // Get an exact fixup power if necessary
89 : 23630 : const int d_expt10 = expt10 - cached_expt10;
90 [ + + ]: 23630 : if (d_expt10) {
91 : 20553 : input = soft_float_multiply(input, soft_float_exact_pow10(d_expt10));
92 [ + + ]: 20553 : if (d_expt10 > uint64_digits10 - n_digits) {
93 : 9899 : error += half_ulp; // Product does not fit in an integer
94 : : }
95 : : }
96 : :
97 : : // Multiply the significand by the power, normalize, and update the error
98 : 23630 : input = soft_float_multiply(input, pow10);
99 : 23630 : error = normalize(&input, product_error(error, half_ulp, half_ulp));
100 : :
101 : : // Get the effective number of significant bits from the order of magnitude
102 : 23630 : const int magnitude = MIN(DBL_MANT_DIG, 64 + input.e - dbl_subnormal_expt);
103 : 23630 : const unsigned n_significant_bits = (unsigned)MAX(0, magnitude);
104 : :
105 : : // Calculate the number of "extra" bits of precision we have
106 [ - + ]: 23630 : assert(n_significant_bits <= 64);
107 : 23630 : unsigned n_extra_bits = 64U - n_significant_bits;
108 [ + + ]: 23630 : if (n_extra_bits + lg_denom >= 64U) {
109 : : // Very small subnormal where extra * denom does not fit in an integer
110 : : // Shift right (and accumulate some more error) to compensate
111 : 4 : const unsigned amount = (n_extra_bits + lg_denom) - 63;
112 : :
113 : 4 : input.f >>= amount;
114 : 4 : input.e += (int)amount;
115 : 4 : error = product_error((error >> amount) + 1U, half_ulp, half_ulp);
116 : 4 : n_extra_bits -= amount;
117 : : }
118 : :
119 : : // Calculate boundaries for the extra bits (with the common denominator)
120 [ - + ]: 23630 : assert(n_extra_bits < 64);
121 : 23630 : const uint64_t extra_mask = (1ULL << n_extra_bits) - 1U;
122 : 23630 : const uint64_t extra_bits = (input.f & extra_mask) * denom;
123 : 23630 : const uint64_t middle = (1ULL << (n_extra_bits - 1U)) * denom;
124 : 23630 : const uint64_t low = middle - error;
125 : 23630 : const uint64_t high = middle + error;
126 : :
127 : : // Round to nearest representable double
128 : 23630 : guess->f = (input.f >> n_extra_bits) + (extra_bits >= high);
129 : 23630 : guess->e = input.e + (int)n_extra_bits;
130 : :
131 : : // Too inaccurate if the extra bits are within the error around the middle
132 [ + + + + ]: 23630 : return extra_bits <= low || extra_bits >= high;
133 : : }
134 : :
135 : : static EXESS_I_NONBLOCKING int
136 : 29 : compare_buffer(const char* buf, const int expt, const ExessSoftFloat upper)
137 : : {
138 : : ExessBigint buf_bigint;
139 : 29 : exess_bigint_set_decimal_string(&buf_bigint, buf);
140 : :
141 : : ExessBigint upper_bigint;
142 : 29 : exess_bigint_set_u64(&upper_bigint, upper.f);
143 : :
144 [ + + ]: 29 : if (expt >= 0) {
145 : 23 : exess_bigint_multiply_pow10(&buf_bigint, (unsigned)expt);
146 : : } else {
147 : 6 : exess_bigint_multiply_pow10(&upper_bigint, (unsigned)-expt);
148 : : }
149 : :
150 [ + + ]: 29 : if (upper.e > 0) {
151 : 23 : exess_bigint_shift_left(&upper_bigint, (unsigned)upper.e);
152 : : } else {
153 : 6 : exess_bigint_shift_left(&buf_bigint, (unsigned)-upper.e);
154 : : }
155 : :
156 : 29 : return exess_bigint_compare(&buf_bigint, &upper_bigint);
157 : : }
158 : :
159 : : EXESS_PURE_FUNC static uint64_t
160 : 30653 : read_fraction(size_t n_digits, const char* const digits)
161 : : {
162 : 30653 : uint64_t frac = 0;
163 : :
164 [ + + ]: 434850 : for (unsigned i = 0U; i < n_digits; ++i) {
165 [ + - - + ]: 404197 : assert(digits[i] >= '0' && digits[i] <= '9');
166 : 404197 : frac = (frac * 10U) + (unsigned)(digits[i] - '0');
167 : : }
168 : :
169 : 30653 : return frac;
170 : : }
171 : :
172 : : EXESS_I_NONBLOCKING double
173 : 31651 : decimal_to_double(const BigDecimal in, const char* const digits)
174 : : {
175 : 31651 : EXESS_CONSTEXPR int n_exact_pow10 = sizeof(POW10) / sizeof(POW10[0]);
176 : :
177 : 31651 : EXESS_CONSTEXPR unsigned max_exact_int_digits = 15; // Digits that fit exactly
178 : 31651 : EXESS_CONSTEXPR int max_decimal_power = 309; // Max finite power
179 : 31651 : EXESS_CONSTEXPR int min_decimal_power = -324; // Min non-zero power
180 : :
181 : 31651 : EXESS_CONSTEXPR double special_values[] = {
182 : : (double)NAN, (double)-INFINITY, (double)INFINITY, -0.0, 0.0};
183 : :
184 [ + + ]: 31651 : if (in.kind < EXESS_NEGATIVE) {
185 : 995 : return special_values[in.kind]; // Special case
186 : : }
187 : :
188 [ + + ]: 30656 : const int sign = in.kind == EXESS_POSITIVE ? 1 : -1;
189 : :
190 [ + + ]: 30656 : if (in.expt > max_decimal_power) {
191 : 2 : return sign * (double)INFINITY; // Infinite magnitude => infinity
192 : : }
193 : :
194 [ + + ]: 30654 : if (in.expt < min_decimal_power) {
195 : 1 : return sign * 0.0; // Effectively zero
196 : : }
197 : :
198 : : // Represent value with integer mantissa and adjusted decimal power
199 : 30653 : const uint64_t frac = read_fraction(in.n_digits, digits);
200 : 30653 : const int power = in.expt - (int)in.n_digits;
201 : :
202 [ + + ]: 30653 : if (in.n_digits < max_exact_int_digits) { // Exact integer value
203 [ + + ]: 12777 : if (power < 0) {
204 [ + + ]: 6681 : if (-power < n_exact_pow10) {
205 : 3182 : return sign * ((double)frac / (double)POW10[-power]);
206 : : }
207 [ + + ]: 6096 : } else if (power < n_exact_pow10) {
208 : 3841 : return sign * ((double)frac * (double)POW10[power]);
209 : : }
210 : : }
211 : :
212 : : // Try to guess the number using only soft floating point (fast path)
213 : 23630 : ExessSoftFloat guess = {0, 0};
214 : 23630 : const bool exact = sftod(frac, power, (int)in.n_digits, &guess);
215 : 23630 : const double g = soft_float_to_double(guess);
216 [ + + ]: 23630 : if (exact) {
217 : 23601 : return sign * g;
218 : : }
219 : :
220 : : // Not sure, guess is either the number or its predecessor (rare slow path)
221 : : // Compare it with the buffer using bigints to find out which
222 : 29 : const ExessSoftFloat upper = {(guess.f * 2) + 1, guess.e - 1};
223 : 29 : const int cmp = compare_buffer(digits, power, upper);
224 [ + + + + : 29 : const bool round_up = (cmp > 0) || (cmp == 0 && (guess.f & 1U) != 0);
+ + ]
225 : :
226 [ + + ]: 29 : return sign * (round_up ? nextafter(g, (double)INFINITY) : g);
227 : : }
|
