1 //! Compare the mantissa to the halfway representation of the float.
2 //!
3 //! Compares the actual significant digits of the mantissa to the
4 //! theoretical digits from `b+h`, scaled into the proper range.
5
6 use super::bignum::*;
7 use super::digit::*;
8 use super::exponent::*;
9 use super::float::*;
10 use super::math::*;
11 use super::num::*;
12 use super::rounding::*;
13 use crate::lib::{cmp, mem};
14
15 // MANTISSA
16
17 /// Parse the full mantissa into a big integer.
18 ///
19 /// Max digits is the maximum number of digits plus one.
parse_mantissa<'a, F, Iter>(mut iter: Iter) -> Bigint where F: Float, Iter: Iterator<Item = &'a u8> + Clone,20 fn parse_mantissa<'a, F, Iter>(mut iter: Iter) -> Bigint
21 where
22 F: Float,
23 Iter: Iterator<Item = &'a u8> + Clone,
24 {
25 // Main loop
26 let small_powers = POW10_LIMB;
27 let step = small_powers.len() - 2;
28 let max_digits = F::MAX_DIGITS - 1;
29 let mut counter = 0;
30 let mut value: Limb = 0;
31 let mut i: usize = 0;
32 let mut result = Bigint::default();
33
34 // Iteratively process all the data in the mantissa.
35 while let Some(&digit) = iter.next() {
36 // We've parsed the max digits using small values, add to bignum
37 if counter == step {
38 result.imul_small(small_powers[counter]);
39 result.iadd_small(value);
40 counter = 0;
41 value = 0;
42 }
43
44 value *= 10;
45 value += as_limb(to_digit(digit).unwrap());
46
47 i += 1;
48 counter += 1;
49 if i == max_digits {
50 break;
51 }
52 }
53
54 // We will always have a remainder, as long as we entered the loop
55 // once, or counter % step is 0.
56 if counter != 0 {
57 result.imul_small(small_powers[counter]);
58 result.iadd_small(value);
59 }
60
61 // If we have any remaining digits after the last value, we need
62 // to add a 1 after the rest of the array, it doesn't matter where,
63 // just move it up. This is good for the worst-possible float
64 // representation. We also need to return an index.
65 // Since we already trimmed trailing zeros, we know there has
66 // to be a non-zero digit if there are any left.
67 let is_consumed = iter.count() == 0;
68 if !is_consumed {
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]
b_extended<F: Float>(f: F) -> ExtendedFloat80 pub(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]
bh_extended<F: Float>(f: F) -> ExtendedFloat86 pub(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]
round_nearest_tie_even(fp: &mut ExtendedFloat, shift: i32, is_truncated: bool)99 fn 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.
large_atof<'a, F, Iter>(iter: Iter, exponent: i32) -> F where F: Float, Iter: Iterator<Item = &'a u8> + Clone,111 fn large_atof<'a, F, Iter>(iter: Iter, exponent: i32) -> F
112 where
113 F: Float,
114 Iter: Iterator<Item = &'a u8> + Clone,
115 {
116 let bits = mem::size_of::<u64>() * 8;
117
118 // Simple, we just need to multiply by the power of the radix.
119 // Now, we can calculate the mantissa and the exponent from this.
120 // The binary exponent is the binary exponent for the mantissa
121 // shifted to the hidden bit.
122 let mut bigmant = parse_mantissa::<F, _>(iter);
123 bigmant.imul_pow10(exponent.as_u32());
124
125 // Get the exact representation of the float from the big integer.
126 let (mant, is_truncated) = bigmant.hi64();
127 let exp = bigmant.bit_length().as_i32() - bits.as_i32();
128 let mut fp = ExtendedFloat {
129 mant,
130 exp,
131 };
132 fp.round_to_native::<F, _>(|fp, shift| round_nearest_tie_even(fp, shift, is_truncated));
133 into_float(fp)
134 }
135
136 /// Calculate the mantissa for a big integer with a negative exponent.
137 ///
138 /// This invokes the comparison with `b+h`.
small_atof<'a, F, Iter>(iter: Iter, exponent: i32, f: F) -> F where F: Float, Iter: Iterator<Item = &'a u8> + Clone,139 fn small_atof<'a, F, Iter>(iter: Iter, exponent: i32, f: F) -> F
140 where
141 F: Float,
142 Iter: Iterator<Item = &'a u8> + Clone,
143 {
144 // Get the significant digits and radix exponent for the real digits.
145 let mut real_digits = parse_mantissa::<F, _>(iter);
146 let real_exp = exponent;
147 debug_assert!(real_exp < 0);
148
149 // Get the significant digits and the binary exponent for `b+h`.
150 let theor = bh_extended(f);
151 let mut theor_digits = Bigint::from_u64(theor.mant);
152 let theor_exp = theor.exp;
153
154 // We need to scale the real digits and `b+h` digits to be the same
155 // order. We currently have `real_exp`, in `radix`, that needs to be
156 // shifted to `theor_digits` (since it is negative), and `theor_exp`
157 // to either `theor_digits` or `real_digits` as a power of 2 (since it
158 // may be positive or negative). Try to remove as many powers of 2
159 // as possible. All values are relative to `theor_digits`, that is,
160 // reflect the power you need to multiply `theor_digits` by.
161
162 // Can remove a power-of-two, since the radix is 10.
163 // Both are on opposite-sides of equation, can factor out a
164 // power of two.
165 //
166 // Example: 10^-10, 2^-10 -> ( 0, 10, 0)
167 // Example: 10^-10, 2^-15 -> (-5, 10, 0)
168 // Example: 10^-10, 2^-5 -> ( 5, 10, 0)
169 // Example: 10^-10, 2^5 -> (15, 10, 0)
170 let binary_exp = theor_exp - real_exp;
171 let halfradix_exp = -real_exp;
172 let radix_exp = 0;
173
174 // Carry out our multiplication.
175 if halfradix_exp != 0 {
176 theor_digits.imul_pow5(halfradix_exp.as_u32());
177 }
178 if radix_exp != 0 {
179 theor_digits.imul_pow10(radix_exp.as_u32());
180 }
181 if binary_exp > 0 {
182 theor_digits.imul_pow2(binary_exp.as_u32());
183 } else if binary_exp < 0 {
184 real_digits.imul_pow2((-binary_exp).as_u32());
185 }
186
187 // Compare real digits to theoretical digits and round the float.
188 match real_digits.compare(&theor_digits) {
189 cmp::Ordering::Greater => f.next_positive(),
190 cmp::Ordering::Less => f,
191 cmp::Ordering::Equal => f.round_positive_even(),
192 }
193 }
194
195 /// Calculate the exact value of the float.
196 ///
197 /// Notes:
198 /// The digits iterator must not have any trailing zeros (true for
199 /// `FloatState2`).
200 /// sci_exponent and digits.size_hint() must not overflow i32.
bhcomp<'a, F, Iter1, Iter2>(b: F, integer: Iter1, fraction: Iter2, exponent: i32) -> F where F: Float, Iter1: Iterator<Item = &'a u8> + Clone, Iter2: Iterator<Item = &'a u8> + Clone,201 pub(crate) fn bhcomp<'a, F, Iter1, Iter2>(b: F, integer: Iter1, fraction: Iter2, exponent: i32) -> F
202 where
203 F: Float,
204 Iter1: Iterator<Item = &'a u8> + Clone,
205 Iter2: Iterator<Item = &'a u8> + Clone,
206 {
207 // Calculate the number of integer digits and use that to determine
208 // where the significant digits start in the fraction.
209 let integer_digits = integer.clone().count();
210 let fraction_digits = fraction.clone().count();
211 let digits_start = match integer_digits {
212 0 => fraction.clone().take_while(|&x| x == &b'0').count(),
213 _ => 0,
214 };
215 let sci_exp = scientific_exponent(exponent, integer_digits, digits_start);
216 let count = F::MAX_DIGITS.min(integer_digits + fraction_digits - digits_start);
217 let scaled_exponent = sci_exp + 1 - count.as_i32();
218
219 // We have a finite conversions number of digits for base10.
220 // We need to then limit the iterator over that number of digits.
221 // Skip all leading zeros (can occur if fraction is empty).
222 let iter = integer.chain(fraction).skip(digits_start);
223
224 if scaled_exponent >= 0 {
225 large_atof(iter, scaled_exponent)
226 } else {
227 small_atof(iter, scaled_exponent, b)
228 }
229 }
230