lexical_parse_float/slow.rs
1//! Slow, fallback cases where we cannot unambiguously round a float.
2//!
3//! This occurs when we cannot determine the exact representation using
4//! both the fast path (native) cases nor the Lemire/Bellerophon algorithms,
5//! and therefore must fallback to a slow, arbitrary-precision representation.
6
7#![doc(hidden)]
8
9use core::cmp;
10
11#[cfg(not(feature = "compact"))]
12use lexical_parse_integer::algorithm;
13use lexical_util::digit::char_to_valid_digit_const;
14#[cfg(feature = "radix")]
15use lexical_util::digit::digit_to_char_const;
16use lexical_util::format::NumberFormat;
17use lexical_util::iterator::{AsBytes, DigitsIter, Iter};
18use lexical_util::num::{AsPrimitive, Integer};
19
20#[cfg(feature = "radix")]
21use crate::bigint::Bigfloat;
22use crate::bigint::{Bigint, Limb};
23use crate::float::{extended_to_float, ExtendedFloat80, RawFloat};
24use crate::limits::{u32_power_limit, u64_power_limit};
25use crate::number::Number;
26use crate::shared;
27
28// ALGORITHM
29// ---------
30
31/// Parse the significant digits and biased, binary exponent of a float.
32///
33/// This is a fallback algorithm that uses a big-integer representation
34/// of the float, and therefore is considerably slower than faster
35/// approximations. However, it will always determine how to round
36/// the significant digits to the nearest machine float, allowing
37/// use to handle near half-way cases.
38///
39/// Near half-way cases are halfway between two consecutive machine floats.
40/// For example, the float `16777217.0` has a bitwise representation of
41/// `100000000000000000000000 1`. Rounding to a single-precision float,
42/// the trailing `1` is truncated. Using round-nearest, tie-even, any
43/// value above `16777217.0` must be rounded up to `16777218.0`, while
44/// any value before or equal to `16777217.0` must be rounded down
45/// to `16777216.0`. These near-halfway conversions therefore may require
46/// a large number of digits to unambiguously determine how to round.
47#[must_use]
48#[allow(clippy::unwrap_used)] // reason = "none is a developer error"
49pub fn slow_radix<F: RawFloat, const FORMAT: u128>(
50 num: Number,
51 fp: ExtendedFloat80,
52) -> ExtendedFloat80 {
53 // Ensure our preconditions are valid:
54 // 1. The significant digits are not shifted into place.
55 debug_assert!(fp.mant & (1 << 63) != 0, "number must be normalized");
56
57 let format = NumberFormat::<{ FORMAT }> {};
58
59 // This assumes the sign bit has already been parsed, and we're
60 // starting with the integer digits, and the float format has been
61 // correctly validated.
62 let sci_exp = scientific_exponent::<FORMAT>(&num);
63
64 // We have 3 major algorithms we use for this:
65 // 1. An algorithm with a finite number of digits and a positive exponent.
66 // 2. An algorithm with a finite number of digits and a negative exponent.
67 // 3. A fallback algorithm with a non-finite number of digits.
68
69 // In order for a float in radix `b` with a finite number of digits
70 // to have a finite representation in radix `y`, `b` should divide
71 // an integer power of `y`. This means for binary, all even radixes
72 // have finite representations, and all odd ones do not.
73 #[cfg(feature = "radix")]
74 {
75 if let Some(max_digits) = F::max_digits(format.radix()) {
76 // Can use our finite number of digit algorithm.
77 digit_comp::<F, FORMAT>(num, fp, sci_exp, max_digits)
78 } else {
79 // Fallback to infinite digits.
80 byte_comp::<F, FORMAT>(num, fp, sci_exp)
81 }
82 }
83
84 #[cfg(not(feature = "radix"))]
85 {
86 // Can use our finite number of digit algorithm.
87 let max_digits = F::max_digits(format.radix()).unwrap();
88 digit_comp::<F, FORMAT>(num, fp, sci_exp, max_digits)
89 }
90}
91
92/// Algorithm that generates the mantissa for a finite representation.
93///
94/// For a positive exponent relative to the significant digits, this
95/// is just a multiplication by an exponent power. For a negative
96/// exponent relative to the significant digits, we scale the real
97/// digits to the theoretical digits for `b` and determine if we
98/// need to round-up.
99#[must_use]
100#[inline(always)]
101#[allow(clippy::cast_possible_wrap)] // reason = "the value range is [-324, 308]"
102pub fn digit_comp<F: RawFloat, const FORMAT: u128>(
103 num: Number,
104 fp: ExtendedFloat80,
105 sci_exp: i32,
106 max_digits: usize,
107) -> ExtendedFloat80 {
108 let (bigmant, digits) = parse_mantissa::<FORMAT>(num, max_digits);
109 // This can't underflow, since `digits` is at most `max_digits`.
110 let exponent = sci_exp + 1 - digits as i32;
111 if exponent >= 0 {
112 positive_digit_comp::<F, FORMAT>(bigmant, exponent)
113 } else {
114 negative_digit_comp::<F, FORMAT>(bigmant, fp, exponent)
115 }
116}
117
118/// Generate the significant digits with a positive exponent relative to
119/// mantissa.
120#[must_use]
121#[inline(always)]
122#[allow(clippy::unwrap_used)] // reason = "none is a developer error"
123#[allow(clippy::cast_possible_wrap)] // reason = "can't wrap in practice: max is ~1000 limbs"
124#[allow(clippy::missing_inline_in_public_items)] // reason = "only public for testing"
125pub fn positive_digit_comp<F: RawFloat, const FORMAT: u128>(
126 mut bigmant: Bigint,
127 exponent: i32,
128) -> ExtendedFloat80 {
129 let format = NumberFormat::<{ FORMAT }> {};
130
131 // Simple, we just need to multiply by the power of the radix.
132 // Now, we can calculate the mantissa and the exponent from this.
133 // The binary exponent is the binary exponent for the mantissa
134 // shifted to the hidden bit.
135 bigmant.pow(format.radix(), exponent as u32).unwrap();
136
137 // Get the exact representation of the float from the big integer.
138 // hi64 checks **all** the remaining bits after the mantissa,
139 // so it will check if **any** truncated digits exist.
140 let (mant, is_truncated) = bigmant.hi64();
141 let exp = bigmant.bit_length() as i32 - 64 + F::EXPONENT_BIAS;
142 let mut fp = ExtendedFloat80 {
143 mant,
144 exp,
145 };
146
147 // Shift the digits into position and determine if we need to round-up.
148 shared::round::<F, _>(&mut fp, |f, s| {
149 shared::round_nearest_tie_even(f, s, |is_odd, is_halfway, is_above| {
150 is_above || (is_halfway && is_truncated) || (is_odd && is_halfway)
151 });
152 });
153 fp
154}
155
156/// Generate the significant digits with a negative exponent relative to
157/// mantissa.
158///
159/// This algorithm is quite simple: we have the significant digits `m1 * b^N1`,
160/// where `m1` is the bigint mantissa, `b` is the radix, and `N1` is the radix
161/// exponent. We then calculate the theoretical representation of `b+h`, which
162/// is `m2 * 2^N2`, where `m2` is the bigint mantissa and `N2` is the binary
163/// exponent. If we had infinite, efficient floating precision, this would be
164/// equal to `m1 / b^-N1` and then compare it to `m2 * 2^N2`.
165///
166/// Since we cannot divide and keep precision, we must multiply the other:
167/// if we want to do `m1 / b^-N1 >= m2 * 2^N2`, we can do
168/// `m1 >= m2 * b^-N1 * 2^N2` Going to the decimal case, we can show and example
169/// and simplify this further: `m1 >= m2 * 2^N2 * 10^-N1`. Since we can remove
170/// a power-of-two, this is `m1 >= m2 * 2^(N2 - N1) * 5^-N1`. Therefore, if
171/// `N2 - N1 > 0`, we need have `m1 >= m2 * 2^(N2 - N1) * 5^-N1`, otherwise,
172/// we have `m1 * 2^(N1 - N2) >= m2 * 5^-N1`, where the resulting exponents
173/// are all positive.
174///
175/// This allows us to compare both floats using integers efficiently
176/// without any loss of precision.
177#[must_use]
178#[inline(always)]
179#[allow(clippy::match_bool)] // reason = "simplifies documentation"
180#[allow(clippy::unwrap_used)] // reason = "unwrap panics if a developer error"
181#[allow(clippy::comparison_chain)] // reason = "logically different conditions for algorithm"
182#[allow(clippy::missing_inline_in_public_items)] // reason = "only exposed for unittesting"
183pub fn negative_digit_comp<F: RawFloat, const FORMAT: u128>(
184 bigmant: Bigint,
185 mut fp: ExtendedFloat80,
186 exponent: i32,
187) -> ExtendedFloat80 {
188 // Ensure our preconditions are valid:
189 // 1. The significant digits are not shifted into place.
190 debug_assert!(fp.mant & (1 << 63) != 0, "the significant digits must be normalized");
191
192 let format = NumberFormat::<FORMAT> {};
193 let radix = format.radix();
194
195 // Get the significant digits and radix exponent for the real digits.
196 let mut real_digits = bigmant;
197 let real_exp = exponent;
198 debug_assert!(real_exp < 0, "algorithm only works with negative numbers");
199
200 // Round down our extended-precision float and calculate `b`.
201 let mut b = fp;
202 shared::round::<F, _>(&mut b, shared::round_down);
203 let b = extended_to_float::<F>(b);
204
205 // Get the significant digits and the binary exponent for `b+h`.
206 let theor = bh(b);
207 let mut theor_digits = Bigint::from_u64(theor.mant);
208 let theor_exp = theor.exp;
209
210 // We need to scale the real digits and `b+h` digits to be the same
211 // order. We currently have `real_exp`, in `radix`, that needs to be
212 // shifted to `theor_digits` (since it is negative), and `theor_exp`
213 // to either `theor_digits` or `real_digits` as a power of 2 (since it
214 // may be positive or negative). Try to remove as many powers of 2
215 // as possible. All values are relative to `theor_digits`, that is,
216 // reflect the power you need to multiply `theor_digits` by.
217 let (binary_exp, halfradix_exp, radix_exp) = match radix.is_even() {
218 // Can remove a power-of-two.
219 // Both are on opposite-sides of equation, can factor out a
220 // power of two.
221 //
222 // Example: 10^-10, 2^-10 -> ( 0, 10, 0)
223 // Example: 10^-10, 2^-15 -> (-5, 10, 0)
224 // Example: 10^-10, 2^-5 -> ( 5, 10, 0)
225 // Example: 10^-10, 2^5 -> (15, 10, 0)
226 true => (theor_exp - real_exp, -real_exp, 0),
227 // Cannot remove a power-of-two.
228 false => (theor_exp, 0, -real_exp),
229 };
230
231 if halfradix_exp != 0 {
232 theor_digits.pow(radix / 2, halfradix_exp as u32).unwrap();
233 }
234 if radix_exp != 0 {
235 theor_digits.pow(radix, radix_exp as u32).unwrap();
236 }
237 if binary_exp > 0 {
238 theor_digits.pow(2, binary_exp as u32).unwrap();
239 } else if binary_exp < 0 {
240 real_digits.pow(2, (-binary_exp) as u32).unwrap();
241 }
242
243 // Compare our theoretical and real digits and round nearest, tie even.
244 let ord = real_digits.data.cmp(&theor_digits.data);
245 shared::round::<F, _>(&mut fp, |f, s| {
246 shared::round_nearest_tie_even(f, s, |is_odd, _, _| {
247 // Can ignore `is_halfway` and `is_above`, since those were
248 // calculates using less significant digits.
249 match ord {
250 cmp::Ordering::Greater => true,
251 cmp::Ordering::Less => false,
252 cmp::Ordering::Equal if is_odd => true,
253 cmp::Ordering::Equal => false,
254 }
255 });
256 });
257 fp
258}
259
260/// Try to parse 8 digits at a time.
261///
262/// - `format` - The numerical format specification as a packed 128-bit integer
263/// - `iter` - An iterator over all bytes in the buffer
264/// - `value` - The currently parsed value.
265/// - `count` - The total number of parsed digits
266/// - `counter` - The number of parsed digits since creating the current u32
267/// - `step` - The maximum number of digits for the radix that can fit in a u32.
268/// - `max_digits` - The maximum number of digits that can affect floating-point
269/// rounding.
270#[cfg(not(feature = "compact"))]
271macro_rules! try_parse_8digits {
272 (
273 $format:ident,
274 $iter:ident,
275 $value:ident,
276 $count:ident,
277 $counter:ident,
278 $step:ident,
279 $max_digits:ident
280 ) => {{
281 let format = NumberFormat::<$format> {};
282 let radix = format.radix() as Limb;
283
284 // Try 8-digit optimizations.
285 if can_try_parse_multidigit!($iter, radix) {
286 debug_assert!(radix < 16);
287 let radix8 = format.radix8() as Limb;
288 while $step - $counter >= 8 && $max_digits - $count >= 8 {
289 if let Some(v) = algorithm::try_parse_8digits::<Limb, _, FORMAT>(&mut $iter) {
290 $value = $value.wrapping_mul(radix8).wrapping_add(v);
291 $counter += 8;
292 $count += 8;
293 } else {
294 break;
295 }
296 }
297 }
298 }};
299}
300
301/// Add a digit to the temporary value.
302///
303/// - `c` - The character to convert to a digit.
304/// - `value` - The currently parsed value.
305/// - `count` - The total number of parsed digits
306/// - `counter` - The number of parsed digits since creating the current u32
307macro_rules! add_digit {
308 ($c:ident, $radix:ident, $value:ident, $counter:ident, $count:ident) => {{
309 let digit = char_to_valid_digit_const($c, $radix);
310 $value *= $radix as Limb;
311 $value += digit as Limb;
312
313 // Increment our counters.
314 $counter += 1;
315 $count += 1;
316 }};
317}
318
319/// Add a temporary value to our mantissa.
320///
321/// - `format` - The numerical format specification as a packed 128-bit integer
322/// - `result` - The big integer,
323/// - `power` - The power to scale the big integer by.
324/// - `value` - The value to add to the big integer,
325/// - `counter` - The number of parsed digits since creating the current u32
326macro_rules! add_temporary {
327 // Multiply by the small power and add the native value.
328 (@mul $result:ident, $power:expr, $value:expr) => {
329 $result.data.mul_small($power).unwrap();
330 $result.data.add_small($value).unwrap();
331 };
332
333 // Add a temporary where we won't read the counter results internally.
334 (@end $format:ident, $result:ident, $counter:ident, $value:ident) => {
335 if $counter != 0 {
336 let small_power = f64::int_pow_fast_path($counter, $format.radix());
337 add_temporary!(@mul $result, small_power as Limb, $value);
338 }
339 };
340
341 // Add the maximum native value.
342 (@max $format:ident, $result:ident, $counter:ident, $value:ident, $max:ident) => {
343 add_temporary!(@mul $result, $max, $value);
344 $counter = 0;
345 $value = 0;
346 };
347}
348
349/// Round-up a truncated value.
350///
351/// - `format` - The numerical format specification as a packed 128-bit integer
352/// - `result` - The big integer,
353/// - `count` - The total number of parsed digits
354macro_rules! round_up_truncated {
355 ($format:ident, $result:ident, $count:ident) => {{
356 // Need to round-up.
357 // Can't just add 1, since this can accidentally round-up
358 // values to a halfway point, which can cause invalid results.
359 add_temporary!(@mul $result, $format.radix() as Limb, 1);
360 $count += 1;
361 }};
362}
363
364/// Check and round-up the fraction if any non-zero digits exist.
365///
366/// - `format` - The numerical format specification as a packed 128-bit integer
367/// - `iter` - An iterator over all bytes in the buffer
368/// - `result` - The big integer,
369/// - `count` - The total number of parsed digits
370macro_rules! round_up_nonzero {
371 ($format:ident, $iter:expr, $result:ident, $count:ident) => {{
372 // NOTE: All digits must already be valid.
373 let mut iter = $iter;
374
375 // First try reading 8-digits at a time.
376 if iter.is_contiguous() {
377 while let Some(value) = iter.peek_u64() {
378 // SAFETY: safe since we have at least 8 bytes in the buffer.
379 unsafe { iter.step_by_unchecked(8) };
380 if value != 0x3030_3030_3030_3030 {
381 // Have non-zero digits, exit early.
382 round_up_truncated!($format, $result, $count);
383 return ($result, $count);
384 }
385 }
386 }
387
388 for &digit in iter {
389 if digit != b'0' {
390 round_up_truncated!($format, $result, $count);
391 return ($result, $count);
392 }
393 }
394 }};
395}
396
397/// Parse the full mantissa into a big integer.
398///
399/// Returns the parsed mantissa and the number of digits in the mantissa.
400/// The max digits is the maximum number of digits plus one.
401#[must_use]
402#[allow(clippy::cognitive_complexity)] // reason = "complexity broken into macros"
403#[allow(clippy::missing_inline_in_public_items)] // reason = "only public for testing"
404pub fn parse_mantissa<const FORMAT: u128>(num: Number, max_digits: usize) -> (Bigint, usize) {
405 let format = NumberFormat::<FORMAT> {};
406 let radix = format.radix();
407
408 // Iteratively process all the data in the mantissa.
409 // We do this via small, intermediate values which once we reach
410 // the maximum number of digits we can process without overflow,
411 // we add the temporary to the big integer.
412 let mut counter: usize = 0;
413 let mut count: usize = 0;
414 let mut value: Limb = 0;
415 let mut result = Bigint::new();
416
417 // Now use our pre-computed small powers iteratively.
418 let step = if Limb::BITS == 32 {
419 u32_power_limit(format.radix())
420 } else {
421 u64_power_limit(format.radix())
422 } as usize;
423 let max_native = (format.radix() as Limb).pow(step as u32);
424
425 // Process the integer digits.
426 let mut integer = num.integer.bytes::<FORMAT>();
427 let mut integer_iter = integer.integer_iter();
428 integer_iter.skip_zeros();
429 'integer: loop {
430 #[cfg(not(feature = "compact"))]
431 try_parse_8digits!(FORMAT, integer_iter, value, count, counter, step, max_digits);
432
433 // Parse a digit at a time, until we reach step.
434 while counter < step && count < max_digits {
435 if let Some(&c) = integer_iter.next() {
436 add_digit!(c, radix, value, counter, count);
437 } else {
438 break 'integer;
439 }
440 }
441
442 // Check if we've exhausted our max digits.
443 if count == max_digits {
444 // Need to check if we're truncated, and round-up accordingly.
445 // SAFETY: safe since `counter <= step`.
446 add_temporary!(@end format, result, counter, value);
447 round_up_nonzero!(format, integer_iter, result, count);
448 if let Some(fraction) = num.fraction {
449 let mut fraction = fraction.bytes::<FORMAT>();
450 round_up_nonzero!(format, fraction.fraction_iter(), result, count);
451 }
452 return (result, count);
453 } else {
454 // Add our temporary from the loop.
455 // SAFETY: safe since `counter <= step`.
456 add_temporary!(@max format, result, counter, value, max_native);
457 }
458 }
459
460 // Process the fraction digits.
461 if let Some(fraction) = num.fraction {
462 let mut fraction = fraction.bytes::<FORMAT>();
463 let mut fraction_iter = fraction.integer_iter();
464 if count == 0 {
465 // No digits added yet, can skip leading fraction zeros too.
466 fraction_iter.skip_zeros();
467 }
468 'fraction: loop {
469 #[cfg(not(feature = "compact"))]
470 try_parse_8digits!(FORMAT, fraction_iter, value, count, counter, step, max_digits);
471
472 // Parse a digit at a time, until we reach step.
473 while counter < step && count < max_digits {
474 if let Some(&c) = fraction_iter.next() {
475 add_digit!(c, radix, value, counter, count);
476 } else {
477 break 'fraction;
478 }
479 }
480
481 // Check if we've exhausted our max digits.
482 if count == max_digits {
483 // SAFETY: safe since `counter <= step`.
484 add_temporary!(@end format, result, counter, value);
485 round_up_nonzero!(format, fraction_iter, result, count);
486 return (result, count);
487 } else {
488 // Add our temporary from the loop.
489 // SAFETY: safe since `counter <= step`.
490 add_temporary!(@max format, result, counter, value, max_native);
491 }
492 }
493 }
494
495 // We will always have a remainder, as long as we entered the loop
496 // once, or counter % step is 0.
497 // SAFETY: safe since `counter <= step`.
498 add_temporary!(@end format, result, counter, value);
499
500 (result, count)
501}
502
503/// Compare actual integer digits to the theoretical digits.
504///
505/// - `iter` - An iterator over all bytes in the buffer
506/// - `num` - The actual digits of the real floating point number.
507/// - `den` - The theoretical digits created by `b+h` to determine if `b` or
508/// `b+1`
509#[cfg(feature = "radix")]
510macro_rules! integer_compare {
511 ($iter:ident, $num:ident, $den:ident, $radix:ident) => {{
512 // Compare the integer digits.
513 while !$num.data.is_empty() {
514 // All digits **must** be valid.
515 let actual = match $iter.next() {
516 Some(&v) => v,
517 // Could have hit the decimal point.
518 _ => break,
519 };
520 let rem = $num.data.quorem(&$den.data) as u32;
521 let expected = digit_to_char_const(rem, $radix);
522 $num.data.mul_small($radix as Limb).unwrap();
523 if actual < expected {
524 return cmp::Ordering::Less;
525 } else if actual > expected {
526 return cmp::Ordering::Greater;
527 }
528 }
529
530 // Still have integer digits, check if any are non-zero.
531 if $num.data.is_empty() {
532 for &digit in $iter {
533 if digit != b'0' {
534 return cmp::Ordering::Greater;
535 }
536 }
537 }
538 }};
539}
540
541/// Compare actual fraction digits to the theoretical digits.
542///
543/// - `iter` - An iterator over all bytes in the buffer
544/// - `num` - The actual digits of the real floating point number.
545/// - `den` - The theoretical digits created by `b+h` to determine if `b` or
546/// `b+1`
547#[cfg(feature = "radix")]
548macro_rules! fraction_compare {
549 ($iter:ident, $num:ident, $den:ident, $radix:ident) => {{
550 // Compare the fraction digits.
551 // We can only be here if we hit a decimal point.
552 while !$num.data.is_empty() {
553 // All digits **must** be valid.
554 let actual = match $iter.next() {
555 Some(&v) => v,
556 // No more actual digits, or hit the exponent.
557 _ => return cmp::Ordering::Less,
558 };
559 let rem = $num.data.quorem(&$den.data) as u32;
560 let expected = digit_to_char_const(rem, $radix);
561 $num.data.mul_small($radix as Limb).unwrap();
562 if actual < expected {
563 return cmp::Ordering::Less;
564 } else if actual > expected {
565 return cmp::Ordering::Greater;
566 }
567 }
568
569 // Still have fraction digits, check if any are non-zero.
570 for &digit in $iter {
571 if digit != b'0' {
572 return cmp::Ordering::Greater;
573 }
574 }
575 }};
576}
577
578/// Compare theoretical digits to halfway point from theoretical digits.
579///
580/// Generates a float representing the halfway point, and generates
581/// theoretical digits as bytes, and compares the generated digits to
582/// the actual input.
583///
584/// Compares the known string to theoretical digits generated on the
585/// fly for `b+h`, where a string representation of a float is between
586/// `b` and `b+u`, where `b+u` is 1 unit in the least-precision. Therefore,
587/// the string must be close to `b+h`.
588///
589/// Adapted from "Bigcomp: Deciding Truncated, Near Halfway Conversions",
590/// available [here](https://www.exploringbinary.com/bigcomp-deciding-truncated-near-halfway-conversions/).
591#[cfg(feature = "radix")]
592#[allow(clippy::unwrap_used)] // reason = "none is a developer error due to shl overflow"
593#[allow(clippy::comparison_chain)] // reason = "logically different conditions for algorithm"
594pub fn byte_comp<F: RawFloat, const FORMAT: u128>(
595 number: Number,
596 mut fp: ExtendedFloat80,
597 sci_exp: i32,
598) -> ExtendedFloat80 {
599 // Ensure our preconditions are valid:
600 // 1. The significant digits are not shifted into place.
601 debug_assert!(fp.mant & (1 << 63) != 0);
602
603 let format = NumberFormat::<FORMAT> {};
604
605 // Round down our extended-precision float and calculate `b`.
606 let mut b = fp;
607 shared::round::<F, _>(&mut b, shared::round_down);
608 let b = extended_to_float::<F>(b);
609
610 // Calculate `b+h` to create a ratio for our theoretical digits.
611 let theor = Bigfloat::from_float(bh::<F>(b));
612
613 // Now, create a scaling factor for the digit count.
614 let mut factor = Bigfloat::from_u32(1);
615 factor.pow(format.radix(), sci_exp.unsigned_abs()).unwrap();
616 let mut num: Bigfloat;
617 let mut den: Bigfloat;
618
619 if sci_exp < 0 {
620 // Need to have the basen factor be the numerator, and the `fp`
621 // be the denominator. Since we assumed that `theor` was the numerator,
622 // if it's the denominator, we need to multiply it into the numerator.
623 num = factor;
624 num.data *= &theor.data;
625 den = Bigfloat::from_u32(1);
626 den.exp = -theor.exp;
627 } else {
628 num = theor;
629 den = factor;
630 }
631
632 // Scale the denominator so it has the number of bits
633 // in the radix as the number of leading zeros.
634 let wlz = integral_binary_factor(format.radix());
635 let nlz = den.leading_zeros().wrapping_sub(wlz) & (32 - 1);
636 if nlz != 0 {
637 den.shl_bits(nlz as usize).unwrap();
638 den.exp -= nlz as i32;
639 }
640
641 // Need to scale the numerator or denominator to the same value.
642 // We don't want to shift the denominator, so...
643 let diff = den.exp - num.exp;
644 let shift = diff.unsigned_abs() as usize;
645 if diff < 0 {
646 // Need to shift the numerator left.
647 num.shl(shift).unwrap();
648 num.exp -= shift as i32;
649 } else if diff > 0 {
650 // Need to shift denominator left, go by a power of Limb::BITS.
651 // After this, the numerator will be non-normalized, and the
652 // denominator will be normalized. We need to add one to the
653 // quotient,since we're calculating the ceiling of the divmod.
654 let (q, r) = shift.ceil_divmod(Limb::BITS as usize);
655 let r = -r;
656 if r != 0 {
657 num.shl_bits(r as usize).unwrap();
658 num.exp -= r;
659 }
660 if q != 0 {
661 den.shl_limbs(q).unwrap();
662 den.exp -= Limb::BITS as i32 * q as i32;
663 }
664 }
665
666 // Compare our theoretical and real digits and round nearest, tie even.
667 let ord = compare_bytes::<FORMAT>(number, num, den);
668 shared::round::<F, _>(&mut fp, |f, s| {
669 shared::round_nearest_tie_even(f, s, |is_odd, _, _| {
670 // Can ignore `is_halfway` and `is_above`, since those were
671 // calculates using less significant digits.
672 match ord {
673 cmp::Ordering::Greater => true,
674 cmp::Ordering::Less => false,
675 cmp::Ordering::Equal if is_odd => true,
676 cmp::Ordering::Equal => false,
677 }
678 });
679 });
680 fp
681}
682
683/// Compare digits between the generated values the ratio and the actual view.
684///
685/// - `number` - The representation of the float as a big number, with the
686/// parsed digits.
687/// - `num` - The actual digits of the real floating point number.
688/// - `den` - The theoretical digits created by `b+h` to determine if `b` or
689/// `b+1`
690#[cfg(feature = "radix")]
691#[allow(clippy::unwrap_used)] // reason = "none is a developer error due to a missing fraction"
692pub fn compare_bytes<const FORMAT: u128>(
693 number: Number,
694 mut num: Bigfloat,
695 den: Bigfloat,
696) -> cmp::Ordering {
697 let format = NumberFormat::<FORMAT> {};
698 let radix = format.radix();
699
700 // Now need to compare the theoretical digits. First, I need to trim
701 // any leading zeros, and will also need to ignore trailing ones.
702 let mut integer = number.integer.bytes::<{ FORMAT }>();
703 let mut integer_iter = integer.integer_iter();
704 integer_iter.skip_zeros();
705 if integer_iter.is_buffer_empty() {
706 // Cannot be empty, since we must have at least **some** significant digits.
707 let mut fraction = number.fraction.unwrap().bytes::<{ FORMAT }>();
708 let mut fraction_iter = fraction.fraction_iter();
709 fraction_iter.skip_zeros();
710 fraction_compare!(fraction_iter, num, den, radix);
711 } else {
712 integer_compare!(integer_iter, num, den, radix);
713 if let Some(fraction) = number.fraction {
714 let mut fraction = fraction.bytes::<{ FORMAT }>();
715 let mut fraction_iter = fraction.fraction_iter();
716 fraction_compare!(fraction_iter, num, den, radix);
717 } else if !num.data.is_empty() {
718 // We had more theoretical digits, but no more actual digits.
719 return cmp::Ordering::Less;
720 }
721 }
722
723 // Exhausted both, must be equal.
724 cmp::Ordering::Equal
725}
726
727// SCALING
728// -------
729
730/// Calculate the scientific exponent from a `Number` value.
731/// Any other attempts would require slowdowns for faster algorithms.
732#[must_use]
733#[inline(always)]
734pub fn scientific_exponent<const FORMAT: u128>(num: &Number) -> i32 {
735 // This has the significant digits and exponent relative to those
736 // digits: therefore, we just need to scale to mantissa to `[1, radix)`.
737 // This doesn't need to be very fast.
738 let format = NumberFormat::<FORMAT> {};
739
740 // Use power reduction to make this faster: we need at least
741 // `F::MANTISSA_SIZE` bits, so we must have at least radix^4 digits.
742 // IF we're using base 3, we can have at most 11 divisions, and
743 // base 36, at most ~4. So, this is reasonably efficient.
744 let radix = format.radix() as u64;
745 let radix2 = radix * radix;
746 let radix4 = radix2 * radix2;
747 let mut mantissa = num.mantissa;
748 let mut exponent = num.exponent;
749 while mantissa >= radix4 {
750 mantissa /= radix4;
751 exponent += 4;
752 }
753 while mantissa >= radix2 {
754 mantissa /= radix2;
755 exponent += 2;
756 }
757 while mantissa >= radix {
758 mantissa /= radix;
759 exponent += 1;
760 }
761 exponent as i32
762}
763
764/// Calculate `b` from a a representation of `b` as a float.
765#[must_use]
766#[inline(always)]
767pub fn b<F: RawFloat>(float: F) -> ExtendedFloat80 {
768 ExtendedFloat80 {
769 mant: float.mantissa().as_u64(),
770 exp: float.exponent(),
771 }
772}
773
774/// Calculate `b+h` from a a representation of `b` as a float.
775#[must_use]
776#[inline(always)]
777pub fn bh<F: RawFloat>(float: F) -> ExtendedFloat80 {
778 let fp = b(float);
779 ExtendedFloat80 {
780 mant: (fp.mant << 1) + 1,
781 exp: fp.exp - 1,
782 }
783}
784
785// NOTE: There will never be binary factors here.
786
787/// Calculate the integral ceiling of the binary factor from a basen number.
788#[must_use]
789#[inline(always)]
790#[cfg(feature = "radix")]
791pub const fn integral_binary_factor(radix: u32) -> u32 {
792 match radix {
793 3 => 2,
794 5 => 3,
795 6 => 3,
796 7 => 3,
797 9 => 4,
798 10 => 4,
799 11 => 4,
800 12 => 4,
801 13 => 4,
802 14 => 4,
803 15 => 4,
804 17 => 5,
805 18 => 5,
806 19 => 5,
807 20 => 5,
808 21 => 5,
809 22 => 5,
810 23 => 5,
811 24 => 5,
812 25 => 5,
813 26 => 5,
814 27 => 5,
815 28 => 5,
816 29 => 5,
817 30 => 5,
818 31 => 5,
819 33 => 6,
820 34 => 6,
821 35 => 6,
822 36 => 6,
823 // Invalid radix
824 _ => 0,
825 }
826}
827
828/// Calculate the integral ceiling of the binary factor from a basen number.
829#[must_use]
830#[inline(always)]
831#[cfg(not(feature = "radix"))]
832pub const fn integral_binary_factor(radix: u32) -> u32 {
833 match radix {
834 10 => 4,
835 // Invalid radix
836 _ => 0,
837 }
838}