rustymath_core/
exact.rs

1//! Exact numbers: rationals of any integer type, such as `BigRational`, which is `Ratio<BigInt>`. These functions give exact results, where the float functions round. [`linalg::exact`](crate::linalg::exact) solves matrices of them, [`polynomial`](crate::polynomial) evaluates polynomials with them, and [`number_theory`](crate::number_theory) works on whole numbers.
2use crate::{BETWEEN_CHECKS, error::Error, interrupt::check};
3use num_integer::Integer;
4use num_rational::Ratio;
5use num_traits::{One, Signed, ToPrimitive, Zero};
6
7/// `y` modulo `x` with the sign of `x`, `y - x⌊y/x⌋`, which APL writes `x|y`. A zero `x` gives `y`.
8pub fn residue<T: Clone + Integer>(x: &Ratio<T>, y: &Ratio<T>) -> Ratio<T> { if x.is_zero() { y.clone() } else { y - x * (y / x).floor() } }
9
10pub(crate) fn boolean(value: bool) -> num_rational::BigRational { Ratio::from_integer(i64::from(value).into()) }
11
12pub fn divide(x: &num_rational::BigRational, y: &num_rational::BigRational, mathish: bool) -> Result<Option<num_rational::BigRational>, Error> {
13    if y.is_zero() { return Ok((mathish && x.is_zero()).then(num_rational::BigRational::zero)); }
14    Ok(Some(x / y))
15}
16
17pub fn floor_divide(x: &num_rational::BigRational, y: &num_rational::BigRational, mathish: bool) -> Result<Option<num_rational::BigRational>, Error> {
18    Ok(divide(x, y, mathish)?.map(|value| value.floor()))
19}
20
21/// The greatest common divisor of `x` and `y`: the largest rational that divides each a whole number of times. It is the gcd of their numerators over the lcm of their denominators. It is never negative, and the gcd of 0 and 0 is 0.
22pub fn gcd<T: Clone + Integer>(x: &Ratio<T>, y: &Ratio<T>) -> Ratio<T> { Ratio::new(x.numer().gcd(y.numer()), x.denom().lcm(y.denom())) }
23
24/// The least common multiple of `x` and `y`, `x y / gcd(x, y)`, with the sign of `x y`. It is 0 where either is 0.
25pub fn lcm<T: Clone + Integer>(x: &Ratio<T>, y: &Ratio<T>) -> Ratio<T> { let g = gcd(x, y); if g.is_zero() { g } else { x / &g * y } }
26
27/// The binomial coefficient C(n, k) for a rational `n` and a whole `k`, with J's extension to negative `k`. For a whole `n` from 0, a `k` above `n` gives 0. A negative `k` gives 0, unless `n` is a whole number with `k ≤ n < 0`, where the result is (-1)^(n-k) C(-k-1, n-k). Otherwise the result is the product of (n-i)/(i+1) for each `i` from 0 below `k`. It checks for an interrupt once every 1,024 steps of that product.
28pub fn binom<T>(n: &Ratio<T>, k: i64) -> Result<Ratio<T>, Error> where T: Clone + Integer + Signed + ToPrimitive + From<i64> {
29    let whole = n.is_integer().then(|| n.to_integer());
30    if k < 0 {
31        let Some(n) = whole else { return Ok(Ratio::zero()) };
32        if !n.is_negative() || T::from(k) > n { return Ok(Ratio::zero()); }
33        let count = (n - T::from(k)).to_i64().expect("n - k lies below -k");
34        let c = binom(&Ratio::from_integer(T::from(-(k + 1))), count)?;
35        return Ok(if count % 2 == 0 { c } else { -c });
36    }
37    let k = match whole {
38        Some(n) if !n.is_negative() => {
39            if T::from(k) > n { return Ok(Ratio::zero()); }
40            let rest = n - T::from(k);
41            if rest < T::from(k) { rest.to_i64().expect("n - k lies below k") } else { k }
42        }
43        _ => k,
44    };
45    let mut value = Ratio::one();
46    for i in 0..k {
47        if i % BETWEEN_CHECKS as i64 == 0 { check()?; }
48        value = value * (n - Ratio::from_integer(T::from(i))) / Ratio::from_integer(T::from(i + 1));
49    }
50    Ok(value)
51}
52
53pub fn factorial(value: &num_rational::BigRational) -> Result<Option<num_rational::BigRational>, Error> {
54    if !value.is_integer() || value.is_negative() { return Ok(None); }
55    let n = value.to_u64().filter(|&n| n <= 100_000).ok_or_else(|| Error::new(crate::error::Kind::Limit, "exact factorial argument is too large"))?;
56    Ok(Some(Ratio::from_integer(crate::number_theory::factorial(n)?.into())))
57}
58
59pub fn power(x: &num_rational::BigRational, y: &num_rational::BigRational) -> Result<Option<num_rational::BigRational>, Error> {
60    if !y.is_integer() { return Ok(None); }
61    let n = y.to_i32().filter(|&n| n.unsigned_abs() <= 1_000_000).ok_or_else(|| Error::new(crate::error::Kind::Limit, "exact exponent is too large"))?;
62    if x.is_zero() && n < 0 { return Ok(None); }
63    Ok(Some(x.pow(n)))
64}
rustymath_core/exact.rs

rustymath_core/
exact.rs
