2
* Copyright (c) 2012 Adam Hraska
5
* Redistribution and use in source and binary forms, with or without
6
* modification, are permitted provided that the following conditions
9
* - Redistributions of source code must retain the above copyright
10
* notice, this list of conditions and the following disclaimer.
11
* - Redistributions in binary form must reproduce the above copyright
12
* notice, this list of conditions and the following disclaimer in the
13
* documentation and/or other materials provided with the distribution.
14
* - The name of the author may not be used to endorse or promote products
15
* derived from this software without specific prior written permission.
17
* THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
18
* IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
19
* OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
20
* IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
21
* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
22
* NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
23
* DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
24
* THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
25
* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
26
* THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
29
#include <double_to_str.h>
31
#include "private/power_of_ten.h"
32
#include <ieee_double.h>
39
* Floating point numbers are converted from their binary representation
40
* into a decimal string using the algorithm described in:
41
* Printing floating-point numbers quickly and accurately with integers
45
/** The computation assumes a significand of 64 bits. */
46
static const int significand_width = 64;
48
/* Scale exponents to interval [alpha, gamma] to simplify conversion. */
49
static const int alpha = -59;
50
static const int gamma = -32;
53
/** Returns true if the most-significant bit of num.significand is set. */
54
static bool is_normalized(fp_num_t num)
56
assert(8*sizeof(num.significand) == significand_width);
58
/* Normalized == most significant bit of the significand is set. */
59
return (num.significand & (1ULL << (significand_width - 1))) != 0;
62
/** Returns a normalized num with the MSbit set. */
63
static fp_num_t normalize(fp_num_t num)
65
const uint64_t top10bits = 0xffc0000000000000ULL;
67
/* num usually comes from ieee_double with top 10 bits zero. */
68
while (0 == (num.significand & top10bits)) {
69
num.significand <<= 10;
73
while (!is_normalized(num)) {
74
num.significand <<= 1;
82
/** Returns x * y with an error of less than 0.5 ulp. */
83
static fp_num_t multiply(fp_num_t x, fp_num_t y)
85
assert(/* is_normalized(x) && */ is_normalized(y));
87
const uint32_t low_bits = -1;
90
a = x.significand >> 32;
91
b = x.significand & low_bits;
92
c = y.significand >> 32;
93
d = y.significand & low_bits;
95
uint64_t bd, ad, bc, ac;
102
/* Denote 32 bit parts of x a y as: x == a b, y == c d. Then:
106
* ad bd .. multiplication of 32bit parts results in 64bit parts
109
* [b|d] .. Depicts 64 bit intermediate results and how
110
* [a|d] the 32 bit parts of these results overlap and
111
* [b|c] contribute to the final result.
117
uint64_t tmp = (bd >> 32) + (ad & low_bits) + (bc & low_bits);
123
ret.significand = ac + (bc >> 32) + (ad >> 32) + (tmp >> 32);
124
ret.exponent = x.exponent + y.exponent + significand_width;
130
/** Returns a - b. Both must have the same exponent. */
131
static fp_num_t subtract(fp_num_t a, fp_num_t b)
133
assert(a.exponent == b.exponent);
134
assert(a.significand >= b.significand);
138
result.significand = a.significand - b.significand;
139
result.exponent = a.exponent;
145
/** Returns the interval [low, high] of numbers that convert to binary val. */
146
static void get_normalized_bounds(ieee_double_t val, fp_num_t *high,
147
fp_num_t *low, fp_num_t *val_dist)
150
* Only works if val comes directly from extract_ieee_double without
151
* being manipulated in any way (eg it must not be normalized).
153
assert(!is_normalized(val.pos_val));
155
high->significand = (val.pos_val.significand << 1) + 1;
156
high->exponent = val.pos_val.exponent - 1;
158
/* val_dist = high - val */
159
val_dist->significand = 1;
160
val_dist->exponent = val.pos_val.exponent - 1;
162
/* Distance from both lower and upper bound is the same. */
163
if (!val.is_accuracy_step) {
164
low->significand = (val.pos_val.significand << 1) - 1;
165
low->exponent = val.pos_val.exponent - 1;
167
low->significand = (val.pos_val.significand << 2) - 1;
168
low->exponent = val.pos_val.exponent - 2;
171
*high = normalize(*high);
174
* Lower bound may not be normalized if subtracting 1 unit
175
* reset the most-significant bit to 0.
177
low->significand = low->significand << (low->exponent - high->exponent);
178
low->exponent = high->exponent;
180
val_dist->significand =
181
val_dist->significand << (val_dist->exponent - high->exponent);
182
val_dist->exponent = high->exponent;
185
/** Determines the interval of numbers that have the binary representation
188
* Numbers in the range [scaled_upper_bound - bounds_delta, scaled_upper_bound]
189
* have the same double binary representation as val.
191
* Bounds are scaled by 10^scale so that alpha <= exponent <= gamma.
192
* Moreover, scaled_upper_bound is normalized.
194
* val_dist is the scaled distance from val to the upper bound, ie
195
* val_dist == (upper_bound - val) * 10^scale
197
static void calc_scaled_bounds(ieee_double_t val, fp_num_t *scaled_upper_bound,
198
fp_num_t *bounds_delta, fp_num_t *val_dist, int *scale)
200
fp_num_t upper_bound, lower_bound;
202
get_normalized_bounds(val, &upper_bound, &lower_bound, val_dist);
204
assert(upper_bound.exponent == lower_bound.exponent);
205
assert(is_normalized(upper_bound));
206
assert(normalize(val.pos_val).exponent == upper_bound.exponent);
209
* Find such a cached normalized power of 10 that if multiplied
210
* by upper_bound the binary exponent of upper_bound almost vanishes,
212
* upper_scaled := upper_bound * 10^scale
213
* alpha <= upper_scaled.exponent <= gamma
214
* alpha <= upper_bound.exponent + pow_10.exponent + 64 <= gamma
216
fp_num_t scaling_power_of_10;
217
int lower_bin_exp = alpha - upper_bound.exponent - significand_width;
219
get_power_of_ten(lower_bin_exp, &scaling_power_of_10, scale);
221
int scale_exp = scaling_power_of_10.exponent;
222
assert(alpha <= upper_bound.exponent + scale_exp + significand_width);
223
assert(upper_bound.exponent + scale_exp + significand_width <= gamma);
225
fp_num_t upper_scaled = multiply(upper_bound, scaling_power_of_10);
226
fp_num_t lower_scaled = multiply(lower_bound, scaling_power_of_10);
227
*val_dist = multiply(*val_dist, scaling_power_of_10);
229
assert(alpha <= upper_scaled.exponent && upper_scaled.exponent <= gamma);
232
* Any value between lower and upper bound would be represented
233
* in binary as the double val originated from. The bounds were
234
* however scaled by an imprecise power of 10 (error less than
235
* 1 ulp) so the scaled bounds have an error of less than 1 ulp.
236
* Conservatively round the lower bound up and the upper bound
237
* down by 1 ulp just to be on the safe side. It avoids pronouncing
238
* produced decimal digits as correct if such a decimal number
239
* is close to the bounds to within 1 ulp.
241
upper_scaled.significand -= 1;
242
lower_scaled.significand += 1;
244
*bounds_delta = subtract(upper_scaled, lower_scaled);
245
*scaled_upper_bound = upper_scaled;
249
/** Rounds the last digit of buf so that it is closest to the converted number.*/
250
static void round_last_digit(uint64_t rest, uint64_t w_dist, uint64_t delta,
251
uint64_t digit_val_diff, char *buf, int len)
254
* | <------- delta -------> |
255
* | | <---- w_dist ----> |
262
* delta = upper - lower .. conservative/safe interval
264
* upper = "number represented by digits in buf" + rest
266
* Changing buf[len - 1] changes the value represented by buf
267
* by digit_val_diff * scaling, where scaling is shared by
272
/* Current number in buf is greater than the double being converted */
273
bool cur_greater_w = rest < w_dist;
274
/* Rounding down by one would keep buf in between bounds (in safe rng). */
275
bool next_in_val_rng = cur_greater_w && (rest + digit_val_diff < delta);
276
/* Rounding down by one would bring buf closer to the processed number. */
277
bool next_closer = next_in_val_rng
278
&& (rest + digit_val_diff < w_dist || rest - w_dist < w_dist - rest);
280
/* Of the shortest strings pick the one that is closest to the actual
281
floating point number. */
282
while (next_closer) {
283
assert('0' < buf[len - 1]);
284
assert(0 < digit_val_diff);
287
rest += digit_val_diff;
289
cur_greater_w = rest < w_dist;
290
next_in_val_rng = cur_greater_w && (rest + digit_val_diff < delta);
291
next_closer = next_in_val_rng
292
&& (rest + digit_val_diff < w_dist || rest - w_dist < w_dist - rest);
297
/** Generates the shortest accurate decimal string representation.
299
* Outputs (mostly) the shortest accurate string representation
300
* for the number scaled_upper - val_dist. Numbers in the interval
301
* [scaled_upper - delta, scaled_upper] have the same binary
302
* floating point representation and will therefore share the
303
* shortest string representation (up to the rounding of the last
304
* digit to bring the shortest string also the closest to the
307
* @param scaled_upper Scaled upper bound of numbers that have the
308
* same binary representation as the converted number.
309
* Scaled by 10^-scale so that alpha <= exponent <= gamma.
310
* @param delta scaled_upper - delta is the lower bound of numbers
311
* that share the same binary representation in double.
312
* @param val_dist scaled_upper - val_dist is the number whose
313
* decimal string we're generating.
314
* @param scale Decimal scaling of the value to convert (ie scaled_upper).
315
* @param buf Buffer to store the string representation. Must be large
316
* enough to store all digits and a null terminator. At most
317
* MAX_DOUBLE_STR_LEN digits will be written (not counting
318
* the null terminator).
319
* @param buf_size Size of buf in bytes.
320
* @param dec_exponent Will be set to the decimal exponent of the number
323
* @return Number of digits; negative on failure (eg buffer too small).
325
static int gen_dec_digits(fp_num_t scaled_upper, fp_num_t delta,
326
fp_num_t val_dist, int scale, char *buf, size_t buf_size, int *dec_exponent)
329
* The integral part of scaled_upper is 5 to 32 bits long while
330
* the remaining fractional part is 59 to 32 bits long because:
331
* -59 == alpha <= scaled_upper.e <= gamma == -32
333
* | <------- delta -------> |
334
* | | <--- val_dist ---> |
335
* | | |<- remainder ->|
342
assert(scaled_upper.significand != 0);
343
assert(alpha <= scaled_upper.exponent && scaled_upper.exponent <= gamma);
344
assert(scaled_upper.exponent == delta.exponent);
345
assert(scaled_upper.exponent == val_dist.exponent);
346
assert(val_dist.significand <= delta.significand);
348
/* We'll produce at least one digit and a null terminator. */
353
/* one is number 1 encoded with the same exponent as scaled_upper */
355
one.significand = ((uint64_t) 1) << (-scaled_upper.exponent);
356
one.exponent = scaled_upper.exponent;
359
* Extract the integral part of scaled_upper.
360
* upper / one == upper >> -one.e
362
uint32_t int_part = (uint32_t)(scaled_upper.significand >> (-one.exponent));
365
* Fractional part of scaled_upper.
366
* upper % one == upper & (one.f - 1)
368
uint64_t frac_part = scaled_upper.significand & (one.significand - 1);
371
* The integral part of upper has at least 5 bits (64 + alpha) and
372
* at most 32 bits (64 + gamma). The integral part has at most 10
373
* decimal digits, so kappa <= 10.
376
uint32_t div = 1000000000;
379
/* Produce decimal digits for the integral part of upper. */
381
int digit = int_part / div;
386
/* Skip leading zeros. */
387
if (digit != 0 || len != 0) {
388
/* Current length + new digit + null terminator <= buf_size */
389
if (len + 2 <= buf_size) {
390
buf[len] = '0' + digit;
398
* Difference between the so far produced decimal number and upper
399
* is calculated as: remaining_int_part * one + frac_part
401
uint64_t remainder = (((uint64_t)int_part) << -one.exponent) + frac_part;
403
/* The produced decimal number would convert back to upper. */
404
if (remainder <= delta.significand) {
405
assert(0 < len && len < buf_size);
406
*dec_exponent = kappa - scale;
409
/* Of the shortest representations choose the numerically closest. */
410
round_last_digit(remainder, val_dist.significand, delta.significand,
411
(uint64_t)div << (-one.exponent), buf, len);
418
/* Generate decimal digits for the fractional part of upper. */
421
* Does not overflow because at least 5 upper bits were
422
* taken by the integral part and are now unused in frac_part.
425
delta.significand *= 10;
426
val_dist.significand *= 10;
428
/* frac_part / one */
429
int digit = (int)(frac_part >> (-one.exponent));
431
/* frac_part %= one */
432
frac_part &= one.significand - 1;
436
/* Skip leading zeros. */
437
if (digit == 0 && len == 0) {
441
/* Current length + new digit + null terminator <= buf_size */
442
if (len + 2 <= buf_size) {
443
buf[len] = '0' + digit;
448
} while (frac_part > delta.significand);
450
assert(0 < len && len < buf_size);
452
*dec_exponent = kappa - scale;
455
/* Of the shortest representations choose the numerically closest one. */
456
round_last_digit(frac_part, val_dist.significand, delta.significand,
457
one.significand, buf, len);
462
/** Produce a string for 0.0 */
463
static int zero_to_str(char *buf, size_t buf_size, int *dec_exponent)
476
/** Converts a non-special double into its shortest accurate string
479
* Produces an accurate string representation, ie the string will
480
* convert back to the same binary double (eg via strtod). In the
481
* vast majority of cases (99%) the string will be the shortest such
482
* string that is also the closest to the value of any shortest
483
* string representations. Therefore, no trailing zeros are ever
486
* Conceptually, the value is: buf * 10^dec_exponent
488
* Never outputs trailing zeros.
490
* @param ieee_val Binary double description to convert. Must be the product
491
* of extract_ieee_double and it must not be a special number.
492
* @param buf Buffer to store the string representation. Must be large
493
* enough to store all digits and a null terminator. At most
494
* MAX_DOUBLE_STR_LEN digits will be written (not counting
495
* the null terminator).
496
* @param buf_size Size of buf in bytes.
497
* @param dec_exponent Will be set to the decimal exponent of the number
500
* @return The number of printed digits. A negative value indicates
501
* an error: buf too small (or ieee_val.is_special).
503
int double_to_short_str(ieee_double_t ieee_val, char *buf, size_t buf_size,
506
/* The whole computation assumes 64bit significand. */
507
assert(sizeof(ieee_val.pos_val.significand) == sizeof(uint64_t));
509
if (ieee_val.is_special) {
513
/* Zero cannot be normalized. Handle it here. */
514
if (0 == ieee_val.pos_val.significand) {
515
return zero_to_str(buf, buf_size, dec_exponent);
518
fp_num_t scaled_upper_bound;
523
calc_scaled_bounds(ieee_val, &scaled_upper_bound,
524
&delta, &val_dist, &scale);
526
int len = gen_dec_digits(scaled_upper_bound, delta, val_dist, scale,
527
buf, buf_size, dec_exponent);
529
assert(len <= MAX_DOUBLE_STR_LEN);
533
/** Generates a fixed number of decimal digits of w_scaled.
535
* double == w_scaled * 10^-scale, where alpha <= w_scaled.e <= gamma
537
* @param w_scaled Scaled number by 10^-scale so that
538
* alpha <= exponent <= gamma
539
* @param scale Decimal scaling of the value to convert (ie w_scaled).
540
* @param signif_d_cnt Maximum number of significant digits to output.
541
* Negative if as many as possible are requested.
542
* @param frac_d_cnt Maximum number of fractional digits to output.
543
* Negative if as many as possible are requested.
544
* Eg. if 2 then 1.234 -> "1.23"; if 2 then 3e-9 -> "0".
545
* @param buf Buffer to store the string representation. Must be large
546
* enough to store all digits and a null terminator. At most
547
* MAX_DOUBLE_STR_LEN digits will be written (not counting
548
* the null terminator).
549
* @param buf_size Size of buf in bytes.
551
* @return Number of digits; negative on failure (eg buffer too small).
553
static int gen_fixed_dec_digits(fp_num_t w_scaled, int scale, int signif_d_cnt,
554
int frac_d_cnt, char *buf, size_t buf_size, int *dec_exponent)
556
/* We'll produce at least one digit and a null terminator. */
557
if (0 == signif_d_cnt || buf_size < 2) {
562
* The integral part of w_scaled is 5 to 32 bits long while the
563
* remaining fractional part is 59 to 32 bits long because:
564
* -59 == alpha <= w_scaled.e <= gamma == -32
567
* | 5..32 bits | 32..59 bits | == w_scaled == w * 10^scale
568
* | int_part | frac_part |
569
* |0 0 .. 0 1|0 0 .. 0 0| == one == 1.0
570
* | 0 |0 0 .. 0 1| == w_err == 1 * 2^w_scaled.e
572
assert(alpha <= w_scaled.exponent && w_scaled.exponent <= gamma);
573
assert(0 != w_scaled.significand);
576
* Scaling the number being converted by 10^scale introduced
577
* an error of less that 1 ulp. The actual value of w_scaled
578
* could lie anywhere between w_scaled.signif +/- w_err.
579
* Scale the error locally as we scale the fractional part
584
/* one is number 1.0 encoded with the same exponent as w_scaled */
586
one.significand = ((uint64_t) 1) << (-w_scaled.exponent);
587
one.exponent = w_scaled.exponent;
589
/* Extract the integral part of w_scaled.
590
w_scaled / one == w_scaled >> -one.e */
591
uint32_t int_part = (uint32_t)(w_scaled.significand >> (-one.exponent));
593
/* Fractional part of w_scaled.
594
w_scaled % one == w_scaled & (one.f - 1) */
595
uint64_t frac_part = w_scaled.significand & (one.significand - 1);
599
* The integral part of w_scaled has at least 5 bits (64 + alpha = 5)
600
* and at most 32 bits (64 + gamma = 32). The integral part has
601
* at most 10 decimal digits, so kappa <= 10.
604
uint32_t div = 1000000000;
606
int rem_signif_d_cnt = signif_d_cnt;
608
(frac_d_cnt >= 0) ? (kappa - scale + frac_d_cnt) : INT_MAX;
610
/* Produce decimal digits for the integral part of w_scaled. */
611
while (kappa > 0 && rem_signif_d_cnt != 0 && rem_frac_d_cnt > 0) {
612
int digit = int_part / div;
619
/* Skip leading zeros. */
620
if (digit == 0 && len == 0) {
624
/* Current length + new digit + null terminator <= buf_size */
625
if (len + 2 <= buf_size) {
626
buf[len] = '0' + digit;
634
/* Generate decimal digits for the fractional part of w_scaled. */
635
while (w_err <= frac_part && rem_signif_d_cnt != 0 && rem_frac_d_cnt > 0) {
637
* Does not overflow because at least 5 upper bits were
638
* taken by the integral part and are now unused in frac_part.
643
/* frac_part / one */
644
int digit = (int)(frac_part >> (-one.exponent));
646
/* frac_part %= one */
647
frac_part &= one.significand - 1;
652
/* Skip leading zeros. */
653
if (digit == 0 && len == 0) {
657
/* Current length + new digit + null terminator <= buf_size */
658
if (len + 2 <= buf_size) {
659
buf[len] = '0' + digit;
667
assert(/* 0 <= len && */ len < buf_size);
670
*dec_exponent = kappa - scale;
671
assert(frac_d_cnt < 0 || -frac_d_cnt <= *dec_exponent);
674
* The number of fractional digits was too limiting to produce
677
assert(rem_frac_d_cnt <= 0 || w_scaled.significand == 0);
683
if (len < buf_size) {
685
assert(signif_d_cnt < 0 || (int)len <= signif_d_cnt);
693
/** Converts a non-special double into its string representation.
695
* Conceptually, the truncated double value is: buf * 10^dec_exponent
697
* Conversion errors are tracked, so all produced digits except the
698
* last one are accurate. Garbage digits are never produced.
699
* If the requested number of digits cannot be produced accurately
700
* due to conversion errors less digits are produced than requested
701
* and the last digit has an error of +/- 1 (so if '7' is the last
702
* converted digit it might have been converted to any of '6'..'8'
703
* had the conversion been completely precise).
705
* If no error occurs at least one digit is output.
707
* The conversion stops once the requested number of significant or
708
* fractional digits is reached or the conversion error is too large
709
* to generate any more digits (whichever happens first).
711
* Any digits following the first (most-significant) digit (this digit
712
* included) are counted as significant digits; eg:
713
* 1.4, 4 signif -> "1400" * 10^-3, ie 1.400
714
* 1000.3, 1 signif -> "1" * 10^3 ie 1000
715
* 0.003, 2 signif -> "30" * 10^-4 ie 0.003
716
* 9.5 1 signif -> "9" * 10^0, ie 9
718
* Any digits following the decimal point are counted as fractional digits.
719
* This includes the zeros that would appear between the decimal point
720
* and the first non-zero fractional digit. If fewer fractional digits
721
* are requested than would allow to place the most-significant digit
722
* a "0" is output. Eg:
723
* 1.4, 3 frac -> "1400" * 10^-3, ie 1.400
724
* 12.34 4 frac -> "123400" * 10^-4, ie 12.3400
725
* 3e-99 4 frac -> "0" * 10^0, ie 0
726
* 0.009 2 frac -> "0" * 10^-2, ie 0
728
* @param ieee_val Binary double description to convert. Must be the product
729
* of extract_ieee_double and it must not be a special number.
730
* @param signif_d_cnt Maximum number of significant digits to produce.
731
* The output is not rounded.
732
* Set to a negative value to generate as many digits
733
* as accurately possible.
734
* @param frac_d_cnt Maximum number of fractional digits to produce including
735
* any zeros immediately trailing the decimal point.
736
* The output is not rounded.
737
* Set to a negative value to generate as many digits
738
* as accurately possible.
739
* @param buf Buffer to store the string representation. Must be large
740
* enough to store all digits and a null terminator. At most
741
* MAX_DOUBLE_STR_LEN digits will be written (not counting
742
* the null terminator).
743
* @param buf_size Size of buf in bytes.
744
* @param dec_exponent Set to the decimal exponent of the number string
747
* @return The number of output digits. A negative value indicates
748
* an error: buf too small (or ieee_val.is_special, or
749
* signif_d_cnt == 0).
751
int double_to_fixed_str(ieee_double_t ieee_val, int signif_d_cnt,
752
int frac_d_cnt, char *buf, size_t buf_size, int *dec_exponent)
754
/* The whole computation assumes 64bit significand. */
755
assert(sizeof(ieee_val.pos_val.significand) == sizeof(uint64_t));
757
if (ieee_val.is_special) {
761
/* Zero cannot be normalized. Handle it here. */
762
if (0 == ieee_val.pos_val.significand) {
763
return zero_to_str(buf, buf_size, dec_exponent);
766
/* Normalize and scale. */
767
fp_num_t w = normalize(ieee_val.pos_val);
769
int lower_bin_exp = alpha - w.exponent - significand_width;
772
fp_num_t scaling_power_of_10;
774
get_power_of_ten(lower_bin_exp, &scaling_power_of_10, &scale);
776
fp_num_t w_scaled = multiply(w, scaling_power_of_10);
778
/* Produce decimal digits from the scaled number. */
779
int len = gen_fixed_dec_digits(w_scaled, scale, signif_d_cnt, frac_d_cnt,
780
buf, buf_size, dec_exponent);
782
assert(len <= MAX_DOUBLE_STR_LEN);
added
removed
uspace/lib/c/generic/double_to_str.c
