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