malachite_float/basic/ulp.rs
1// Copyright © 2025 Mikhail Hogrefe
2//
3// This file is part of Malachite.
4//
5// Malachite is free software: you can redistribute it and/or modify it under the terms of the GNU
6// Lesser General Public License (LGPL) as published by the Free Software Foundation; either version
7// 3 of the License, or (at your option) any later version. See <https://www.gnu.org/licenses/>.
8
9use crate::InnerFloat::Finite;
10use crate::{Float, significand_bits};
11use malachite_base::num::arithmetic::traits::{DivisibleByPowerOf2, NegAssign, PowerOf2};
12use malachite_base::num::basic::integers::PrimitiveInt;
13use malachite_base::num::basic::traits::{Infinity, Zero};
14use malachite_base::num::conversion::traits::{ExactFrom, WrappingFrom};
15use malachite_base::num::logic::traits::SignificantBits;
16use malachite_nz::platform::Limb;
17
18impl Float {
19 /// Gets a [`Float`]'s ulp (unit in last place, or unit of least precision).
20 ///
21 /// If the [`Float`] is positive, its ulp is the distance to the next-largest [`Float`] with the
22 /// same precision; if it is negative, the next-smallest. (This definition works even if the
23 /// [`Float`] is the largest in its binade. If the [`Float`] is the largest in its binade and
24 /// has the maximum exponent, we can define its ulp to be the distance to the next-smallest
25 /// [`Float`] with the same precision if positive, and to the next-largest [`Float`] with the
26 /// same precision if negative.)
27 ///
28 /// If the [`Float`] is NaN, infinite, or zero, then `None` is returned.
29 ///
30 /// This function does not overflow or underflow, technically. But it is possible that a
31 /// [`Float`]'s ulp is too small to represent, for example if the [`Float`] has the minimum
32 /// exponent and its precision is greater than 1, or if the precision is extremely large in
33 /// general. In such cases, `None` is returned.
34 ///
35 /// $$
36 /// f(\text{NaN}) = f(\pm\infty) = f(\pm 0.0) = \text{None},
37 /// $$
38 ///
39 /// and, if $x$ is finite and nonzero,
40 ///
41 /// $$
42 /// f(x) = \operatorname{Some}(2^{\lfloor \log_2 x \rfloor-p+1}),
43 /// $$
44 /// where $p$ is the precision of $x$.
45 ///
46 /// # Worst-case complexity
47 /// $T(n) = O(n)$
48 ///
49 /// $M(n) = O(n)$
50 ///
51 /// where $T$ is time, $M$ is additional memory, and $n$ is `self.significant_bits()`.
52 ///
53 /// # Examples
54 /// ```
55 /// use malachite_base::num::arithmetic::traits::PowerOf2;
56 /// use malachite_base::num::basic::traits::{Infinity, NaN, NegativeOne, One, Zero};
57 /// use malachite_float::Float;
58 ///
59 /// assert_eq!(Float::NAN.ulp(), None);
60 /// assert_eq!(Float::INFINITY.ulp(), None);
61 /// assert_eq!(Float::ZERO.ulp(), None);
62 ///
63 /// let s = Float::ONE.ulp().map(|x| x.to_string());
64 /// assert_eq!(s.as_ref().map(|s| s.as_str()), Some("1.0"));
65 ///
66 /// let s = Float::one_prec(100).ulp().map(|x| x.to_string());
67 /// assert_eq!(s.as_ref().map(|s| s.as_str()), Some("2.0e-30"));
68 ///
69 /// let s = Float::from(std::f64::consts::PI)
70 /// .ulp()
71 /// .map(|x| x.to_string());
72 /// assert_eq!(s.as_ref().map(|s| s.as_str()), Some("4.0e-15"));
73 ///
74 /// let s = Float::power_of_2(100u64).ulp().map(|x| x.to_string());
75 /// assert_eq!(s.as_ref().map(|s| s.as_str()), Some("1.0e30"));
76 ///
77 /// let s = Float::power_of_2(-100i64).ulp().map(|x| x.to_string());
78 /// assert_eq!(s.as_ref().map(|s| s.as_str()), Some("8.0e-31"));
79 ///
80 /// let s = Float::NEGATIVE_ONE.ulp().map(|x| x.to_string());
81 /// assert_eq!(s.as_ref().map(|s| s.as_str()), Some("1.0"));
82 /// ```
83 pub fn ulp(&self) -> Option<Float> {
84 match self {
85 Float(Finite {
86 exponent,
87 precision,
88 ..
89 }) => {
90 let ulp_exponent =
91 i64::from(*exponent).checked_sub(i64::try_from(*precision).ok()?)?;
92 if i32::try_from(ulp_exponent).ok()? >= Float::MIN_EXPONENT - 1 {
93 Some(Float::power_of_2(ulp_exponent))
94 } else {
95 None
96 }
97 }
98 _ => None,
99 }
100 }
101
102 /// Increments a [`Float`] by its ulp. See [`Float::ulp`] for details.
103 ///
104 /// If the [`Float`] is positive and is the largest [`Float`] in its binade with its precision,
105 /// then
106 /// - If its exponent is not the maximum exponent, it will become the smallest [`Float`] in the
107 /// next-higher binade, and its precision will increase by 1 (so that its ulp remains the
108 /// same);
109 /// - If its exponent is the maximum exponent, it will become $\infty$.
110 ///
111 /// If the [`Float`] is negative and is closer to zero than any other [`Float`] in its binade
112 /// with its precision, then
113 /// - If its precision is 1, it will become negative zero.
114 /// - If its precision is greater than 1 and its exponent is not the minimum exponent, it will
115 /// become the farthest-from-zero [`Float`] in the next-lower binade, and its precision will
116 /// decrease by 1 (so that its ulp remains the same).
117 /// - If its precision is greater than 1 and its exponent is the minimum exponent, it will
118 /// become negative zero.
119 ///
120 /// # Worst-case complexity
121 /// $T(n) = O(n)$
122 ///
123 /// $M(n) = O(n)$
124 ///
125 /// where $T$ is time, $M$ is additional memory, and $n$ is `self.significant_bits()`.
126 ///
127 /// # Panics
128 /// Panics if `self` is NaN, infinite, or zero.
129 ///
130 /// # Examples
131 /// ```
132 /// use malachite_base::num::arithmetic::traits::PowerOf2;
133 /// use malachite_base::num::basic::traits::{NegativeOne, One};
134 /// use malachite_float::Float;
135 ///
136 /// let mut x = Float::ONE;
137 /// assert_eq!(x.to_string(), "1.0");
138 /// x.increment();
139 /// assert_eq!(x.to_string(), "2.0");
140 ///
141 /// let mut x = Float::one_prec(100);
142 /// assert_eq!(x.to_string(), "1.0");
143 /// x.increment();
144 /// assert_eq!(x.to_string(), "1.000000000000000000000000000002");
145 ///
146 /// let mut x = Float::from(std::f64::consts::PI);
147 /// assert_eq!(x.to_string(), "3.141592653589793");
148 /// x.increment();
149 /// assert_eq!(x.to_string(), "3.141592653589797");
150 ///
151 /// let mut x = Float::power_of_2(100u64);
152 /// assert_eq!(x.to_string(), "1.0e30");
153 /// x.increment();
154 /// assert_eq!(x.to_string(), "3.0e30");
155 ///
156 /// let mut x = Float::power_of_2(-100i64);
157 /// assert_eq!(x.to_string(), "8.0e-31");
158 /// x.increment();
159 /// assert_eq!(x.to_string(), "1.6e-30");
160 ///
161 /// let mut x = Float::NEGATIVE_ONE;
162 /// assert_eq!(x.to_string(), "-1.0");
163 /// x.increment();
164 /// assert_eq!(x.to_string(), "-0.0");
165 /// ```
166 pub fn increment(&mut self) {
167 if self.is_sign_negative() {
168 self.neg_assign();
169 self.decrement();
170 self.neg_assign();
171 } else if let Float(Finite {
172 exponent,
173 precision,
174 significand,
175 ..
176 }) = self
177 {
178 let ulp = Limb::power_of_2(significand_bits(significand) - *precision);
179 let limb_count = significand.limb_count();
180 significand.add_assign_at_limb(
181 usize::wrapping_from(limb_count)
182 - 1
183 - usize::exact_from((*precision - 1) >> Limb::LOG_WIDTH),
184 ulp,
185 );
186 if significand.limb_count() > limb_count {
187 if *exponent == Float::MAX_EXPONENT {
188 *self = Float::INFINITY;
189 return;
190 }
191 *significand >>= 1;
192 *exponent += 1;
193 if precision.divisible_by_power_of_2(Limb::LOG_WIDTH) {
194 *significand <<= Limb::WIDTH;
195 }
196 *precision = precision.checked_add(1).unwrap();
197 }
198 } else {
199 panic!("Cannot increment float is non-finite or zero");
200 }
201 }
202
203 /// Decrements a [`Float`] by its ulp. See [`Float::ulp`] for details.
204 ///
205 /// If the [`Float`] is negative and is the largest [`Float`] in its binade with its precision,
206 /// then
207 /// - If its exponent is not the maximum exponent, it will become the closest-to-zero [`Float`]
208 /// in the next-higher binade, and its precision will increase by 1 (so that its ulp remains
209 /// the same);
210 /// - If its exponent is the maximum exponent, it will become $-\infty$.
211 ///
212 /// If the [`Float`] is positive and is smaller than any other [`Float`] in its binade with its
213 /// precision, then
214 /// - If its precision is 1, it will become positive zero.
215 /// - If its precision is greater than 1 and its exponent is not the minimum exponent, it will
216 /// become the largest [`Float`] in the next-lower binade, and its precision will decrease by
217 /// 1 (so that its ulp remains the same).
218 /// - If its precision is greater than 1 and its exponent is the minimum exponent, it will
219 /// become positive zero.
220 ///
221 /// # Worst-case complexity
222 /// $T(n) = O(n)$
223 ///
224 /// $M(n) = O(n)$
225 ///
226 /// where $T$ is time, $M$ is additional memory, and $n$ is `self.significant_bits()`.
227 ///
228 /// # Panics
229 /// Panics if `self` is NaN, infinite, or zero.
230 ///
231 /// # Examples
232 /// ```
233 /// use malachite_base::num::arithmetic::traits::PowerOf2;
234 /// use malachite_base::num::basic::traits::{NegativeOne, One};
235 /// use malachite_float::Float;
236 ///
237 /// let mut x = Float::ONE;
238 /// assert_eq!(x.to_string(), "1.0");
239 /// x.decrement();
240 /// assert_eq!(x.to_string(), "0.0");
241 ///
242 /// let mut x = Float::one_prec(100);
243 /// assert_eq!(x.to_string(), "1.0");
244 /// x.decrement();
245 /// assert_eq!(x.to_string(), "0.999999999999999999999999999998");
246 ///
247 /// let mut x = Float::from(std::f64::consts::PI);
248 /// assert_eq!(x.to_string(), "3.141592653589793");
249 /// x.decrement();
250 /// assert_eq!(x.to_string(), "3.14159265358979");
251 ///
252 /// let mut x = Float::power_of_2(100u64);
253 /// assert_eq!(x.to_string(), "1.0e30");
254 /// x.decrement();
255 /// assert_eq!(x.to_string(), "0.0");
256 ///
257 /// let mut x = Float::power_of_2(-100i64);
258 /// assert_eq!(x.to_string(), "8.0e-31");
259 /// x.decrement();
260 /// assert_eq!(x.to_string(), "0.0");
261 ///
262 /// let mut x = Float::NEGATIVE_ONE;
263 /// assert_eq!(x.to_string(), "-1.0");
264 /// x.decrement();
265 /// assert_eq!(x.to_string(), "-2.0");
266 /// ```
267 pub fn decrement(&mut self) {
268 if self.is_sign_negative() {
269 self.neg_assign();
270 self.increment();
271 self.neg_assign();
272 } else if let Float(Finite {
273 exponent,
274 precision,
275 significand,
276 ..
277 }) = self
278 {
279 let bits = significand_bits(significand);
280 let ulp = Limb::power_of_2(bits - *precision);
281 significand.sub_assign_at_limb(
282 usize::wrapping_from(significand.limb_count())
283 - 1
284 - usize::exact_from((*precision - 1) >> Limb::LOG_WIDTH),
285 ulp,
286 );
287 if *significand == 0u32 {
288 *self = Float::ZERO;
289 } else if significand.significant_bits() < bits {
290 if *exponent == Float::MIN_EXPONENT {
291 *self = Float::ZERO;
292 return;
293 }
294 *significand <<= 1;
295 *exponent -= 1;
296 *precision = precision.checked_sub(1).unwrap();
297 if bits - *precision == Limb::WIDTH {
298 *significand >>= Limb::WIDTH;
299 }
300 }
301 } else {
302 panic!("Cannot decrement float that is non-finite or zero");
303 }
304 }
305}