1// Adapted from https://github.com/Alexhuszagh/rust-lexical.
2
3//! Compare the mantissa to the halfway representation of the float.
4//!
5//! Compares the actual significant digits of the mantissa to the
6//! theoretical digits from `b+h`, scaled into the proper range.
7
8use super::bignum::*;
9use super::digit::*;
10use super::exponent::*;
11use super::float::*;
12use super::math::*;
13use super::num::*;
14use super::rounding::*;
15use core::{cmp, mem};
16
17// MANTISSA
18
19/// Parse the full mantissa into a big integer.
20///
21/// Max digits is the maximum number of digits plus one.
22fn parse_mantissa<F>(integer: &[u8], fraction: &[u8]) -> Bigint
23where
24 F: Float,
25{
26 // Main loop
27 let small_powers = POW10_LIMB;
28 let step = small_powers.len() - 2;
29 let max_digits = F::MAX_DIGITS - 1;
30 let mut counter = 0;
31 let mut value: Limb = 0;
32 let mut i: usize = 0;
33 let mut result = Bigint::default();
34
35 // Iteratively process all the data in the mantissa.
36 for &digit in integer.iter().chain(fraction) {
37 // We've parsed the max digits using small values, add to bignum
38 if counter == step {
39 result.imul_small(small_powers[counter]);
40 result.iadd_small(value);
41 counter = 0;
42 value = 0;
43 }
44
45 value *= 10;
46 value += as_limb(to_digit(digit).unwrap());
47
48 i += 1;
49 counter += 1;
50 if i == max_digits {
51 break;
52 }
53 }
54
55 // We will always have a remainder, as long as we entered the loop
56 // once, or counter % step is 0.
57 if counter != 0 {
58 result.imul_small(small_powers[counter]);
59 result.iadd_small(value);
60 }
61
62 // If we have any remaining digits after the last value, we need
63 // to add a 1 after the rest of the array, it doesn't matter where,
64 // just move it up. This is good for the worst-possible float
65 // representation. We also need to return an index.
66 // Since we already trimmed trailing zeros, we know there has
67 // to be a non-zero digit if there are any left.
68 if i < integer.len() + fraction.len() {
69 result.imul_small(10);
70 result.iadd_small(1);
71 }
72
73 result
74}
75
76// FLOAT OPS
77
78/// Calculate `b` from a a representation of `b` as a float.
79#[inline]
80pub(super) fn b_extended<F: Float>(f: F) -> ExtendedFloat {
81 ExtendedFloat::from_float(f)
82}
83
84/// Calculate `b+h` from a a representation of `b` as a float.
85#[inline]
86pub(super) fn bh_extended<F: Float>(f: F) -> ExtendedFloat {
87 // None of these can overflow.
88 let b = b_extended(f);
89 ExtendedFloat {
90 mant: (b.mant << 1) + 1,
91 exp: b.exp - 1,
92 }
93}
94
95// ROUNDING
96
97/// Custom round-nearest, tie-event algorithm for bhcomp.
98#[inline]
99fn round_nearest_tie_even(fp: &mut ExtendedFloat, shift: i32, is_truncated: bool) {
100 let (mut is_above, mut is_halfway) = round_nearest(fp, shift);
101 if is_halfway && is_truncated {
102 is_above = true;
103 is_halfway = false;
104 }
105 tie_even(fp, is_above, is_halfway);
106}
107
108// BHCOMP
109
110/// Calculate the mantissa for a big integer with a positive exponent.
111fn large_atof<F>(mantissa: Bigint, exponent: i32) -> F
112where
113 F: Float,
114{
115 let bits = mem::size_of::<u64>() * 8;
116
117 // Simple, we just need to multiply by the power of the radix.
118 // Now, we can calculate the mantissa and the exponent from this.
119 // The binary exponent is the binary exponent for the mantissa
120 // shifted to the hidden bit.
121 let mut bigmant = mantissa;
122 bigmant.imul_pow10(exponent as u32);
123
124 // Get the exact representation of the float from the big integer.
125 let (mant, is_truncated) = bigmant.hi64();
126 let exp = bigmant.bit_length() as i32 - bits as i32;
127 let mut fp = ExtendedFloat { mant, exp };
128 fp.round_to_native::<F, _>(|fp, shift| round_nearest_tie_even(fp, shift, is_truncated));
129 into_float(fp)
130}
131
132/// Calculate the mantissa for a big integer with a negative exponent.
133///
134/// This invokes the comparison with `b+h`.
135fn small_atof<F>(mantissa: Bigint, exponent: i32, f: F) -> F
136where
137 F: Float,
138{
139 // Get the significant digits and radix exponent for the real digits.
140 let mut real_digits = mantissa;
141 let real_exp = exponent;
142 debug_assert!(real_exp < 0);
143
144 // Get the significant digits and the binary exponent for `b+h`.
145 let theor = bh_extended(f);
146 let mut theor_digits = Bigint::from_u64(theor.mant);
147 let theor_exp = theor.exp;
148
149 // We need to scale the real digits and `b+h` digits to be the same
150 // order. We currently have `real_exp`, in `radix`, that needs to be
151 // shifted to `theor_digits` (since it is negative), and `theor_exp`
152 // to either `theor_digits` or `real_digits` as a power of 2 (since it
153 // may be positive or negative). Try to remove as many powers of 2
154 // as possible. All values are relative to `theor_digits`, that is,
155 // reflect the power you need to multiply `theor_digits` by.
156
157 // Can remove a power-of-two, since the radix is 10.
158 // Both are on opposite-sides of equation, can factor out a
159 // power of two.
160 //
161 // Example: 10^-10, 2^-10 -> ( 0, 10, 0)
162 // Example: 10^-10, 2^-15 -> (-5, 10, 0)
163 // Example: 10^-10, 2^-5 -> ( 5, 10, 0)
164 // Example: 10^-10, 2^5 -> (15, 10, 0)
165 let binary_exp = theor_exp - real_exp;
166 let halfradix_exp = -real_exp;
167 let radix_exp = 0;
168
169 // Carry out our multiplication.
170 if halfradix_exp != 0 {
171 theor_digits.imul_pow5(halfradix_exp as u32);
172 }
173 if radix_exp != 0 {
174 theor_digits.imul_pow10(radix_exp as u32);
175 }
176 if binary_exp > 0 {
177 theor_digits.imul_pow2(binary_exp as u32);
178 } else if binary_exp < 0 {
179 real_digits.imul_pow2(-binary_exp as u32);
180 }
181
182 // Compare real digits to theoretical digits and round the float.
183 match real_digits.compare(&theor_digits) {
184 cmp::Ordering::Greater => f.next_positive(),
185 cmp::Ordering::Less => f,
186 cmp::Ordering::Equal => f.round_positive_even(),
187 }
188}
189
190/// Calculate the exact value of the float.
191///
192/// Note: fraction must not have trailing zeros.
193pub(crate) fn bhcomp<F>(b: F, integer: &[u8], mut fraction: &[u8], exponent: i32) -> F
194where
195 F: Float,
196{
197 // Calculate the number of integer digits and use that to determine
198 // where the significant digits start in the fraction.
199 let integer_digits = integer.len();
200 let fraction_digits = fraction.len();
201 let digits_start = if integer_digits == 0 {
202 let start = fraction.iter().take_while(|&x| *x == b'0').count();
203 fraction = &fraction[start..];
204 start
205 } else {
206 0
207 };
208 let sci_exp = scientific_exponent(exponent, integer_digits, digits_start);
209 let count = F::MAX_DIGITS.min(integer_digits + fraction_digits - digits_start);
210 let scaled_exponent = sci_exp + 1 - count as i32;
211
212 let mantissa = parse_mantissa::<F>(integer, fraction);
213 if scaled_exponent >= 0 {
214 large_atof(mantissa, scaled_exponent)
215 } else {
216 small_atof(mantissa, scaled_exponent, b)
217 }
218}
219