8000 py/formatfloat: Calculate reference values with pow(). · dpwe/micropython@4923165 · GitHub
[go: up one dir, main page]
Skip to content

Commit 4923165

Browse files
dpwedpwe
committed
py/formatfloat: Calculate reference values with pow().
Formerly, ```formatfloat``` struggled with formatting exact powers of 10 as created by ```parsenum``` because of small differences in how those powers of 10 were calculated: ```formatfloat``` multiplied together a few values of 10^(2^Y) from a lookup table, but ```parsenum``` used ```pow(10, X)```. This was mitigated in micropython#8905 by adding 2eps of extra margin when comparing values ("aggressive rounding up"). This patch instead eliminates the discrepency by using ```pow()``` in ```formatfloat``` also. This eliminates the need for the 2eps margin, as well as the lookup tables. It is likely slower to run, however. In addition, this patch directly estimates the power-of-10 exponent from the power-of-2 exponent in the floating-point representation. This change surfaced a previously-noted problem with parsing numbers including trailing zeros. The fix to that problem, micropython#8980, is included in this PR to ensure the tests pass. Signed-off-by: Dan Ellis <dan.ellis@gmail.com>
1 parent 963e599 commit 4923165

File tree

4 files changed

+73
-93
lines changed

4 files changed

+73
-93
lines changed

‎py/formatfloat.c

Lines changed: 35 additions & 84 deletions
Original file line numberDiff line numberDiff line change
@@ -65,23 +65,19 @@
6565
#define FLT_EXP_MASK 0x7F800000
6666
#define FLT_MAN_MASK 0x007FFFFF
6767

68-
union floatbits {
69-
float f;
70-
uint32_t u;
71-
};
7268
static inline int fp_signbit(float x) {
73-
union floatbits fb = {x};
74-
return fb.u & FLT_SIGN_MASK;
69+
mp_float_union_t fb = {x};
70+
return fb.i & FLT_SIGN_MASK;
7571
}
7672
#define fp_isnan(x) isnan(x)
7773
#define fp_isinf(x) isinf(x)
7874
static inline int fp_iszero(float x) {
79-
union floatbits fb = {x};
80-
return fb.u == 0;
75+
mp_float_union_t fb = {x};
76+
return fb.i == 0;
8177
}
8278
static inline int fp_isless1(float x) {
83-
union floatbits fb = {x};
84-
return fb.u < 0x3f800000;
79+
mp_float_union_t fb = {x};
80+
return fb.i < 0x3f800000;
8581
}
8682

8783
#elif MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_DOUBLE
@@ -99,28 +95,6 @@ static inline int fp_isless1(float x) {
9995

10096
#endif // MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_FLOAT/DOUBLE
10197

102-
static inline int fp_ge_eps(FPTYPE x, FPTYPE y) {
103-
mp_float_union_t fb_y = {y};
104-
// Back off 2 eps.
105-
// This is valid for almost all values, but in practice
106-
// it's only used when y = 1eX for X>=0.
107-
fb_y.i -= 2;
108-
return x >= fb_y.f;
109-
}
110-
111-
static const FPTYPE g_pos_pow[] = {
112-
#if FPDECEXP > 32
113-
MICROPY_FLOAT_CONST(1e256), MICROPY_FLOAT_CONST(1e128), MICROPY_FLOAT_CONST(1e64),
114-
#endif
115-
MICROPY_FLOAT_CONST(1e32), MICROPY_FLOAT_CONST(1e16), MICROPY_FLOAT_CONST(1e8), MICROPY_FLOAT_CONST(1e4), MICROPY_FLOAT_CONST(1e2), MICROPY_FLOAT_CONST(1e1)
116-
};
117-
static const FPTYPE g_neg_pow[] = {
118-
#if FPDECEXP > 32
119-
MICROPY_FLOAT_CONST(1e-256), MICROPY_FLOAT_CONST(1e-128), MICROPY_FLOAT_CONST(1e-64),
120-
#endif
121-
MICROPY_FLOAT_CONST(1e-32), MICROPY_FLOAT_CONST(1e-16), MICROPY_FLOAT_CONST(1e-8), MICROPY_FLOAT_CONST(1e-4), MICROPY_FLOAT_CONST(1e-2), MICROPY_FLOAT_CONST(1e-1)
122-
};
123-
12498
int mp_format_float(FPTYPE f, char *buf, size_t buf_size, char fmt, int prec, char sign) {
12599

126100
char *s = buf;
@@ -177,14 +151,20 @@ int mp_format_float(FPTYPE f, char *buf, size_t buf_size, char fmt, int prec, ch
177151
if (fmt == 'g' && prec == 0) {
178152
prec = 1;
179153
}
180-
int e, e1;
154+
int e;
181155
int dec = 0;
182156
char e_sign = '\0';
183157
int num_digits = 0;
184-
const FPTYPE *pos_pow = g_pos_pow;
185-
const FPTYPE *neg_pow = g_neg_pow;
186158
int signed_e = 0;
187159

160+
// Approximate (?) power of 10 exponent from binary exponent.
161+
// It seems to be right all the time, but to be sure, we try several values
162+
// of 10 from one smaller.
163+
mp_float_union_t fb = {f};
164+
FPTYPE e_float = (fb.p.exp - MP_FLOAT_EXP_OFFSET) *
165+
FPCONST(0.3010299956639812); // 1/log2(10).
166+
int e_guess = (int)((FPCONST(0.5) * (1 - 2 * (e_float < 0))) + e_float);
167+
188168
if (fp_iszero(f)) {
189169
e = 0;
190170
if (fmt == 'f') {
@@ -203,24 +183,25 @@ int mp_format_float(FPTYPE f, char *buf, size_t buf_size, char fmt, int prec, ch
203183
}
204184
}
205185
} else if (fp_isless1(f)) {
206-
FPTYPE f_mod = f;
186+
FPTYPE f_entry = f; // Save f in case we go to 'f' format.
207187
// Build negative exponent
208-
for (e = 0, e1 = FPDECEXP; e1; e1 >>= 1, pos_pow++, neg_pow++) {
209-
if (*neg_pow > f_mod) {
210-
e += e1;
211-
f_mod *= *pos_pow;
212-
}
188+
e = -e_guess - 1;
189+
FPTYPE u_base = MICROPY_FLOAT_C_FUN(pow)(10, -e);
190+
while (u_base > f) {
191+
++e;
192+
u_base = MICROPY_FLOAT_C_FUN(pow)(10, -e);
213193
}
194+
f *= MICROPY_FLOAT_C_FUN(pow)(10, e);
214195

215196
char e_sign_char = '-';
216-
if (fp_isless1(f_mod) && f_mod >= FPROUND_TO_ONE) {
217-
f_mod = FPCONST(1.0);
197+
if (fp_isless1(f) && f >= FPROUND_TO_ONE) {
198+
f = FPCONST(1.0);
218199
if (e == 0) {
219200
e_sign_char = '+';
220201
}
221-
} else if (fp_isless1(f_mod)) {
202+
} else if (fp_isless1(f)) {
222203
e++;
223-
f_mod *= FPCONST(10.0);
204+
f *= FPCONST(10.0);
224205
}
225206

226207
// If the user specified 'g' format, and e is <= 4, then we'll switch
@@ -241,6 +222,7 @@ int mp_format_float(FPTYPE f, char *buf, size_t buf_size, char fmt, int prec, ch
241222

242223
num_digits = prec;
243224
signed_e = 0;
225+
f = f_entry;
244226
++num_digits;
245227
} else {
246228
// For e & g formats, we'll be printing the exponent, so set the
@@ -262,19 +244,11 @@ int mp_format_float(FPTYPE f, char *buf, size_t buf_size, char fmt, int prec, ch
262244
// scaling it. Instead, we find the product of powers of 10
263245
// that is not greater than it, and use that to start the
264246
// mantissa.
265-
FPTYPE u_base = FPCONST(1.0);
266-
for (e = 0, e1 = FPDECEXP; e1; e1 >>= 1, pos_pow++) {
267-
FPTYPE next_u = u_base * *pos_pow;
268-
// fp_ge_eps performs "f >= (next_u - 2eps)" so that if, for
269-
// numerical reasons, f is very close to a power of ten but
270-
// not strictly equal, we still treat it as that power of 10.
271-
// The comparison was failing for maybe 10% of 1eX values, but
272-
// although rounding fixed many of them, there were still some
273-
// rendering as 9.99999998e(X-1).
274-
if (fp_ge_eps(f, next_u)) {
275-
u_base = next_u;
276-
e += e1;
277-
}
247+
e = e_guess - 1;
248+
FPTYPE next_u = MICROPY_FLOAT_C_FUN(pow)(10, e + 1);
249+
while (f >= next_u) {
250+
++e;
251+
next_u = MICROPY_FLOAT_C_FUN(pow)(10, e + 1);
278252
}
279253

280254
// If the user specified fixed format (fmt == 'f') and e makes the
@@ -341,39 +315,16 @@ int mp_format_float(FPTYPE f, char *buf, size_t buf_size, char fmt, int prec, ch
341315
num_digits = prec;
342316
}
343317

344-
if (signed_e < 0) {
345-
// The algorithm below treats numbers smaller than 1 by scaling them
346-
// repeatedly by 10 to bring the new digit to the top. Our input number
347-
// was smaller than 1, so scale it up to be 1 <= f < 10.
348-
FPTYPE u_base = FPCONST(1.0);
349-
const FPTYPE *pow_u = g_pos_pow;
350-
for (int m = FPDECEXP; m; m >>= 1, pow_u++) {
351-
if (m & e) {
352-
u_base *= *pow_u;
353-
}
354-
}
355-
f *= u_base;
356-
}
357-
358318
int d = 0;
359319
int num_digits_left = num_digits;
360320
for (int digit_index = signed_e; num_digits_left >= 0; --digit_index) {
361321
FPTYPE u_base = FPCONST(1.0);
362322
if (digit_index > 0) {
363323
// Generate 10^digit_index for positive digit_index.
364-
const FPTYPE *pow_u = g_pos_pow;
365-
int target_index = digit_index;
366-
for (int m = FPDECEXP; m; m >>= 1, pow_u++) {
367-
if (m & target_index) {
368-
u_base *= *pow_u;
369-
}
370-
}
324+
u_base = MICROPY_FLOAT_C_FUN(pow)(10, digit_index);
371325
}
372326
for (d = 0; d < 9; ++d) {
373-
// This is essentially "if (f < u_base)", but with 2eps margin
374-
// so that if f is just a tiny bit smaller, we treat it as
375-
// equal (and accept the additional digit value).
376-
if (!fp_ge_eps(f, u_base)) {
327+
if (f < u_base) {
377328
break;
378329
}
379330
f -= u_base;
@@ -390,7 +341,7 @@ int mp_format_float(FPTYPE f, char *buf, size_t buf_size, char fmt, int prec, ch
390341
--dec;
391342
--num_digits_left;
392343
if (digit_index <= 0) {
393-
// Once we get below 1.0, we scale up f instead of calculting
344+
// Once we get below 1.0, we scale up f instead of calculating
394345
// negative powers of 10 in u_base. This provides better
395346
// renditions of exact decimals like 1/16 etc.
396347
f *= FPCONST(10.0);

‎py/misc.h

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -243,10 +243,12 @@ extern mp_uint_t mp_verbose_flag;
243243

244244
#if MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_DOUBLE
245245
#define MP_FLOAT_EXP_BITS (11)
246+
#define MP_FLOAT_EXP_OFFSET (1024)
246247
#define MP_FLOAT_FRAC_BITS (52)
247248
typedef uint64_t mp_float_uint_t;
248249
#elif MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_FLOAT
249250
#define MP_FLOAT_EXP_BITS (8)
251+
#define MP_FLOAT_EXP_OFFSET (128)
250252
#define MP_FLOAT_FRAC_BITS (23)
251253
typedef uint32_t mp_float_uint_t;
252254
#endif

‎py/parsenum.c

Lines changed: 33 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -255,6 +255,7 @@ mp_obj_t mp_parse_num_float(const char *str, size_t len, bool allow_imag, mp_lex
255255
bool exp_neg = false;
256256
int exp_val = 0;
257257
int exp_extra = 0;
258+
int trailing_zeros = 0, trailing_zeros_frac = 0;
258259
while (str < top) {
259260
unsigned int dig = *str++;
260261
if ('0' <= dig && dig <= '9') {
@@ -267,17 +268,40 @@ mp_obj_t mp_parse_num_float(const char *str, size_t len, bool allow_imag, mp_lex
267268
exp_val = 10 * exp_val + dig 93C6 ;
268269
}
269270
} else {
270-
if (dec_val < DEC_VAL_MAX) {
271-
// dec_val won't overflow so keep accumulating
272-
dec_val = 10 * dec_val + dig;
271+
if (dig == 0) {
272+
// Defer treatment of zeros in fractional part. If nothing comes afterwards, ignore them.
273+
++trailing_zeros;
273274
if (in == PARSE_DEC_IN_FRAC) {
274-
--exp_extra;
275+
++trailing_zeros_frac;
275276
}
276277
} else {
277-
// dec_val might overflow and we anyway can't represent more digits
278-
// of precision, so ignore the digit and just adjust the exponent
279-
if (in == PARSE_DEC_IN_INTG) {
280-
++exp_extra;
278+
// Time to un-defer any trailing zeros.
279+
if (trailing_zeros) {
280+
while (trailing_zeros) {
281+
if (dec_val < DEC_VAL_MAX) {
282+
dec_val = 10 * dec_val;
283+
if (trailing_zeros_frac) {
284+
--exp_extra;
285+
--trailing_zeros_frac;
286+
}
287+
}
288+
--trailing_zeros;
289+
}
290+
trailing_zeros_frac = 0;
291+
}
292+
293+
if (dec_val < DEC_VAL_MAX) {
294+
// dec_val won't overflow so keep accumulating
295+
dec_val = 10 * dec_val + dig;
296+
if (in == PARSE_DEC_IN_FRAC) {
297+
--exp_extra;
298+
}
299+
} else {
300+
// dec_val might overflow and we anyway can't represent more digits
301+
// of precision, so ignore the digit and just adjust the exponent
302+
if (in == PARSE_DEC_IN_INTG) {
303+
++exp_extra;
304+
}
281305
}
282306
}
283307
}
@@ -311,7 +335,7 @@ mp_obj_t mp_parse_num_float(const char *str, size_t len, bool allow_imag, mp_lex
311335
}
312336

313337
// apply the exponent, making sure it's not a subnormal value
314-
exp_val += exp_extra;
338+
exp_val += (exp_extra + trailing_zeros - trailing_zeros_frac);
315339
if (exp_val < SMALL_NORMAL_EXP) {
316340
exp_val -= SMALL_NORMAL_EXP;
317341
dec_val *= SMALL_NORMAL_VAL;

‎tests/float/float_format_ints_doubleprec.py

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -13,3 +13,6 @@
1313
v2 = 0x6974E718D7D7625A # 1e200
1414
print("{:.12e}".format(array.array("d", v1.to_bytes(8, sys.byteorder))[0]))
1515
print("{:.12e}".format(array.array("d", v2.to_bytes(8, sys.byteorder))[0]))
16+
17+
for i in range(300):
18+
print(float("1e" + str(i)))

0 commit comments

Comments
 (0)
0