zenith_float_num/
rational.rs

1//! Exact rationals `num/den` with `ExactNum` limbs, reduced to lowest terms.
2
3use crate::defs::RoundingMode;
4use crate::defs::DEFAULT_P;
5use crate::defs::WORD_BIT_SIZE;
6use crate::ExactInt;
7use crate::ExactNum;
8use crate::Radix;
9use crate::NAN;
10use alloc::string::String;
11use alloc::vec::Vec;
12use core::cmp::Ordering;
13
14/// Exact rational `num/den`. Distinct from [`ExactNum`], which is floating-point.
15///
16/// Both parts are integer-valued `ExactNum`s after [`Self::new`]. The
17/// denominator is positive. A zero denominator or a non-finite part is stored
18/// as `NaN/1`.
19#[derive(Clone, Debug)]
20pub struct ExactRational {
21    num: ExactNum,
22    den: ExactNum,
23}
24
25fn rat_one() -> ExactNum {
26    ExactNum::from_u8(1, DEFAULT_P)
27}
28
29fn work_p(xs: &[&ExactNum]) -> usize {
30    let mut p = DEFAULT_P;
31    for x in xs {
32        if let Some(px) = x.precision() {
33            p = p.max(px);
34        }
35    }
36    p.saturating_add(WORD_BIT_SIZE)
37}
38
39fn is_int(x: &ExactNum) -> bool {
40    if x.is_nan() || x.is_inf() {
41        return false;
42    }
43    x.fract().is_zero()
44}
45
46fn gcd_int(mut a: ExactNum, mut b: ExactNum) -> ExactNum {
47    a = a.abs();
48    b = b.abs();
49    while !b.is_zero() {
50        let r = a.rem(&b);
51        a = b;
52        b = r;
53    }
54    if a.is_zero() {
55        rat_one()
56    } else {
57        a
58    }
59}
60
61fn reduce(mut num: ExactNum, mut den: ExactNum) -> ExactRational {
62    if num.is_nan() || den.is_nan() || num.is_inf() || den.is_inf() || den.is_zero() {
63        return ExactRational {
64            num: NAN,
65            den: rat_one(),
66        };
67    }
68    if den.is_negative() {
69        num = num.neg();
70        den = den.neg();
71    }
72    if is_int(&num) && is_int(&den) {
73        let g = gcd_int(num.clone(), den.clone());
74        if !g.is_zero() && g.cmp(&rat_one()) != Some(0) {
75            let p = work_p(&[&num, &den, &g]);
76            num = num.div(&g, p, RoundingMode::ToEven);
77            den = den.div(&g, p, RoundingMode::ToEven);
78        }
79    }
80    ExactRational { num, den }
81}
82
83impl ExactRational {
84    /// `num/den` reduced to lowest terms. Negative `den` moves the sign to `num`.
85    /// Zero `den` or a non-finite part is `NaN`.
86    pub fn new(num: ExactNum, den: ExactNum) -> Self {
87        reduce(num, den)
88    }
89
90    /// Integer ratio `n/d` at default limb precision.
91    pub fn from_i64(n: i64, d: i64) -> Self {
92        Self::new(
93            ExactNum::from_i64(n, DEFAULT_P),
94            ExactNum::from_i64(d, DEFAULT_P),
95        )
96    }
97
98    /// Numerator (integer-valued after construction).
99    pub fn num(&self) -> &ExactNum {
100        &self.num
101    }
102
103    /// Denominator (positive integer-valued after construction).
104    pub fn den(&self) -> &ExactNum {
105        &self.den
106    }
107
108    /// True if a part is NaN or the denominator was zero.
109    pub fn is_nan(&self) -> bool {
110        self.num.is_nan() || self.den.is_nan()
111    }
112
113    /// True if the denominator is `1`.
114    pub fn is_integer(&self) -> bool {
115        !self.is_nan() && self.den.cmp(&rat_one()) == Some(0)
116    }
117
118    /// `self + rhs` as an exact rational.
119    pub fn add(&self, rhs: &Self) -> Self {
120        if self.is_nan() || rhs.is_nan() {
121            return Self::new(NAN, rat_one());
122        }
123        let p = work_p(&[&self.num, &self.den, &rhs.num, &rhs.den]);
124        let rm = RoundingMode::None;
125        let ad = self.num.mul(&rhs.den, p, rm);
126        let bc = rhs.num.mul(&self.den, p, rm);
127        let num = ad.add(&bc, p, rm);
128        let den = self.den.mul(&rhs.den, p, rm);
129        Self::new(num, den)
130    }
131
132    /// `self - rhs`.
133    pub fn sub(&self, rhs: &Self) -> Self {
134        self.add(&Self::new(rhs.num.neg(), rhs.den.clone()))
135    }
136
137    /// `self * rhs`.
138    pub fn mul(&self, rhs: &Self) -> Self {
139        if self.is_nan() || rhs.is_nan() {
140            return Self::new(NAN, rat_one());
141        }
142        let p = work_p(&[&self.num, &self.den, &rhs.num, &rhs.den]);
143        let rm = RoundingMode::None;
144        Self::new(self.num.mul(&rhs.num, p, rm), self.den.mul(&rhs.den, p, rm))
145    }
146
147    /// `self / rhs`. Zero `rhs` is `NaN`.
148    pub fn div(&self, rhs: &Self) -> Self {
149        if self.is_nan() || rhs.is_nan() || rhs.num.is_zero() {
150            return Self::new(NAN, rat_one());
151        }
152        let p = work_p(&[&self.num, &self.den, &rhs.num, &rhs.den]);
153        let rm = RoundingMode::None;
154        Self::new(self.num.mul(&rhs.den, p, rm), self.den.mul(&rhs.num, p, rm))
155    }
156
157    /// Convert to an `ExactNum` at `(p, rm)`.
158    pub fn to_exact_num(&self, p: usize, rm: RoundingMode) -> ExactNum {
159        if self.is_nan() {
160            return NAN;
161        }
162        self.num.div(&self.den, p, rm)
163    }
164
165    fn quot(&self) -> ExactNum {
166        let p = work_p(&[&self.num, &self.den]);
167        self.num.div(&self.den, p, RoundingMode::None)
168    }
169
170    /// Greatest integer `≤ self` as a rational with denominator `1`.
171    pub fn floor(&self) -> Self {
172        if self.is_nan() {
173            return Self::new(NAN, rat_one());
174        }
175        Self::new(self.quot().floor(), rat_one())
176    }
177
178    /// Least integer `≥ self` as a rational with denominator `1`.
179    pub fn ceil(&self) -> Self {
180        if self.is_nan() {
181            return Self::new(NAN, rat_one());
182        }
183        Self::new(self.quot().ceil(), rat_one())
184    }
185
186    /// Nearest integer, ties away from zero, as a rational with denominator `1`.
187    pub fn round(&self) -> Self {
188        if self.is_nan() {
189            return Self::new(NAN, rat_one());
190        }
191        Self::new(self.quot().round(0, RoundingMode::FromZero), rat_one())
192    }
193
194    /// Build `num/den` from limb integers.
195    pub fn from_ints(num: ExactInt, den: ExactInt) -> Self {
196        if den.is_zero() {
197            return Self::new(NAN, rat_one());
198        }
199        let p = num
200            .bit_length()
201            .max(den.bit_length())
202            .max(DEFAULT_P)
203            .saturating_add(WORD_BIT_SIZE);
204        Self::new(
205            num.to_exact_num(p, RoundingMode::None),
206            den.to_exact_num(p, RoundingMode::None),
207        )
208    }
209
210    /// Parse a decimal (optional `e`/`E` exponent) as an exact rational.
211    ///
212    /// `0.1` is `1/10`, not a binary float. Invalid syntax is `None`.
213    pub fn parse_exact(s: &str) -> Option<Self> {
214        let s = s.trim();
215        if s.is_empty() {
216            return None;
217        }
218        let (neg, rest) = match s.as_bytes()[0] {
219            b'+' => (false, &s[1..]),
220            b'-' => (true, &s[1..]),
221            _ => (false, s),
222        };
223        if rest.is_empty() {
224            return None;
225        }
226        let (mant, exp) = split_dec_exp(rest)?;
227        let (int_part, frac_part) = split_dot(mant);
228        if int_part.is_empty() && frac_part.is_empty() {
229            return None;
230        }
231        if !int_part.bytes().all(|c| c.is_ascii_digit())
232            || !frac_part.bytes().all(|c| c.is_ascii_digit())
233        {
234            return None;
235        }
236        let mut digits = String::new();
237        digits.push_str(int_part);
238        digits.push_str(frac_part);
239        if digits.is_empty() || digits.bytes().all(|c| c == b'0') {
240            return Some(Self::from_i64(0, 1));
241        }
242        let num = dec_digits_to_int(&digits)?;
243        let exp_adj = exp - frac_part.len() as i32;
244        let ten = ExactInt::from_i64(10);
245        let (n, d) = if exp_adj >= 0 {
246            (num.mul(&ten.pow(exp_adj as u64)), ExactInt::one())
247        } else {
248            (num, ten.pow((-exp_adj) as u64))
249        };
250        let mut r = Self::from_ints(n, d);
251        if neg {
252            r = Self::new(r.num.neg(), r.den.clone());
253        }
254        Some(r)
255    }
256
257    /// Minimum digits in `rdx` for an exact terminating representation. `None` if the
258    /// denominator has a prime factor that does not divide the radix.
259    pub fn format_exact(&self, rdx: Radix) -> Option<String> {
260        if self.is_nan() {
261            return None;
262        }
263        let num = ExactInt::from_exact_num(&self.num)?;
264        let den = ExactInt::from_exact_num(&self.den)?;
265        format_terminating(&num, &den, rdx)
266    }
267
268    /// Exact comparison via cross-multiplication. `None` if either value is NaN.
269    pub fn partial_cmp(&self, other: &Self) -> Option<Ordering> {
270        if self.is_nan() || other.is_nan() {
271            return None;
272        }
273        let p = work_p(&[&self.num, &self.den, &other.num, &other.den]);
274        let rm = RoundingMode::None;
275        let left = self.num.mul(&other.den, p, rm);
276        let right = other.num.mul(&self.den, p, rm);
277        match left.cmp(&right) {
278            Some(c) if c > 0 => Some(Ordering::Greater),
279            Some(c) if c < 0 => Some(Ordering::Less),
280            Some(_) => Some(Ordering::Equal),
281            None => None,
282        }
283    }
284}
285
286fn split_dec_exp(s: &str) -> Option<(&str, i32)> {
287    let bytes = s.as_bytes();
288    let mut epos = None;
289    for (i, c) in bytes.iter().enumerate() {
290        if *c == b'e' || *c == b'E' {
291            epos = Some(i);
292            break;
293        }
294    }
295    match epos {
296        None => Some((s, 0)),
297        Some(i) => {
298            if i == 0 {
299                return None;
300            }
301            let exp_s = &s[i + 1..];
302            if exp_s.is_empty() {
303                return None;
304            }
305            let exp: i32 = exp_s.parse().ok()?;
306            Some((&s[..i], exp))
307        }
308    }
309}
310
311fn split_dot(s: &str) -> (&str, &str) {
312    match s.find('.') {
313        Some(i) => (&s[..i], &s[i + 1..]),
314        None => (s, ""),
315    }
316}
317
318fn dec_digits_to_int(s: &str) -> Option<ExactInt> {
319    let ten = ExactInt::from_i64(10);
320    let mut v = ExactInt::zero();
321    for c in s.bytes() {
322        if !c.is_ascii_digit() {
323            return None;
324        }
325        v = v.mul(&ten).add(&ExactInt::from_u64((c - b'0') as u64));
326    }
327    Some(v)
328}
329
330fn format_terminating(num: &ExactInt, den: &ExactInt, rdx: Radix) -> Option<String> {
331    if den.is_zero() {
332        return None;
333    }
334    let radix = ExactInt::from_u64(rdx.value() as u64);
335    let mut pk = ExactInt::one();
336    let max_k = den.bit_length().saturating_add(8);
337    for k in 0..=max_k {
338        if let Some((scale, rem)) = pk.div_rem(den) {
339            if rem.is_zero() {
340                let digits = num.mul(&scale);
341                return Some(format_fixed(&digits, k, rdx));
342            }
343        }
344        pk = pk.mul(&radix);
345    }
346    None
347}
348
349fn format_fixed(n: &ExactInt, frac_digits: usize, rdx: Radix) -> String {
350    let neg = n.is_negative();
351    let abs = if neg { n.neg() } else { n.clone() };
352    let mut digits = to_radix_digits(&abs, rdx);
353    let mut k = frac_digits;
354    while k > 0 && digits.ends_with('0') {
355        digits.pop();
356        k -= 1;
357    }
358    if digits.is_empty() {
359        digits.push('0');
360    }
361    let mut out = String::new();
362    if neg && digits != "0" {
363        out.push('-');
364    }
365    if k == 0 {
366        out.push_str(&digits);
367        return out;
368    }
369    if digits.len() <= k {
370        out.push('0');
371        out.push('.');
372        for _ in 0..(k - digits.len()) {
373            out.push('0');
374        }
375        out.push_str(&digits);
376    } else {
377        let split = digits.len() - k;
378        out.push_str(&digits[..split]);
379        out.push('.');
380        out.push_str(&digits[split..]);
381    }
382    out
383}
384
385fn to_radix_digits(n: &ExactInt, rdx: Radix) -> String {
386    if n.is_zero() {
387        return String::from("0");
388    }
389    let r = ExactInt::from_u64(rdx.value() as u64);
390    let mut v = if n.is_negative() { n.neg() } else { n.clone() };
391    let mut digits = Vec::new();
392    while !v.is_zero() {
393        let (q, rem) = v
394            .div_rem(&r)
395            .unwrap_or((ExactInt::zero(), ExactInt::zero()));
396        let d = rem.low_word() as u8;
397        digits.push(if d < 10 { b'0' + d } else { b'a' + (d - 10) });
398        v = q;
399    }
400    digits.reverse();
401    String::from_utf8(digits).unwrap_or_else(|_| String::from("0"))
402}
403
404fn exact_num_to_rational(x: &ExactNum) -> Option<ExactRational> {
405    if x.is_nan() || x.is_inf() {
406        return None;
407    }
408    if x.is_zero() {
409        return Some(ExactRational::from_i64(0, 1));
410    }
411    let (m, _n, s, e, _) = x.as_raw_parts()?;
412    let pbuf = m.len() * WORD_BIT_SIZE;
413    let mut num = ExactInt::from_le_words(s, m);
414    let shift = e as i64 - pbuf as i64;
415    if shift >= 0 {
416        num = num.shl(shift as usize);
417        Some(ExactRational::from_ints(num, ExactInt::one()))
418    } else {
419        let den = ExactInt::one().shl((-shift) as usize);
420        Some(ExactRational::from_ints(num, den))
421    }
422}
423
424impl ExactNum {
425    /// Parse a decimal that is exact in binary (a dyadic rational). `None` if the
426    /// value is not a finite dyadic (so `0.1` is `None` here; use [`ExactRational::parse_exact`]).
427    pub fn parse_exact(s: &str) -> Option<Self> {
428        let r = ExactRational::parse_exact(s)?;
429        let den = ExactInt::from_exact_num(r.den())?;
430        if !den.is_one() && !is_pow2(&den) {
431            return None;
432        }
433        let p = den.bit_length().max(WORD_BIT_SIZE);
434        Some(r.to_exact_num(p, RoundingMode::None))
435    }
436
437    /// Minimum digits in `rdx` that recover `self` exactly when the value is a
438    /// terminating expansion in that radix.
439    pub fn format_exact(&self, rdx: Radix) -> Option<String> {
440        exact_num_to_rational(self)?.format_exact(rdx)
441    }
442}
443
444fn is_pow2(n: &ExactInt) -> bool {
445    if n.is_zero() || n.is_negative() {
446        return false;
447    }
448    n.sub(&ExactInt::one()).bit_length() < n.bit_length()
449}
450
451impl PartialEq for ExactRational {
452    fn eq(&self, other: &Self) -> bool {
453        matches!(self.partial_cmp(other), Some(Ordering::Equal))
454    }
455}
456
457#[cfg(test)]
458mod tests {
459    use super::*;
460    use crate::Consts;
461    use crate::Radix;
462
463    fn r(n: i64, d: i64) -> ExactRational {
464        ExactRational::new(
465            ExactNum::from_i64(n, DEFAULT_P),
466            ExactNum::from_i64(d, DEFAULT_P),
467        )
468    }
469
470    #[test]
471    fn rational_reduce_add_sign_to_float() {
472        let half = r(1, 2);
473        let two_four = r(2, 4);
474        assert_eq!(two_four, half);
475        assert_eq!(ExactRational::from_i64(2, 4), ExactRational::from_i64(1, 2));
476        assert!(!two_four.is_integer());
477        assert_eq!(
478            two_four.num().cmp(&ExactNum::from_u8(1, DEFAULT_P)),
479            Some(0)
480        );
481        assert_eq!(
482            two_four.den().cmp(&ExactNum::from_u8(2, DEFAULT_P)),
483            Some(0)
484        );
485
486        let s = r(1, 3).add(&r(1, 6));
487        assert_eq!(s, half);
488
489        let pos = ExactRational::from_i64(-3, -4);
490        let want = ExactRational::from_i64(3, 4);
491        assert_eq!(pos, want);
492        assert!(pos.den().is_positive());
493
494        let p = 256;
495        let rm = RoundingMode::ToEven;
496        let got = ExactRational::from_i64(1, 3).to_exact_num(p, rm);
497        let want_f = ExactNum::from_i64(1, p).div(&ExactNum::from_i64(3, p), p, rm);
498        assert_eq!(got.cmp(&want_f), Some(0));
499        let mut cc = Consts::new().unwrap();
500        let parsed = ExactNum::parse(
501            "0.33333333333333333333333333333333333333333333333333333333333333333333333333333333",
502            Radix::Dec,
503            p,
504            rm,
505            &mut cc,
506        );
507        assert_eq!(got.cmp(&parsed), Some(0));
508
509        assert!(ExactRational::from_i64(1, 3).partial_cmp(&half) == Some(Ordering::Less));
510        assert!(ExactRational::from_i64(2, 2).is_integer());
511        assert_eq!(
512            ExactRational::from_i64(5, 3).floor(),
513            ExactRational::from_i64(1, 1)
514        );
515        assert_eq!(
516            ExactRational::from_i64(5, 3).ceil(),
517            ExactRational::from_i64(2, 1)
518        );
519        assert!(ExactRational::from_i64(1, 0).is_nan());
520    }
521
522    #[test]
523    fn parse_exact_tenth_half_eighth_sci() {
524        let tenth = ExactRational::parse_exact("0.1").unwrap();
525        assert_eq!(
526            tenth.mul(&ExactRational::from_i64(10, 1)),
527            ExactRational::from_i64(1, 1)
528        );
529
530        let half = ExactNum::parse_exact("0.5").unwrap();
531        let p = half.precision().unwrap();
532        let want = ExactNum::from_u8(1, p).div(&ExactNum::from_u8(2, p), p, RoundingMode::ToEven);
533        assert_eq!(half.cmp(&want), Some(0));
534        assert_eq!(half.ilogb(), Some(-1));
535        let (_m, _n, s, _e, _) = half.as_raw_parts().unwrap();
536        assert_eq!(s, crate::Sign::Pos);
537
538        assert_eq!(
539            ExactRational::parse_exact("0.125")
540                .unwrap()
541                .format_exact(Radix::Dec)
542                .as_deref(),
543            Some("0.125")
544        );
545        assert_eq!(
546            ExactNum::parse_exact("0.125")
547                .unwrap()
548                .format_exact(Radix::Dec)
549                .as_deref(),
550            Some("0.125")
551        );
552        assert!(ExactNum::parse_exact("0.1").is_none());
553
554        let mut cc = Consts::new().unwrap();
555        let sci = ExactNum::parse("1.5e3", Radix::Dec, 256, RoundingMode::ToEven, &mut cc);
556        assert_eq!(sci.cmp(&ExactNum::from_i32(1500, 256)), Some(0));
557    }
558}
zenith_float_num/rational.rs

zenith_float_num/
rational.rs
