1#include "mupdf/fitz.h"
2
3#include <assert.h>
4
5/*
6 Convert IEEE single precision numbers into decimal ASCII strings, while
7 satisfying the following two properties:
8 1) Calling strtof or '(float) strtod' on the result must produce the
9 original float, independent of the rounding mode used by strtof/strtod.
10 2) Minimize the number of produced decimal digits. E.g. the float 0.7f
11 should convert to "0.7", not "0.69999999".
12
13 To solve this we use a dedicated single precision version of
14 Florian Loitsch's Grisu2 algorithm. See
15 http://florian.loitsch.com/publications/dtoa-pldi2010.pdf?attredirects=0
16
17 The code below is derived from Loitsch's C code, which
18 implements the same algorithm for IEEE double precision. See
19 http://florian.loitsch.com/publications/bench.tar.gz?attredirects=0
20*/
21
22/*
23 Copyright (c) 2009 Florian Loitsch
24
25 Permission is hereby granted, free of charge, to any person
26 obtaining a copy of this software and associated documentation
27 files (the "Software"), to deal in the Software without
28 restriction, including without limitation the rights to use,
29 copy, modify, merge, publish, distribute, sublicense, and/or sell
30 copies of the Software, and to permit persons to whom the
31 Software is furnished to do so, subject to the following
32 conditions:
33
34 The above copyright notice and this permission notice shall be
35 included in all copies or substantial portions of the Software.
36
37 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
38 EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
39 OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
40 NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
41 HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
42 WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
43 FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
44 OTHER DEALINGS IN THE SOFTWARE.
45*/
46
47static uint32_t
48float_to_uint32(float d)
49{
50 union
51 {
52 float d;
53 uint32_t n;
54 } tmp;
55 tmp.d = d;
56 return tmp.n;
57}
58
59typedef struct
60{
61 uint64_t f;
62 int e;
63} diy_fp_t;
64
65#define DIY_SIGNIFICAND_SIZE 64
66#define DIY_LEADING_BIT ((uint64_t) 1 << (DIY_SIGNIFICAND_SIZE - 1))
67
68static diy_fp_t
69minus(diy_fp_t x, diy_fp_t y)
70{
71 diy_fp_t result = {x.f - y.f, x.e};
72 assert(x.e == y.e && x.f >= y.f);
73 return result;
74}
75
76static diy_fp_t
77multiply(diy_fp_t x, diy_fp_t y)
78{
79 uint64_t a, b, c, d, ac, bc, ad, bd, tmp;
80 int half = DIY_SIGNIFICAND_SIZE / 2;
81 diy_fp_t r; uint64_t mask = ((uint64_t) 1 << half) - 1;
82 a = x.f >> half; b = x.f & mask;
83 c = y.f >> half; d = y.f & mask;
84 ac = a * c; bc = b * c; ad = a * d; bd = b * d;
85 tmp = (bd >> half) + (ad & mask) + (bc & mask);
86 tmp += ((uint64_t)1U) << (half - 1); /* Round. */
87 r.f = ac + (ad >> half) + (bc >> half) + (tmp >> half);
88 r.e = x.e + y.e + half * 2;
89 return r;
90}
91
92#define SP_SIGNIFICAND_SIZE 23
93#define SP_EXPONENT_BIAS (127 + SP_SIGNIFICAND_SIZE)
94#define SP_MIN_EXPONENT (-SP_EXPONENT_BIAS)
95#define SP_EXPONENT_MASK 0x7f800000
96#define SP_SIGNIFICAND_MASK 0x7fffff
97#define SP_HIDDEN_BIT 0x800000 /* 2^23 */
98
99/* Does not normalize the result. */
100static diy_fp_t
101float2diy_fp(float d)
102{
103 uint32_t d32 = float_to_uint32(d);
104 int biased_e = (d32 & SP_EXPONENT_MASK) >> SP_SIGNIFICAND_SIZE;
105 uint32_t significand = d32 & SP_SIGNIFICAND_MASK;
106 diy_fp_t res;
107
108 if (biased_e != 0)
109 {
110 res.f = significand + SP_HIDDEN_BIT;
111 res.e = biased_e - SP_EXPONENT_BIAS;
112 }
113 else
114 {
115 res.f = significand;
116 res.e = SP_MIN_EXPONENT + 1;
117 }
118 return res;
119}
120
121static diy_fp_t
122normalize_boundary(diy_fp_t in)
123{
124 diy_fp_t res = in;
125 /* The original number could have been a denormal. */
126 while (! (res.f & (SP_HIDDEN_BIT << 1)))
127 {
128 res.f <<= 1;
129 res.e--;
130 }
131 /* Do the final shifts in one go. */
132 res.f <<= (DIY_SIGNIFICAND_SIZE - SP_SIGNIFICAND_SIZE - 2);
133 res.e = res.e - (DIY_SIGNIFICAND_SIZE - SP_SIGNIFICAND_SIZE - 2);
134 return res;
135}
136
137static void
138normalized_boundaries(float f, diy_fp_t* lower_ptr, diy_fp_t* upper_ptr)
139{
140 diy_fp_t v = float2diy_fp(f);
141 diy_fp_t upper, lower;
142 int significand_is_zero = v.f == SP_HIDDEN_BIT;
143
144 upper.f = (v.f << 1) + 1; upper.e = v.e - 1;
145 upper = normalize_boundary(upper);
146 if (significand_is_zero)
147 {
148 lower.f = (v.f << 2) - 1;
149 lower.e = v.e - 2;
150 }
151 else
152 {
153 lower.f = (v.f << 1) - 1;
154 lower.e = v.e - 1;
155 }
156 lower.f <<= lower.e - upper.e;
157 lower.e = upper.e;
158
159 /* Adjust to double boundaries, so that we can also read the numbers with '(float) strtod'. */
160 upper.f -= 1 << 10;
161 lower.f += 1 << 10;
162
163 *upper_ptr = upper;
164 *lower_ptr = lower;
165}
166
167static int
168k_comp(int n)
169{
170 /* Avoid ceil and floating point multiplication for better
171 * performance and portability. Instead use the approximation
172 * log10(2) ~ 1233/(2^12). Tests show that this gives the correct
173 * result for all values of n in the range -500..500. */
174 int tmp = n + DIY_SIGNIFICAND_SIZE - 1;
175 int k = (tmp * 1233) / (1 << 12);
176 return tmp > 0 ? k + 1 : k;
177}
178
179/* Cached powers of ten from 10**-37..10**46. Produced using GNU MPFR's mpfr_pow_si. */
180
181/* Significands. */
182static uint64_t powers_ten[84] = {
183 0x881cea14545c7575ull, 0xaa242499697392d3ull, 0xd4ad2dbfc3d07788ull,
184 0x84ec3c97da624ab5ull, 0xa6274bbdd0fadd62ull, 0xcfb11ead453994baull,
185 0x81ceb32c4b43fcf5ull, 0xa2425ff75e14fc32ull, 0xcad2f7f5359a3b3eull,
186 0xfd87b5f28300ca0eull, 0x9e74d1b791e07e48ull, 0xc612062576589ddbull,
187 0xf79687aed3eec551ull, 0x9abe14cd44753b53ull, 0xc16d9a0095928a27ull,
188 0xf1c90080baf72cb1ull, 0x971da05074da7befull, 0xbce5086492111aebull,
189 0xec1e4a7db69561a5ull, 0x9392ee8e921d5d07ull, 0xb877aa3236a4b449ull,
190 0xe69594bec44de15bull, 0x901d7cf73ab0acd9ull, 0xb424dc35095cd80full,
191 0xe12e13424bb40e13ull, 0x8cbccc096f5088ccull, 0xafebff0bcb24aaffull,
192 0xdbe6fecebdedd5bfull, 0x89705f4136b4a597ull, 0xabcc77118461cefdull,
193 0xd6bf94d5e57a42bcull, 0x8637bd05af6c69b6ull, 0xa7c5ac471b478423ull,
194 0xd1b71758e219652cull, 0x83126e978d4fdf3bull, 0xa3d70a3d70a3d70aull,
195 0xcccccccccccccccdull, 0x8000000000000000ull, 0xa000000000000000ull,
196 0xc800000000000000ull, 0xfa00000000000000ull, 0x9c40000000000000ull,
197 0xc350000000000000ull, 0xf424000000000000ull, 0x9896800000000000ull,
198 0xbebc200000000000ull, 0xee6b280000000000ull, 0x9502f90000000000ull,
199 0xba43b74000000000ull, 0xe8d4a51000000000ull, 0x9184e72a00000000ull,
200 0xb5e620f480000000ull, 0xe35fa931a0000000ull, 0x8e1bc9bf04000000ull,
201 0xb1a2bc2ec5000000ull, 0xde0b6b3a76400000ull, 0x8ac7230489e80000ull,
202 0xad78ebc5ac620000ull, 0xd8d726b7177a8000ull, 0x878678326eac9000ull,
203 0xa968163f0a57b400ull, 0xd3c21bcecceda100ull, 0x84595161401484a0ull,
204 0xa56fa5b99019a5c8ull, 0xcecb8f27f4200f3aull, 0x813f3978f8940984ull,
205 0xa18f07d736b90be5ull, 0xc9f2c9cd04674edfull, 0xfc6f7c4045812296ull,
206 0x9dc5ada82b70b59eull, 0xc5371912364ce305ull, 0xf684df56c3e01bc7ull,
207 0x9a130b963a6c115cull, 0xc097ce7bc90715b3ull, 0xf0bdc21abb48db20ull,
208 0x96769950b50d88f4ull, 0xbc143fa4e250eb31ull, 0xeb194f8e1ae525fdull,
209 0x92efd1b8d0cf37beull, 0xb7abc627050305aeull, 0xe596b7b0c643c719ull,
210 0x8f7e32ce7bea5c70ull, 0xb35dbf821ae4f38cull, 0xe0352f62a19e306full,
211};
212
213/* Exponents. */
214static int powers_ten_e[84] = {
215 -186, -183, -180, -176, -173, -170, -166, -163, -160, -157, -153,
216 -150, -147, -143, -140, -137, -133, -130, -127, -123, -120, -117,
217 -113, -110, -107, -103, -100, -97, -93, -90, -87, -83, -80,
218 -77, -73, -70, -67, -63, -60, -57, -54, -50, -47, -44,
219 -40, -37, -34, -30, -27, -24, -20, -17, -14, -10, -7,
220 -4, 0, 3, 6, 10, 13, 16, 20, 23, 26, 30,
221 33, 36, 39, 43, 46, 49, 53, 56, 59, 63, 66,
222 69, 73, 76, 79, 83, 86, 89
223};
224
225static diy_fp_t
226cached_power(int i)
227{
228 diy_fp_t result;
229
230 assert (i >= -37 && i <= 46);
231 result.f = powers_ten[i + 37];
232 result.e = powers_ten_e[i + 37];
233 return result;
234}
235
236/* Returns buffer length. */
237static int
238digit_gen_mix_grisu2(diy_fp_t D_upper, diy_fp_t delta, char* buffer, int* K)
239{
240 int kappa;
241 diy_fp_t one = {(uint64_t) 1 << -D_upper.e, D_upper.e};
242 unsigned char p1 = D_upper.f >> -one.e;
243 uint64_t p2 = D_upper.f & (one.f - 1);
244 unsigned char div = 10;
245 uint64_t mask = one.f - 1;
246 int len = 0;
247 for (kappa = 2; kappa > 0; --kappa)
248 {
249 unsigned char digit = p1 / div;
250 if (digit || len)
251 buffer[len++] = '0' + digit;
252 p1 %= div; div /= 10;
253 if ((((uint64_t) p1) << -one.e) + p2 <= delta.f)
254 {
255 *K += kappa - 1;
256 return len;
257 }
258 }
259 do
260 {
261 p2 *= 10;
262 buffer[len++] = '0' + (p2 >> -one.e);
263 p2 &= mask;
264 kappa--;
265 delta.f *= 10;
266 }
267 while (p2 > delta.f);
268 *K += kappa;
269 return len;
270}
271
272/*
273 Compute decimal integer m, exp such that:
274 f = m * 10^exp
275 m is as short as possible without losing exactness
276 Assumes special cases (0, NaN, +Inf, -Inf) have been handled.
277*/
278int
279fz_grisu(float v, char* buffer, int* K)
280{
281 diy_fp_t w_lower, w_upper, D_upper, D_lower, c_mk, delta;
282 int length, mk, alpha = -DIY_SIGNIFICAND_SIZE + 4;
283
284 normalized_boundaries(v, &w_lower, &w_upper);
285 mk = k_comp(alpha - w_upper.e - DIY_SIGNIFICAND_SIZE);
286 c_mk = cached_power(mk);
287
288 D_upper = multiply(w_upper, c_mk);
289 D_lower = multiply(w_lower, c_mk);
290
291 D_upper.f--;
292 D_lower.f++;
293
294 delta = minus(D_upper, D_lower);
295
296 *K = -mk;
297 length = digit_gen_mix_grisu2(D_upper, delta, buffer, K);
298
299 buffer[length] = 0;
300 return length;
301}
302