Coverage Report

Created: 2025-01-28 06:17

/src/mupdf/source/fitz/ftoa.c
Line
Count
Source (jump to first uncovered line)
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
47
static uint32_t
48
float_to_uint32(float d)
49
267k
{
50
267k
  union
51
267k
  {
52
267k
    float d;
53
267k
    uint32_t n;
54
267k
  } tmp;
55
267k
  tmp.d = d;
56
267k
  return tmp.n;
57
267k
}
58
59
typedef struct
60
{
61
  uint64_t f;
62
  int e;
63
} diy_fp_t;
64
65
1.87M
#define DIY_SIGNIFICAND_SIZE 64
66
#define DIY_LEADING_BIT ((uint64_t) 1 << (DIY_SIGNIFICAND_SIZE - 1))
67
68
static diy_fp_t
69
minus(diy_fp_t x, diy_fp_t y)
70
267k
{
71
267k
  diy_fp_t result = {x.f - y.f, x.e};
72
267k
  assert(x.e == y.e && x.f >= y.f);
73
267k
  return result;
74
267k
}
75
76
static diy_fp_t
77
multiply(diy_fp_t x, diy_fp_t y)
78
535k
{
79
535k
  uint64_t a, b, c, d, ac, bc, ad, bd, tmp;
80
535k
  int half = DIY_SIGNIFICAND_SIZE / 2;
81
535k
  diy_fp_t r; uint64_t mask = ((uint64_t) 1 << half) - 1;
82
535k
  a = x.f >> half; b = x.f & mask;
83
535k
  c = y.f >> half; d = y.f & mask;
84
535k
  ac = a * c; bc = b * c; ad = a * d; bd = b * d;
85
535k
  tmp = (bd >> half) + (ad & mask) + (bc & mask);
86
535k
  tmp += ((uint64_t)1U) << (half - 1); /* Round. */
87
535k
  r.f = ac + (ad >> half) + (bc >> half) + (tmp >> half);
88
535k
  r.e = x.e + y.e + half * 2;
89
535k
  return r;
90
535k
}
91
92
1.07M
#define SP_SIGNIFICAND_SIZE 23
93
267k
#define SP_EXPONENT_BIAS (127 + SP_SIGNIFICAND_SIZE)
94
0
#define SP_MIN_EXPONENT (-SP_EXPONENT_BIAS)
95
267k
#define SP_EXPONENT_MASK 0x7f800000
96
267k
#define SP_SIGNIFICAND_MASK 0x7fffff
97
802k
#define SP_HIDDEN_BIT 0x800000 /* 2^23 */
98
99
/* Does not normalize the result. */
100
static diy_fp_t
101
float2diy_fp(float d)
102
267k
{
103
267k
  uint32_t d32 = float_to_uint32(d);
104
267k
  int biased_e = (d32 & SP_EXPONENT_MASK) >> SP_SIGNIFICAND_SIZE;
105
267k
  uint32_t significand = d32 & SP_SIGNIFICAND_MASK;
106
267k
  diy_fp_t res;
107
108
267k
  if (biased_e != 0)
109
267k
  {
110
267k
    res.f = significand + SP_HIDDEN_BIT;
111
267k
    res.e = biased_e - SP_EXPONENT_BIAS;
112
267k
  }
113
0
  else
114
0
  {
115
0
    res.f = significand;
116
0
    res.e = SP_MIN_EXPONENT + 1;
117
0
  }
118
267k
  return res;
119
267k
}
120
121
static diy_fp_t
122
normalize_boundary(diy_fp_t in)
123
267k
{
124
267k
  diy_fp_t res = in;
125
  /* The original number could have been a denormal. */
126
267k
  while (! (res.f & (SP_HIDDEN_BIT << 1)))
127
0
  {
128
0
    res.f <<= 1;
129
0
    res.e--;
130
0
  }
131
  /* Do the final shifts in one go. */
132
267k
  res.f <<= (DIY_SIGNIFICAND_SIZE - SP_SIGNIFICAND_SIZE - 2);
133
267k
  res.e = res.e - (DIY_SIGNIFICAND_SIZE - SP_SIGNIFICAND_SIZE - 2);
134
267k
  return res;
135
267k
}
136
137
static void
138
normalized_boundaries(float f, diy_fp_t* lower_ptr, diy_fp_t* upper_ptr)
139
267k
{
140
267k
  diy_fp_t v = float2diy_fp(f);
141
267k
  diy_fp_t upper, lower;
142
267k
  int significand_is_zero = v.f == SP_HIDDEN_BIT;
143
144
267k
  upper.f = (v.f << 1) + 1; upper.e = v.e - 1;
145
267k
  upper = normalize_boundary(upper);
146
267k
  if (significand_is_zero)
147
86.8k
  {
148
86.8k
    lower.f = (v.f << 2) - 1;
149
86.8k
    lower.e = v.e - 2;
150
86.8k
  }
151
180k
  else
152
180k
  {
153
180k
    lower.f = (v.f << 1) - 1;
154
180k
    lower.e = v.e - 1;
155
180k
  }
156
267k
  lower.f <<= lower.e - upper.e;
157
267k
  lower.e = upper.e;
158
159
  /* Adjust to double boundaries, so that we can also read the numbers with '(float) strtod'. */
160
267k
  upper.f -= 1 << 10;
161
267k
  lower.f += 1 << 10;
162
163
267k
  *upper_ptr = upper;
164
267k
  *lower_ptr = lower;
165
267k
}
166
167
static int
168
k_comp(int n)
169
267k
{
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
267k
  int tmp = n + DIY_SIGNIFICAND_SIZE - 1;
175
267k
  int k = (tmp * 1233) / (1 << 12);
176
267k
  return tmp > 0 ? k + 1 : k;
177
267k
}
178
179
/* Cached powers of ten from 10**-37..10**46. Produced using GNU MPFR's mpfr_pow_si. */
180
181
/* Significands. */
182
static 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. */
214
static 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
225
static diy_fp_t
226
cached_power(int i)
227
267k
{
228
267k
  diy_fp_t result;
229
230
267k
  assert (i >= -37 && i <= 46);
231
267k
  result.f = powers_ten[i + 37];
232
267k
  result.e = powers_ten_e[i + 37];
233
267k
  return result;
234
267k
}
235
236
/* Returns buffer length. */
237
static int
238
digit_gen_mix_grisu2(diy_fp_t D_upper, diy_fp_t delta, char* buffer, int* K)
239
267k
{
240
267k
  int kappa;
241
267k
  diy_fp_t one = {(uint64_t) 1 << -D_upper.e, D_upper.e};
242
267k
  unsigned char p1 = D_upper.f >> -one.e;
243
267k
  uint64_t p2 = D_upper.f & (one.f - 1);
244
267k
  unsigned char div = 10;
245
267k
  uint64_t mask = one.f - 1;
246
267k
  int len = 0;
247
580k
  for (kappa = 2; kappa > 0; --kappa)
248
444k
  {
249
444k
    unsigned char digit = p1 / div;
250
444k
    if (digit || len)
251
413k
      buffer[len++] = '0' + digit;
252
444k
    p1 %= div; div /= 10;
253
444k
    if ((((uint64_t) p1) << -one.e) + p2 <= delta.f)
254
131k
    {
255
131k
      *K += kappa - 1;
256
131k
      return len;
257
131k
    }
258
444k
  }
259
136k
  do
260
620k
  {
261
620k
    p2 *= 10;
262
620k
    buffer[len++] = '0' + (p2 >> -one.e);
263
620k
    p2 &= mask;
264
620k
    kappa--;
265
620k
    delta.f *= 10;
266
620k
  }
267
620k
  while (p2 > delta.f);
268
136k
  *K += kappa;
269
136k
  return len;
270
267k
}
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
*/
278
int
279
fz_grisu(float v, char* buffer, int* K)
280
267k
{
281
267k
  diy_fp_t w_lower, w_upper, D_upper, D_lower, c_mk, delta;
282
267k
  int length, mk, alpha = -DIY_SIGNIFICAND_SIZE + 4;
283
284
267k
  normalized_boundaries(v, &w_lower, &w_upper);
285
267k
  mk = k_comp(alpha - w_upper.e - DIY_SIGNIFICAND_SIZE);
286
267k
  c_mk = cached_power(mk);
287
288
267k
  D_upper = multiply(w_upper, c_mk);
289
267k
  D_lower = multiply(w_lower, c_mk);
290
291
267k
  D_upper.f--;
292
267k
  D_lower.f++;
293
294
267k
  delta = minus(D_upper, D_lower);
295
296
267k
  *K = -mk;
297
267k
  length = digit_gen_mix_grisu2(D_upper, delta, buffer, K);
298
299
267k
  buffer[length] = 0;
300
267k
  return length;
301
267k
}