1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362
//! Implementations of methods related to simplification.
//!
//! Note that these methods are only implemented for [RBig] but not for [Relaxed][crate::Relaxed]
//! because the latter one is naturally not in the simplest form.
use crate::{error::panic_divide_by_0, rbig::RBig, repr::Repr};
use core::{cmp::Ordering, mem};
use dashu_base::{AbsOrd, Approximation, DivRem, UnsignedAbs};
use dashu_int::{IBig, Sign, UBig};
impl Repr {
/// Find the simplest rational number in the open interval `(lower, upper)`.
/// See [RBig::simplest_in()] and <https://stackoverflow.com/q/66980340/5960776>.
pub fn simplest_in(mut lower: Self, mut upper: Self) -> Self {
let sign = if lower.numerator.sign() != upper.numerator.sign() {
// if lower < 0 < upper, then 0 is the simplest
return Self::zero();
} else {
lower.numerator.sign()
};
lower = lower.abs();
upper = upper.abs();
match lower.cmp(&upper) {
// swap so that lower is less than upper
Ordering::Greater => mem::swap(&mut lower, &mut upper),
Ordering::Equal => return sign * lower,
Ordering::Less => {}
}
let Repr {
numerator: mut num_l,
denominator: den_l,
} = lower;
let Repr {
numerator: mut num_r,
denominator: den_r,
} = upper;
// negative values might exist during the calculation
let (mut den_l, mut den_r) = (IBig::from(den_l), IBig::from(den_r));
// use continued fraction expansion to find this float
let (mut n0, mut d0) = (IBig::ONE, IBig::ZERO);
let (mut n1, mut d1) = (IBig::ZERO, IBig::ONE);
let (num, den) = loop {
let (q, r1) = num_l.div_rem(&den_l);
n1 += &q * &n0;
mem::swap(&mut n0, &mut n1);
d1 += &q * &d0;
mem::swap(&mut d0, &mut d1);
let r2 = mem::take(&mut num_r) - q * &den_r;
num_l = mem::replace(&mut den_r, r1);
num_r = mem::replace(&mut den_l, r2);
if num_l < den_l {
break (n0 + n1, d0 + d1);
}
};
debug_assert!(num.sign() == den.sign());
Repr {
numerator: num.unsigned_abs() * sign,
denominator: den.unsigned_abs(),
}
}
}
/// Implementation of simplest_from_f32, simplest_from_f64
macro_rules! impl_simplest_from_float {
($f:ident) => {{
if $f.is_infinite() || $f.is_nan() {
return None;
} else if $f == 0. {
return Some(Self::ZERO);
}
// get the range (f - ulp/2, f + ulp/2)
// if f is negative, then range will be flipped by simplest_in()
let mut est = Repr::try_from($f).unwrap();
est.numerator <<= 1;
est.denominator <<= 1;
let left = Self(
Repr {
numerator: &est.numerator + IBig::ONE,
denominator: est.denominator.clone(),
}
.reduce(),
);
let right = Self(
Repr {
numerator: est.numerator - IBig::ONE,
denominator: est.denominator,
}
.reduce(),
);
// find the simplest float in the range
let mut simplest = Self::simplest_in(left.clone(), right.clone());
if $f.to_bits() & 1 == 0 {
// consider boundry values when last bit is 0 (because ties to even)
if left.is_simpler_than(&simplest) {
simplest = left;
}
if right.is_simpler_than(&simplest) {
simplest = right;
}
}
Some(simplest)
}};
}
impl RBig {
/// Determine if this rational number is simpler than the other number.
///
/// This method only make sense for canonicalized ratios.
#[inline]
pub fn is_simpler_than(&self, other: &Self) -> bool {
(self.denominator() < other.denominator()) // first compare denominator
&& self.numerator().abs_cmp(other.numerator()).is_le() // then compare numerator
&& self.sign() > other.sign() // then compare sign
}
/// Find the simplest rational number in the rounding interval of the [f32] number.
///
/// This method returns None when the floating point value is not representable by a rational number,
/// such as infinities or nans.
///
/// See [RBig::simplest_in] for the definition of `simplicity`.
///
/// The rounding interval of a [f32] value is an interval such that all numbers in this
/// range will rounded to this [f32] value. For example the rounding interval for `1f32`
/// is `[1. - f32::EPSILON / 2, 1. + f32::EPSILON / 2]`. That is, the error of result value will
/// be less than 1/2 ULP.
///
/// This method can be used to recover the original fraction represented as a division of [f32].
/// See the examples below.
///
/// # Examples
///
/// ```
/// # use dashu_ratio::RBig;
/// assert_eq!(
/// RBig::simplest_from_f32(1e-1).unwrap(),
/// RBig::from_parts(1.into(), 10u8.into())
/// );
/// assert_eq!(
/// RBig::simplest_from_f32(22./7.).unwrap(),
/// RBig::from_parts(22.into(), 7u8.into())
/// );
/// ```
pub fn simplest_from_f32(f: f32) -> Option<Self> {
impl_simplest_from_float!(f)
}
/// Find the simplest rational number in the rounding interval of the [f64] number.
///
/// This method returns None when the floating point value is not representable by a rational number,
/// such as infinities or nans.
///
/// See [RBig::simplest_in] for the definition of `simplicity`.
///
/// The rounding interval of a [f64] value is an interval such that all numbers in this
/// range will rounded to this [f64] value. For example the rounding interval for `1f64`
/// is `[1. - f64::EPSILON / 2, 1. + f64::EPSILON / 2]`. That is, the error of result value will
/// be less than 1/2 ULP.
///
/// This method can be used to recover the original fraction represented as a division of [f64].
/// See the examples below.
///
/// # Examples
///
/// ```
/// # use dashu_ratio::RBig;
/// assert_eq!(
/// RBig::simplest_from_f64(1e-1).unwrap(),
/// RBig::from_parts(1.into(), 10u8.into())
/// );
/// assert_eq!(
/// RBig::simplest_from_f64(22./7.).unwrap(),
/// RBig::from_parts(22.into(), 7u8.into())
/// );
pub fn simplest_from_f64(f: f64) -> Option<Self> {
impl_simplest_from_float!(f)
}
/// Find the simplest rational number in the open interval `(lower, upper)`.
///
/// A rational `n₁/d₁` is simpler than another rational number `n₂/d₂` if:
/// * `d₁ < d₂` (compare denominator)
/// * or `|n₁| < |n₂|` (then compare the magnitude of numerator)
/// * or `n₂ < 0 < n₁` (then compare the sign)
///
/// `lower` and `upper` will be swapped if necessary. If `lower` and `upper` are
/// the same number, then this number will be directly returned.
///
/// # Examples
///
/// ```
/// # use dashu_ratio::RBig;
/// let a = RBig::from_parts(1234.into(), 5678u16.into());
/// let b = RBig::from_parts(1235.into(), 5679u16.into());
/// let s = RBig::simplest_in(a, b);
/// // 1234/5678 < 5/23 < 1235/5679
/// assert_eq!(s, RBig::from_parts(5.into(), 23u8.into()));
/// ```
///
#[inline]
pub fn simplest_in(lower: Self, upper: Self) -> Self {
Self(Repr::simplest_in(lower.0, upper.0).reduce())
}
/// Find the previous and next value in farey sequence (from 0 to 1) with `limit` as the order.
///
/// This function requires `-1 < x < 1` and `x.denominator` > `limit`
fn farey_neighbors(x: &Self, limit: &UBig) -> (Self, Self) {
debug_assert!(x.denominator() > limit);
debug_assert!(!x.numerator().is_zero());
debug_assert!(x.numerator().abs_cmp(x.denominator()).is_le());
let (mut left, mut right) = match x.sign() {
Sign::Positive => (Repr::zero(), Repr::one()),
Sign::Negative => (Repr::neg_one(), Repr::zero()),
};
// the Farey neighbors can be found directly by adding the
// numerator and denominator together, see <https://en.wikipedia.org/wiki/Farey_sequence#Farey_neighbours>.
loop {
let mut next = Repr {
numerator: &left.numerator + &right.numerator,
denominator: &left.denominator + &right.denominator,
};
// test if the denominator has exceeded the limit
if &next.denominator > limit {
next = next.reduce();
if &next.denominator > limit {
return (Self(left), Self(right));
}
}
// tighten the bounds
if next > x.0 {
right = next;
} else {
left = next;
}
}
}
/// Find the closest rational number to this number with a limit of the denominator.
///
/// If the denominator of this number is larger than the limit, then it returns the closest one
/// between `self.next_up()` and `self.next_down()` to `self`. If the denominator of this number
/// is already less than or equal to the limit, then `Exact(self)` will be returned.
///
/// The error `|self - self.nearest()|` will be less than `1/(2*limit)`, and the sign of
/// the error `self - self.nearest()` will be returned if the result is not [Exact][Approximation::Exact].
///
/// # Examples
///
/// ```
/// # use dashu_base::{Approximation::*, Sign};
/// # use dashu_ratio::RBig;
/// let a: RBig = 3.141592653.try_into().unwrap();
/// assert_eq!(a.nearest(&10u8.into()), Inexact(
/// RBig::from_parts(22.into(), 7u8.into()),
/// Sign::Positive // 22/7 > 3.141592653
/// ));
/// ```
pub fn nearest(&self, limit: &UBig) -> Approximation<Self, Sign> {
if limit.is_zero() {
panic_divide_by_0()
}
// directly return this number if it's already simple enough
if self.denominator() <= limit {
return Approximation::Exact(self.clone());
}
let (trunc, r) = self.clone().split_at_point();
let (left, right) = Self::farey_neighbors(&r, limit);
// find the closest one (compare r - left and right - r)
let mut mid = (&left + &right).0;
mid.denominator <<= 1;
if r.0 > mid {
Approximation::Inexact(trunc + right, Sign::Positive)
} else {
Approximation::Inexact(trunc + left, Sign::Negative)
}
}
/// Find the closest rational number that is greater than this number and has a denominator less than `limit`.
///
/// It's equivalent to finding the next element in Farey sequence of order `limit`. The error
/// `|self - self.next_up()|` will be less than `1/limit`.
///
/// # Examples
///
/// ```
/// # use dashu_ratio::RBig;
/// let a: RBig = 3.141592653.try_into().unwrap();
/// assert_eq!(a.next_up(&10u8.into()), RBig::from_parts(22.into(), 7u8.into()));
/// ```
pub fn next_up(&self, limit: &UBig) -> Self {
if limit.is_zero() {
panic_divide_by_0()
}
let (trunc, fract) = self.clone().split_at_point();
let up = if self.denominator() <= limit {
// If the denominator of the number is already small enough, increase the number a little
// bit before finding the farey neighbors. Note that the distance between two adjacent
// numbers in a farey sequence is at least limit^-2, so we just increase limit^-2
let target = fract
+ Self(Repr {
numerator: IBig::ONE,
denominator: limit.sqr(),
});
Self::farey_neighbors(&target, limit).1
} else {
// otherwise we can directly find the next value by finding farey bounds
Self::farey_neighbors(&fract, limit).1
};
trunc + up
}
/// Find the closest rational number that is less than this number and has a denominator less than `limit`.
///
/// It's equivalent to finding the previous element in Farey sequence of order `limit`. The error
/// `|self - self.next_down()|` will be less than `1/limit`.
///
/// # Examples
///
/// ```
/// # use dashu_ratio::RBig;
/// let a: RBig = 3.141592653.try_into().unwrap();
/// assert_eq!(a.next_down(&10u8.into()), RBig::from_parts(25.into(), 8u8.into()));
/// ```
pub fn next_down(&self, limit: &UBig) -> Self {
if limit.is_zero() {
panic_divide_by_0()
}
// similar to next_up()
let (trunc, fract) = self.clone().split_at_point();
let down = if self.denominator() <= limit {
let target = fract
- Self(Repr {
numerator: IBig::ONE,
denominator: limit.sqr(),
});
Self::farey_neighbors(&target, limit).0
} else {
Self::farey_neighbors(&fract, limit).0
};
trunc + down
}
}