improvise/crates/numtower/src/ratio.rs

//! Exact rationals in bigints: the widened intermediate of question 13 (b), and
//! the one correctly rounded conversion of a rational to f64 (clauses 4 and 5).

use std::cmp::Ordering;

use num_bigint::{BigInt, BigUint, Sign};
use num_integer::Integer;
use num_traits::{One, Signed, ToPrimitive, Zero};

/// A reduced fraction with a positive denominator.
#[derive(Clone, Debug, PartialEq, Eq)]
pub(crate) struct BigQ {
    pub(crate) num: BigInt,
    pub(crate) den: BigInt,
}

impl BigQ {
    pub(crate) fn new(num: BigInt, den: BigInt) -> BigQ {
        let (num, den) = match den.is_negative() {
            true => (-num, -den),
            false => (num, den),
        };
        let g = num.gcd(&den);
        match g.is_zero() || g.is_one() {
            true => BigQ { num, den },
            false => BigQ {
                num: num / &g,
                den: den / &g,
            },
        }
    }

    pub(crate) fn from_i128(num: i128, den: i128) -> BigQ {
        BigQ::new(BigInt::from(num), BigInt::from(den))
    }

    pub(crate) fn add(&self, o: &BigQ) -> BigQ {
        BigQ::new(&self.num * &o.den + &o.num * &self.den, &self.den * &o.den)
    }

    pub(crate) fn neg(&self) -> BigQ {
        BigQ {
            num: -&self.num,
            den: self.den.clone(),
        }
    }

    pub(crate) fn mul(&self, o: &BigQ) -> BigQ {
        BigQ::new(&self.num * &o.num, &self.den * &o.den)
    }

    /// `None` when `o` is zero.
    pub(crate) fn div(&self, o: &BigQ) -> Option<BigQ> {
        (!o.num.is_zero()).then(|| BigQ::new(&self.num * &o.den, &self.den * &o.num))
    }

    pub(crate) fn cmp_value(&self, o: &BigQ) -> Ordering {
        (&self.num * &o.den).cmp(&(&o.num * &self.den))
    }

    pub(crate) fn bits(&self) -> u64 {
        self.num.bits().max(self.den.bits())
    }

    pub(crate) fn to_f64(&self) -> f64 {
        ratio_to_f64(&self.num, &self.den)
    }

    /// The exact value of a finite f64.
    pub(crate) fn from_f64(x: f64) -> BigQ {
        let bits = x.to_bits();
        let biased = ((bits >> 52) & 0x7ff) as i64;
        let fraction = bits & ((1u64 << 52) - 1);
        let (mant, exp) = match biased {
            0 => (fraction, -1074),
            e => (fraction | (1 << 52), e - 1075),
        };
        let signed = match x.is_sign_negative() {
            true => -BigInt::from(mant),
            false => BigInt::from(mant),
        };
        match exp >= 0 {
            true => BigQ::new(signed << exp as usize, BigInt::one()),
            false => BigQ::new(signed, BigInt::one() << (-exp) as usize),
        }
    }
}

/// `num / den` rounded to the nearest f64, ties to even, subnormals and overflow included.
pub(crate) fn ratio_to_f64(num: &BigInt, den: &BigInt) -> f64 {
    let negative = (num.sign() == Sign::Minus) != (den.sign() == Sign::Minus);
    let a = num.magnitude();
    let b = den.magnitude();
    let magnitude = match a.is_zero() {
        true => 0.0,
        false => magnitude_to_f64(a, b),
    };
    match negative {
        true => -magnitude,
        false => magnitude,
    }
}

fn magnitude_to_f64(a: &BigUint, b: &BigUint) -> f64 {
    let shift = 54 - (a.bits() as i64 - b.bits() as i64);
    let (n, d) = match shift >= 0 {
        true => (a << shift as usize, b.clone()),
        false => (a.clone(), b << (-shift) as usize),
    };
    let (q, r) = n.div_rem(&d);
    let qbits = q.bits() as i64;
    let drop = (qbits - 53).max(shift - 1074);
    if drop > qbits + 1 {
        return 0.0;
    }
    let kept = (&q >> drop as usize).to_u64().unwrap_or(0);
    let round = q.bit(drop as u64 - 1);
    let sticky = !r.is_zero() || q.trailing_zeros().is_some_and(|t| (t as i64) < drop - 1);
    let up = round && (sticky || kept & 1 == 1);
    ldexp(kept + up as u64, drop - shift)
}

/// `m * 2^e` for `m <= 2^53`, exact unless it overflows to infinity.
fn ldexp(m: u64, e: i64) -> f64 {
    let top = e + 63 - m.leading_zeros() as i64;
    match (m == 0, top > 1023) {
        (true, _) => 0.0,
        (false, true) => f64::INFINITY,
        (false, false) => {
            let half = e / 2;
            m as f64 * pow2(half) * pow2(e - half)
        }
    }
}

fn pow2(e: i64) -> f64 {
    match e {
        e if e >= -1022 => f64::from_bits(((e + 1023) as u64) << 52),
        e => f64::from_bits(1u64 << (e + 1074)),
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn rounding_to_f64_reaches_subnormals_ties_and_overflow() {
        let one = BigInt::one();
        let min_sub = f64::from_bits(1);
        assert_eq!(ratio_to_f64(&one, &(BigInt::one() << 1074)), min_sub);
        assert_eq!(ratio_to_f64(&one, &(BigInt::one() << 1075)), 0.0);
        let three = BigInt::from(3);
        assert_eq!(ratio_to_f64(&three, &(BigInt::one() << 1076)), min_sub);
        assert_eq!(ratio_to_f64(&one, &BigInt::from(3)), 1.0 / 3.0);
        let tie = (BigInt::one() << 53) + 1;
        assert_eq!(ratio_to_f64(&tie, &one), 9007199254740992.0);
        let big = BigInt::one() << 1024;
        assert_eq!(ratio_to_f64(&big, &one), f64::INFINITY);
        assert_eq!(ratio_to_f64(&-big, &one), f64::NEG_INFINITY);
        let max = BigQ::from_f64(f64::MAX);
        assert_eq!(max.to_f64(), f64::MAX);
    }

    #[test]
    fn every_f64_converts_to_its_exact_value_and_back() {
        for x in [0.1, -2.5e-310, 1e300, f64::MIN_POSITIVE, 123456.789, -7.0] {
            assert_eq!(BigQ::from_f64(x).to_f64().to_bits(), x.to_bits(), "{x}");
        }
        assert!(BigQ::from_f64(-0.0).num.is_zero());
    }
}