diff options
Diffstat (limited to 'source/fitz/ftoa.c')
| -rw-r--r-- | source/fitz/ftoa.c | 323 |
1 files changed, 305 insertions, 18 deletions
diff --git a/source/fitz/ftoa.c b/source/fitz/ftoa.c index 92ee9edf..eb6f750b 100644 --- a/source/fitz/ftoa.c +++ b/source/fitz/ftoa.c @@ -1,42 +1,329 @@ #include "mupdf/fitz.h" /* + Convert IEEE single precison numbers into decimal ASCII strings, while + satisfying the following two properties: + 1) Calling strtof or '(float) strtod' on the result must produce the + original float, independent of the rounding mode used by strtof/strtod. + 2) Minimize the number of produced decimal digits. E.g. the float 0.7f + should convert to "0.7", not "0.69999999". + + To solve this we use a dedicated single precision version of + Florian Loitsch's Grisu2 algorithm. See + http://florian.loitsch.com/publications/dtoa-pldi2010.pdf?attredirects=0 + + The code below is derived from Loitsch's C code, which + implements the same algorithm for IEEE double precision. See + http://florian.loitsch.com/publications/bench.tar.gz?attredirects=0 +*/ + +static int grisu2(float v, char* buffer, int* K); + +/* * compute decimal integer m, exp such that: * f = m*10^exp - * m is as short as possible with losing exactness + * m is as short as possible without losing exactness * assumes special cases (NaN, +Inf, -Inf) have been handled. */ void fz_ftoa(float f, char *s, int *exp, int *neg, int *ns) { - char buf[40], *p = buf; - int i; + /* Handle zero. */ + if (f == 0) + { + /* f is either +0 or -0. */ + *neg = 0; + *exp = 0; + *ns = 1; + s[0] = '0'; + s[1] = 0; + return; + } + + if (f > 0) + *neg = 0; + else + *neg = 1; + + *ns = grisu2(f, s, exp); + return; +} + +/* + Copyright (c) 2009 Florian Loitsch + + Permission is hereby granted, free of charge, to any person + obtaining a copy of this software and associated documentation + files (the "Software"), to deal in the Software without + restriction, including without limitation the rights to use, + copy, modify, merge, publish, distribute, sublicense, and/or sell + copies of the Software, and to permit persons to whom the + Software is furnished to do so, subject to the following + conditions: + + The above copyright notice and this permission notice shall be + included in all copies or substantial portions of the Software. + + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, + EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES + OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND + NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT + HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, + WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING + FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR + OTHER DEALINGS IN THE SOFTWARE. +*/ + +typedef union +{ + float d; + uint32_t n; +} converter_t; + +static uint32_t +float_to_uint32(float d) +{ + converter_t tmp; + tmp.d = d; + return tmp.n; +} + +typedef struct +{ + uint64_t f; + int e; +} diy_fp_t; + +#define DIY_SIGNIFICAND_SIZE 64 +#define DIY_LEADING_BIT ((uint64_t) 1 << (DIY_SIGNIFICAND_SIZE - 1)) - for (i = 0; i < 10; ++i) +static diy_fp_t +minus(diy_fp_t x, diy_fp_t y) +{ + diy_fp_t result = {x.f - y.f, x.e}; + assert(x.e == y.e && x.f >= y.f); + return result; +} + +static diy_fp_t +multiply(diy_fp_t x, diy_fp_t y) +{ + uint64_t a, b, c, d, ac, bc, ad, bd, tmp; + int half = DIY_SIGNIFICAND_SIZE / 2; + diy_fp_t r; uint64_t mask = ((uint64_t) 1 << half) - 1; + a = x.f >> half; b = x.f & mask; + c = y.f >> half; d = y.f & mask; + ac = a * c; bc = b * c; ad = a * d; bd = b * d; + tmp = (bd >> half) + (ad & mask) + (bc & mask); + tmp += 1U << (half - 1); /* Round. */ + r.f = ac + (ad >> half) + (bc >> half) + (tmp >> half); + r.e = x.e + y.e + half * 2; + return r; +} + +#define SP_SIGNIFICAND_SIZE 23 +#define SP_EXPONENT_BIAS (127 + SP_SIGNIFICAND_SIZE) +#define SP_MIN_EXPONENT (-SP_EXPONENT_BIAS) +#define SP_EXPONENT_MASK 0x7f800000 +#define SP_SIGNIFICAND_MASK 0x7fffff +#define SP_HIDDEN_BIT 0x800000 /* 2^23 */ + +/* Does not normalize the result. */ +static diy_fp_t +float2diy_fp(float d) +{ + uint32_t d32 = float_to_uint32(d); + int biased_e = (d32 & SP_EXPONENT_MASK) >> SP_SIGNIFICAND_SIZE; + uint32_t significand = d32 & SP_SIGNIFICAND_MASK; + diy_fp_t res; + + if (biased_e != 0) { - sprintf(buf, "%.*e", i, f); - if (fz_atof(buf) == f) - break; + res.f = significand + SP_HIDDEN_BIT; + res.e = biased_e - SP_EXPONENT_BIAS; } + else + { + res.f = significand; + res.e = SP_MIN_EXPONENT + 1; + } + return res; +} - if (*p == '-') +static diy_fp_t +normalize_boundary(diy_fp_t in) +{ + diy_fp_t res = in; + /* The original number could have been a denormal. */ + while (! (res.f & (SP_HIDDEN_BIT << 1))) { - *neg = 1; - ++p; + res.f <<= 1; + res.e--; + } + /* Do the final shifts in one go. */ + res.f <<= (DIY_SIGNIFICAND_SIZE - SP_SIGNIFICAND_SIZE - 2); + res.e = res.e - (DIY_SIGNIFICAND_SIZE - SP_SIGNIFICAND_SIZE - 2); + return res; +} + +static void +normalized_boundaries(float f, diy_fp_t* lower_ptr, diy_fp_t* upper_ptr) +{ + diy_fp_t v = float2diy_fp(f); + diy_fp_t upper, lower; + int significand_is_zero = v.f == SP_HIDDEN_BIT; + + upper.f = (v.f << 1) + 1; upper.e = v.e - 1; + upper = normalize_boundary(upper); + if (significand_is_zero) + { + lower.f = (v.f << 2) - 1; + lower.e = v.e - 2; } else - *neg = 0; + { + lower.f = (v.f << 1) - 1; + lower.e = v.e - 1; + } + lower.f <<= lower.e - upper.e; + lower.e = upper.e; + + /* Adjust to double boundaries, so that we can also read the numbers + with '(float) strtod'. */ + upper.f -= 1 << 10; + lower.f += 1 << 10; + + *upper_ptr = upper; + *lower_ptr = lower; +} + +static int +k_comp(int n) +{ + /* Avoid ceil and floating point multiplication for better + * performace and portability. Instead use the approximation + * log10(2) ~ 1233/(2^12). Tests show that this gives the correct + * result for all values of n in the range -500..500. */ + int tmp = n + DIY_SIGNIFICAND_SIZE - 1; + int k = (tmp * 1233) / (1 << 12); + return tmp > 0 ? k + 1 : k; +} - *ns = 0; - while (*p && *p != 'e') +/* Cached powers of ten from 10**-37..10**46. + Produced using GNU MPFR's mpfr_pow_si. */ + +/* Significands. */ +static uint64_t powers_ten[84] = { + 0x881cea14545c7575, 0xaa242499697392d3, 0xd4ad2dbfc3d07788, + 0x84ec3c97da624ab5, 0xa6274bbdd0fadd62, 0xcfb11ead453994ba, + 0x81ceb32c4b43fcf5, 0xa2425ff75e14fc32, 0xcad2f7f5359a3b3e, + 0xfd87b5f28300ca0e, 0x9e74d1b791e07e48, 0xc612062576589ddb, + 0xf79687aed3eec551, 0x9abe14cd44753b53, 0xc16d9a0095928a27, + 0xf1c90080baf72cb1, 0x971da05074da7bef, 0xbce5086492111aeb, + 0xec1e4a7db69561a5, 0x9392ee8e921d5d07, 0xb877aa3236a4b449, + 0xe69594bec44de15b, 0x901d7cf73ab0acd9, 0xb424dc35095cd80f, + 0xe12e13424bb40e13, 0x8cbccc096f5088cc, 0xafebff0bcb24aaff, + 0xdbe6fecebdedd5bf, 0x89705f4136b4a597, 0xabcc77118461cefd, + 0xd6bf94d5e57a42bc, 0x8637bd05af6c69b6, 0xa7c5ac471b478423, + 0xd1b71758e219652c, 0x83126e978d4fdf3b, 0xa3d70a3d70a3d70a, + 0xcccccccccccccccd, 0x8000000000000000, 0xa000000000000000, + 0xc800000000000000, 0xfa00000000000000, 0x9c40000000000000, + 0xc350000000000000, 0xf424000000000000, 0x9896800000000000, + 0xbebc200000000000, 0xee6b280000000000, 0x9502f90000000000, + 0xba43b74000000000, 0xe8d4a51000000000, 0x9184e72a00000000, + 0xb5e620f480000000, 0xe35fa931a0000000, 0x8e1bc9bf04000000, + 0xb1a2bc2ec5000000, 0xde0b6b3a76400000, 0x8ac7230489e80000, + 0xad78ebc5ac620000, 0xd8d726b7177a8000, 0x878678326eac9000, + 0xa968163f0a57b400, 0xd3c21bcecceda100, 0x84595161401484a0, + 0xa56fa5b99019a5c8, 0xcecb8f27f4200f3a, 0x813f3978f8940984, + 0xa18f07d736b90be5, 0xc9f2c9cd04674edf, 0xfc6f7c4045812296, + 0x9dc5ada82b70b59e, 0xc5371912364ce305, 0xf684df56c3e01bc7, + 0x9a130b963a6c115c, 0xc097ce7bc90715b3, 0xf0bdc21abb48db20, + 0x96769950b50d88f4, 0xbc143fa4e250eb31, 0xeb194f8e1ae525fd, + 0x92efd1b8d0cf37be, 0xb7abc627050305ae, 0xe596b7b0c643c719, + 0x8f7e32ce7bea5c70, 0xb35dbf821ae4f38c, 0xe0352f62a19e306f, +}; + +/* Exponents. */ +static uint64_t powers_ten_e[84] = { + -186, -183, -180, -176, -173, -170, -166, -163, -160, -157, -153, + -150, -147, -143, -140, -137, -133, -130, -127, -123, -120, -117, + -113, -110, -107, -103, -100, -97, -93, -90, -87, -83, -80, + -77, -73, -70, -67, -63, -60, -57, -54, -50, -47, -44, + -40, -37, -34, -30, -27, -24, -20, -17, -14, -10, -7, + -4, 0, 3, 6, 10, 13, 16, 20, 23, 26, 30, + 33, 36, 39, 43, 46, 49, 53, 56, 59, 63, 66, + 69, 73, 76, 79, 83, 86, 89 +}; + +static diy_fp_t +cached_power(int i) +{ + diy_fp_t result; + + assert (i >= -37 && i <= 46); + result.f = powers_ten[i + 37]; + result.e = powers_ten_e[i + 37]; + return result; +} + +/* Returns buffer length. */ +static int +digit_gen_mix_grisu2(diy_fp_t D_upper, diy_fp_t delta, char* buffer, int* K) +{ + int kappa; + diy_fp_t one = {(uint64_t) 1 << -D_upper.e, D_upper.e}; + unsigned char p1 = D_upper.f >> -one.e; + uint64_t p2 = D_upper.f & (one.f - 1); + unsigned char div = 10; + uint64_t mask = one.f - 1; + int len = 0; + for (kappa = 2; kappa > 0; --kappa) { - if (*p >= '0' && *p <= '9') + unsigned char digit = p1 / div; + if (digit || len) + buffer[len++] = '0' + digit; + p1 %= div; div /= 10; + if ((((uint64_t) p1) << -one.e) + p2 <= delta.f) { - *ns += 1; - *s++ = *p; + *K += kappa - 1; + return len; } - ++p; } + do + { + p2 *= 10; + buffer[len++] = '0' + (p2 >> -one.e); + p2 &= mask; + kappa--; + delta.f *= 10; + } + while (p2 > delta.f); + *K += kappa; + return len; +} + +static int +grisu2(float v, char* buffer, int* K) +{ + diy_fp_t w_lower, w_upper, D_upper, D_lower, c_mk, delta; + int length, mk, alpha = -DIY_SIGNIFICAND_SIZE + 4; + + normalized_boundaries(v, &w_lower, &w_upper); + mk = k_comp(alpha - w_upper.e - DIY_SIGNIFICAND_SIZE); + c_mk = cached_power(mk); + + D_upper = multiply(w_upper, c_mk); + D_lower = multiply(w_lower, c_mk); + + D_upper.f--; + D_lower.f++; + + delta = minus(D_upper, D_lower); + + *K = -mk; + length = digit_gen_mix_grisu2(D_upper, delta, buffer, K); - *exp = fz_atoi(p+1) - (*ns) + 1; + buffer[length] = 0; + return length; } |