]>
Commit | Line | Data |
---|---|---|
8faf50e0 XL |
1 | //! Almost direct (but slightly optimized) Rust translation of Figure 3 of "Printing |
2 | //! Floating-Point Numbers Quickly and Accurately"[^1]. | |
3 | //! | |
4 | //! [^1]: Burger, R. G. and Dybvig, R. K. 1996. Printing floating-point numbers | |
5 | //! quickly and accurately. SIGPLAN Not. 31, 5 (May. 1996), 108-116. | |
d9579d0f | 6 | |
48663c56 | 7 | use crate::cmp::Ordering; |
d9579d0f | 8 | |
48663c56 | 9 | use crate::num::bignum::Big32x40 as Big; |
60c5eb7d XL |
10 | use crate::num::bignum::Digit32 as Digit; |
11 | use crate::num::flt2dec::estimator::estimate_scaling_factor; | |
12 | use crate::num::flt2dec::{round_up, Decoded, MAX_SIG_DIGITS}; | |
d9579d0f | 13 | |
60c5eb7d XL |
14 | static POW10: [Digit; 10] = |
15 | [1, 10, 100, 1000, 10000, 100000, 1000000, 10000000, 100000000, 1000000000]; | |
16 | static TWOPOW10: [Digit; 10] = | |
17 | [2, 20, 200, 2000, 20000, 200000, 2000000, 20000000, 200000000, 2000000000]; | |
d9579d0f AL |
18 | |
19 | // precalculated arrays of `Digit`s for 10^(2^n) | |
20 | static POW10TO16: [Digit; 2] = [0x6fc10000, 0x2386f2]; | |
21 | static POW10TO32: [Digit; 4] = [0, 0x85acef81, 0x2d6d415b, 0x4ee]; | |
22 | static POW10TO64: [Digit; 7] = [0, 0, 0xbf6a1f01, 0x6e38ed64, 0xdaa797ed, 0xe93ff9f4, 0x184f03]; | |
60c5eb7d XL |
23 | static POW10TO128: [Digit; 14] = [ |
24 | 0, 0, 0, 0, 0x2e953e01, 0x3df9909, 0xf1538fd, 0x2374e42f, 0xd3cff5ec, 0xc404dc08, 0xbccdb0da, | |
25 | 0xa6337f19, 0xe91f2603, 0x24e, | |
26 | ]; | |
27 | static POW10TO256: [Digit; 27] = [ | |
28 | 0, 0, 0, 0, 0, 0, 0, 0, 0x982e7c01, 0xbed3875b, 0xd8d99f72, 0x12152f87, 0x6bde50c6, 0xcf4a6e70, | |
29 | 0xd595d80f, 0x26b2716e, 0xadc666b0, 0x1d153624, 0x3c42d35a, 0x63ff540e, 0xcc5573c0, 0x65f9ef17, | |
30 | 0x55bc28f2, 0x80dcc7f7, 0xf46eeddc, 0x5fdcefce, 0x553f7, | |
31 | ]; | |
d9579d0f AL |
32 | |
33 | #[doc(hidden)] | |
e9174d1e | 34 | pub fn mul_pow10(x: &mut Big, n: usize) -> &mut Big { |
d9579d0f | 35 | debug_assert!(n < 512); |
60c5eb7d XL |
36 | if n & 7 != 0 { |
37 | x.mul_small(POW10[n & 7]); | |
38 | } | |
39 | if n & 8 != 0 { | |
40 | x.mul_small(POW10[8]); | |
41 | } | |
42 | if n & 16 != 0 { | |
43 | x.mul_digits(&POW10TO16); | |
44 | } | |
45 | if n & 32 != 0 { | |
46 | x.mul_digits(&POW10TO32); | |
47 | } | |
48 | if n & 64 != 0 { | |
49 | x.mul_digits(&POW10TO64); | |
50 | } | |
51 | if n & 128 != 0 { | |
52 | x.mul_digits(&POW10TO128); | |
53 | } | |
54 | if n & 256 != 0 { | |
55 | x.mul_digits(&POW10TO256); | |
56 | } | |
d9579d0f AL |
57 | x |
58 | } | |
59 | ||
e9174d1e | 60 | fn div_2pow10(x: &mut Big, mut n: usize) -> &mut Big { |
d9579d0f AL |
61 | let largest = POW10.len() - 1; |
62 | while n > largest { | |
63 | x.div_rem_small(POW10[largest]); | |
64 | n -= largest; | |
65 | } | |
66 | x.div_rem_small(TWOPOW10[n]); | |
67 | x | |
68 | } | |
69 | ||
70 | // only usable when `x < 16 * scale`; `scaleN` should be `scale.mul_small(N)` | |
60c5eb7d XL |
71 | fn div_rem_upto_16<'a>( |
72 | x: &'a mut Big, | |
73 | scale: &Big, | |
74 | scale2: &Big, | |
75 | scale4: &Big, | |
76 | scale8: &Big, | |
77 | ) -> (u8, &'a mut Big) { | |
d9579d0f | 78 | let mut d = 0; |
60c5eb7d XL |
79 | if *x >= *scale8 { |
80 | x.sub(scale8); | |
81 | d += 8; | |
82 | } | |
83 | if *x >= *scale4 { | |
84 | x.sub(scale4); | |
85 | d += 4; | |
86 | } | |
87 | if *x >= *scale2 { | |
88 | x.sub(scale2); | |
89 | d += 2; | |
90 | } | |
91 | if *x >= *scale { | |
92 | x.sub(scale); | |
93 | d += 1; | |
94 | } | |
d9579d0f AL |
95 | debug_assert!(*x < *scale); |
96 | (d, x) | |
97 | } | |
98 | ||
99 | /// The shortest mode implementation for Dragon. | |
100 | pub fn format_shortest(d: &Decoded, buf: &mut [u8]) -> (/*#digits*/ usize, /*exp*/ i16) { | |
101 | // the number `v` to format is known to be: | |
102 | // - equal to `mant * 2^exp`; | |
103 | // - preceded by `(mant - 2 * minus) * 2^exp` in the original type; and | |
104 | // - followed by `(mant + 2 * plus) * 2^exp` in the original type. | |
105 | // | |
106 | // obviously, `minus` and `plus` cannot be zero. (for infinities, we use out-of-range values.) | |
0731742a | 107 | // also we assume that at least one digit is generated, i.e., `mant` cannot be zero too. |
d9579d0f AL |
108 | // |
109 | // this also means that any number between `low = (mant - minus) * 2^exp` and | |
110 | // `high = (mant + plus) * 2^exp` will map to this exact floating point number, | |
0731742a | 111 | // with bounds included when the original mantissa was even (i.e., `!mant_was_odd`). |
d9579d0f AL |
112 | |
113 | assert!(d.mant > 0); | |
114 | assert!(d.minus > 0); | |
115 | assert!(d.plus > 0); | |
116 | assert!(d.mant.checked_add(d.plus).is_some()); | |
117 | assert!(d.mant.checked_sub(d.minus).is_some()); | |
118 | assert!(buf.len() >= MAX_SIG_DIGITS); | |
119 | ||
120 | // `a.cmp(&b) < rounding` is `if d.inclusive {a <= b} else {a < b}` | |
60c5eb7d | 121 | let rounding = if d.inclusive { Ordering::Greater } else { Ordering::Equal }; |
d9579d0f AL |
122 | |
123 | // estimate `k_0` from original inputs satisfying `10^(k_0-1) < high <= 10^(k_0+1)`. | |
124 | // the tight bound `k` satisfying `10^(k-1) < high <= 10^k` is calculated later. | |
125 | let mut k = estimate_scaling_factor(d.mant + d.plus, d.exp); | |
126 | ||
127 | // convert `{mant, plus, minus} * 2^exp` into the fractional form so that: | |
128 | // - `v = mant / scale` | |
129 | // - `low = (mant - minus) / scale` | |
130 | // - `high = (mant + plus) / scale` | |
131 | let mut mant = Big::from_u64(d.mant); | |
132 | let mut minus = Big::from_u64(d.minus); | |
133 | let mut plus = Big::from_u64(d.plus); | |
134 | let mut scale = Big::from_small(1); | |
135 | if d.exp < 0 { | |
136 | scale.mul_pow2(-d.exp as usize); | |
137 | } else { | |
138 | mant.mul_pow2(d.exp as usize); | |
139 | minus.mul_pow2(d.exp as usize); | |
140 | plus.mul_pow2(d.exp as usize); | |
141 | } | |
142 | ||
143 | // divide `mant` by `10^k`. now `scale / 10 < mant + plus <= scale * 10`. | |
144 | if k >= 0 { | |
145 | mul_pow10(&mut scale, k as usize); | |
146 | } else { | |
147 | mul_pow10(&mut mant, -k as usize); | |
148 | mul_pow10(&mut minus, -k as usize); | |
149 | mul_pow10(&mut plus, -k as usize); | |
150 | } | |
151 | ||
152 | // fixup when `mant + plus > scale` (or `>=`). | |
153 | // we are not actually modifying `scale`, since we can skip the initial multiplication instead. | |
154 | // now `scale < mant + plus <= scale * 10` and we are ready to generate digits. | |
155 | // | |
156 | // note that `d[0]` *can* be zero, when `scale - plus < mant < scale`. | |
157 | // in this case rounding-up condition (`up` below) will be triggered immediately. | |
158 | if scale.cmp(mant.clone().add(&plus)) < rounding { | |
159 | // equivalent to scaling `scale` by 10 | |
160 | k += 1; | |
161 | } else { | |
162 | mant.mul_small(10); | |
163 | minus.mul_small(10); | |
164 | plus.mul_small(10); | |
165 | } | |
166 | ||
167 | // cache `(2, 4, 8) * scale` for digit generation. | |
60c5eb7d XL |
168 | let mut scale2 = scale.clone(); |
169 | scale2.mul_pow2(1); | |
170 | let mut scale4 = scale.clone(); | |
171 | scale4.mul_pow2(2); | |
172 | let mut scale8 = scale.clone(); | |
173 | scale8.mul_pow2(3); | |
d9579d0f AL |
174 | |
175 | let mut down; | |
176 | let mut up; | |
177 | let mut i = 0; | |
178 | loop { | |
179 | // invariants, where `d[0..n-1]` are digits generated so far: | |
180 | // - `v = mant / scale * 10^(k-n-1) + d[0..n-1] * 10^(k-n)` | |
181 | // - `v - low = minus / scale * 10^(k-n-1)` | |
182 | // - `high - v = plus / scale * 10^(k-n-1)` | |
183 | // - `(mant + plus) / scale <= 10` (thus `mant / scale < 10`) | |
184 | // where `d[i..j]` is a shorthand for `d[i] * 10^(j-i) + ... + d[j-1] * 10 + d[j]`. | |
185 | ||
186 | // generate one digit: `d[n] = floor(mant / scale) < 10`. | |
187 | let (d, _) = div_rem_upto_16(&mut mant, &scale, &scale2, &scale4, &scale8); | |
188 | debug_assert!(d < 10); | |
189 | buf[i] = b'0' + d; | |
190 | i += 1; | |
191 | ||
192 | // this is a simplified description of the modified Dragon algorithm. | |
193 | // many intermediate derivations and completeness arguments are omitted for convenience. | |
194 | // | |
195 | // start with modified invariants, as we've updated `n`: | |
196 | // - `v = mant / scale * 10^(k-n) + d[0..n-1] * 10^(k-n)` | |
197 | // - `v - low = minus / scale * 10^(k-n)` | |
198 | // - `high - v = plus / scale * 10^(k-n)` | |
199 | // | |
200 | // assume that `d[0..n-1]` is the shortest representation between `low` and `high`, | |
0731742a | 201 | // i.e., `d[0..n-1]` satisfies both of the following but `d[0..n-2]` doesn't: |
d9579d0f AL |
202 | // - `low < d[0..n-1] * 10^(k-n) < high` (bijectivity: digits round to `v`); and |
203 | // - `abs(v / 10^(k-n) - d[0..n-1]) <= 1/2` (the last digit is correct). | |
204 | // | |
205 | // the second condition simplifies to `2 * mant <= scale`. | |
206 | // solving invariants in terms of `mant`, `low` and `high` yields | |
207 | // a simpler version of the first condition: `-plus < mant < minus`. | |
208 | // since `-plus < 0 <= mant`, we have the correct shortest representation | |
209 | // when `mant < minus` and `2 * mant <= scale`. | |
210 | // (the former becomes `mant <= minus` when the original mantissa is even.) | |
211 | // | |
212 | // when the second doesn't hold (`2 * mant > scale`), we need to increase the last digit. | |
213 | // this is enough for restoring that condition: we already know that | |
214 | // the digit generation guarantees `0 <= v / 10^(k-n) - d[0..n-1] < 1`. | |
215 | // in this case, the first condition becomes `-plus < mant - scale < minus`. | |
216 | // since `mant < scale` after the generation, we have `scale < mant + plus`. | |
217 | // (again, this becomes `scale <= mant + plus` when the original mantissa is even.) | |
218 | // | |
219 | // in short: | |
220 | // - stop and round `down` (keep digits as is) when `mant < minus` (or `<=`). | |
221 | // - stop and round `up` (increase the last digit) when `scale < mant + plus` (or `<=`). | |
222 | // - keep generating otherwise. | |
223 | down = mant.cmp(&minus) < rounding; | |
224 | up = scale.cmp(mant.clone().add(&plus)) < rounding; | |
60c5eb7d XL |
225 | if down || up { |
226 | break; | |
227 | } // we have the shortest representation, proceed to the rounding | |
d9579d0f AL |
228 | |
229 | // restore the invariants. | |
230 | // this makes the algorithm always terminating: `minus` and `plus` always increases, | |
231 | // but `mant` is clipped modulo `scale` and `scale` is fixed. | |
232 | mant.mul_small(10); | |
233 | minus.mul_small(10); | |
234 | plus.mul_small(10); | |
235 | } | |
236 | ||
237 | // rounding up happens when | |
238 | // i) only the rounding-up condition was triggered, or | |
239 | // ii) both conditions were triggered and tie breaking prefers rounding up. | |
240 | if up && (!down || *mant.mul_pow2(1) >= scale) { | |
241 | // if rounding up changes the length, the exponent should also change. | |
242 | // it seems that this condition is very hard to satisfy (possibly impossible), | |
243 | // but we are just being safe and consistent here. | |
244 | if let Some(c) = round_up(buf, i) { | |
245 | buf[i] = c; | |
246 | i += 1; | |
247 | k += 1; | |
248 | } | |
249 | } | |
250 | ||
251 | (i, k) | |
252 | } | |
253 | ||
254 | /// The exact and fixed mode implementation for Dragon. | |
255 | pub fn format_exact(d: &Decoded, buf: &mut [u8], limit: i16) -> (/*#digits*/ usize, /*exp*/ i16) { | |
256 | assert!(d.mant > 0); | |
257 | assert!(d.minus > 0); | |
258 | assert!(d.plus > 0); | |
259 | assert!(d.mant.checked_add(d.plus).is_some()); | |
260 | assert!(d.mant.checked_sub(d.minus).is_some()); | |
261 | ||
262 | // estimate `k_0` from original inputs satisfying `10^(k_0-1) < v <= 10^(k_0+1)`. | |
263 | let mut k = estimate_scaling_factor(d.mant, d.exp); | |
264 | ||
265 | // `v = mant / scale`. | |
266 | let mut mant = Big::from_u64(d.mant); | |
267 | let mut scale = Big::from_small(1); | |
268 | if d.exp < 0 { | |
269 | scale.mul_pow2(-d.exp as usize); | |
270 | } else { | |
271 | mant.mul_pow2(d.exp as usize); | |
272 | } | |
273 | ||
274 | // divide `mant` by `10^k`. now `scale / 10 < mant <= scale * 10`. | |
275 | if k >= 0 { | |
276 | mul_pow10(&mut scale, k as usize); | |
277 | } else { | |
278 | mul_pow10(&mut mant, -k as usize); | |
279 | } | |
280 | ||
281 | // fixup when `mant + plus >= scale`, where `plus / scale = 10^-buf.len() / 2`. | |
282 | // in order to keep the fixed-size bignum, we actually use `mant + floor(plus) >= scale`. | |
283 | // we are not actually modifying `scale`, since we can skip the initial multiplication instead. | |
284 | // again with the shortest algorithm, `d[0]` can be zero but will be eventually rounded up. | |
285 | if *div_2pow10(&mut scale.clone(), buf.len()).add(&mant) >= scale { | |
286 | // equivalent to scaling `scale` by 10 | |
287 | k += 1; | |
288 | } else { | |
289 | mant.mul_small(10); | |
290 | } | |
291 | ||
292 | // if we are working with the last-digit limitation, we need to shorten the buffer | |
293 | // before the actual rendering in order to avoid double rounding. | |
294 | // note that we have to enlarge the buffer again when rounding up happens! | |
295 | let mut len = if k < limit { | |
296 | // oops, we cannot even produce *one* digit. | |
297 | // this is possible when, say, we've got something like 9.5 and it's being rounded to 10. | |
298 | // we return an empty buffer, with an exception of the later rounding-up case | |
299 | // which occurs when `k == limit` and has to produce exactly one digit. | |
300 | 0 | |
301 | } else if ((k as i32 - limit as i32) as usize) < buf.len() { | |
302 | (k - limit) as usize | |
303 | } else { | |
304 | buf.len() | |
305 | }; | |
306 | ||
307 | if len > 0 { | |
308 | // cache `(2, 4, 8) * scale` for digit generation. | |
309 | // (this can be expensive, so do not calculate them when the buffer is empty.) | |
60c5eb7d XL |
310 | let mut scale2 = scale.clone(); |
311 | scale2.mul_pow2(1); | |
312 | let mut scale4 = scale.clone(); | |
313 | scale4.mul_pow2(2); | |
314 | let mut scale8 = scale.clone(); | |
315 | scale8.mul_pow2(3); | |
d9579d0f AL |
316 | |
317 | for i in 0..len { | |
60c5eb7d XL |
318 | if mant.is_zero() { |
319 | // following digits are all zeroes, we stop here | |
d9579d0f | 320 | // do *not* try to perform rounding! rather, fill remaining digits. |
60c5eb7d XL |
321 | for c in &mut buf[i..len] { |
322 | *c = b'0'; | |
323 | } | |
d9579d0f AL |
324 | return (len, k); |
325 | } | |
326 | ||
327 | let mut d = 0; | |
60c5eb7d XL |
328 | if mant >= scale8 { |
329 | mant.sub(&scale8); | |
330 | d += 8; | |
331 | } | |
332 | if mant >= scale4 { | |
333 | mant.sub(&scale4); | |
334 | d += 4; | |
335 | } | |
336 | if mant >= scale2 { | |
337 | mant.sub(&scale2); | |
338 | d += 2; | |
339 | } | |
340 | if mant >= scale { | |
341 | mant.sub(&scale); | |
342 | d += 1; | |
343 | } | |
d9579d0f AL |
344 | debug_assert!(mant < scale); |
345 | debug_assert!(d < 10); | |
346 | buf[i] = b'0' + d; | |
347 | mant.mul_small(10); | |
348 | } | |
349 | } | |
350 | ||
351 | // rounding up if we stop in the middle of digits | |
352 | // if the following digits are exactly 5000..., check the prior digit and try to | |
0731742a | 353 | // round to even (i.e., avoid rounding up when the prior digit is even). |
d9579d0f | 354 | let order = mant.cmp(scale.mul_small(5)); |
60c5eb7d XL |
355 | if order == Ordering::Greater |
356 | || (order == Ordering::Equal && (len == 0 || buf[len - 1] & 1 == 1)) | |
357 | { | |
d9579d0f AL |
358 | // if rounding up changes the length, the exponent should also change. |
359 | // but we've been requested a fixed number of digits, so do not alter the buffer... | |
360 | if let Some(c) = round_up(buf, len) { | |
361 | // ...unless we've been requested the fixed precision instead. | |
362 | // we also need to check that, if the original buffer was empty, | |
363 | // the additional digit can only be added when `k == limit` (edge case). | |
364 | k += 1; | |
365 | if k > limit && len < buf.len() { | |
366 | buf[len] = c; | |
367 | len += 1; | |
368 | } | |
369 | } | |
370 | } | |
371 | ||
372 | (len, k) | |
373 | } |