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