retex_arith/

lib.rs

//! TeX-exact fixed-point arithmetic.
//!
//! TeX does all dimension arithmetic with integers. A dimension is a *scaled*
//! value: an `i32` counting units of 2⁻¹⁶ pt (`sp`). Results must match
//! TeX bit for bit, so every routine here follows the specification in
//! *TeX: The Program* (§99–§108) and the e-TeX expression routines
//! (`\numexpr`, `\dimexpr`) exactly, including the overflow rules.
//!
//! Overflow is reported as [`Overflow`]. TeX sets its global `arith_error`
//! flag in that case and the caller decides what to do; callers that need
//! TeX's exact fallback values get them from the documented error variants.

#![forbid(unsafe_code)]

use std::fmt;

pub mod random;

/// One point in scaled units (2¹⁶ sp).
pub const UNITY: i32 = 0x1_0000;
/// Two points in scaled units, used by decimal rounding (§102).
pub const TWO: i32 = 0x2_0000;
/// The largest legal dimension, `\maxdimen` = 16383.99999 pt (§421).
pub const MAX_DIMEN: i32 = 0o7777777777;
/// The largest legal integer, 2³¹ − 1 (TeX's `infinity`).
pub const INFINITY: i32 = 0o17777777777;
/// Badness of an infinitely bad box (§108).
pub const INF_BAD: i32 = 10000;

/// An arithmetic overflow, equivalent to TeX setting `arith_error`.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub struct Overflow;

impl fmt::Display for Overflow {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        f.write_str("arithmetic overflow")
    }
}

impl std::error::Error for Overflow {}

/// A scaled dimension: an integer number of sp (2⁻¹⁶ pt).
#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash, Default)]
pub struct Scaled(pub i32);

impl Scaled {
    /// Zero.
    pub const ZERO: Scaled = Scaled(0);
    /// `\maxdimen`.
    pub const MAX: Scaled = Scaled(MAX_DIMEN);

    /// Formats the value the way TeX's `print_scaled` does (§103), without
    /// the `pt` unit. For example, 65536 sp becomes `"1.0"`.
    pub fn to_tex_string(self) -> String {
        let mut out = String::new();
        print_scaled(&mut out, self.0);
        out
    }
}

impl fmt::Display for Scaled {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        let mut s = String::new();
        print_scaled(&mut s, self.0);
        f.write_str(&s)
    }
}

/// `half(x)` (§100): x/2 rounded, with odd values rounded towards +∞.
pub fn half(x: i32) -> i32 {
    if x & 1 != 0 {
        // Pascal `div` truncates towards zero, like Rust's `/`.
        (x + 1) / 2
    } else {
        x / 2
    }
}

/// `round_decimals` (§102): converts the decimal fraction `.d₀d₁…dₖ₋₁`
/// to scaled units, rounded. `digits` holds the digit values (0–9). TeX
/// only looks at the first 17 digits; callers must truncate to that.
pub fn round_decimals(digits: &[u8]) -> i32 {
    let mut a: i32 = 0;
    for &d in digits.iter().rev() {
        a = (a + i32::from(d) * TWO) / 10;
    }
    (a + 1) / 2
}

/// Appends TeX's decimal form of a scaled value to `out` (§103). It prints
/// the shortest decimal that rounds back to the same scaled value, always
/// with at least one digit after the point.
pub fn print_scaled(out: &mut String, s: i32) {
    if s == i32::MIN {
        // TeX never produces this value, but Lua can store it (e.g.
        // `pdf.setlinkmargin(2^31)`); LuaTeX's negation wraps around and it
        // prints a doubled minus sign.
        out.push_str("--32768.0");
        return;
    }
    let mut s = i64::from(s);
    if s < 0 {
        out.push('-');
        s = -s;
    }
    let unity = i64::from(UNITY);
    out.push_str(&(s / unity).to_string());
    out.push('.');
    s = 10 * (s % unity) + 5;
    let mut delta: i64 = 10;
    loop {
        if delta > unity {
            // Round the last digit.
            s += 0o100000 - 50000;
        }
        let digit = (s / unity) as u8; // always 0..=9
        out.push(char::from(b'0' + digit));
        s = 10 * (s % unity);
        delta *= 10;
        if s <= delta {
            break;
        }
    }
}

/// `mult_and_add` (§105): computes `n*x + y`, failing if the absolute
/// value would exceed `max_answer`.
pub fn mult_and_add(mut n: i32, mut x: i32, y: i32, max_answer: i32) -> Result<i32, Overflow> {
    if n < 0 {
        x = x.wrapping_neg();
        n = -n;
    }
    if n == 0 {
        return Ok(y);
    }
    // TeX's test, done in i64 to avoid intermediate overflow in Rust.
    let (n64, x64, y64, m) = (
        i64::from(n),
        i64::from(x),
        i64::from(y),
        i64::from(max_answer),
    );
    if x64 <= (m - y64) / n64 && -x64 <= (m + y64) / n64 {
        Ok((n64 * x64 + y64) as i32)
    } else {
        Err(Overflow)
    }
}

/// `nx_plus_y(n, x, y)` (§105) for dimensions: the result must be at most
/// 2³⁰ − 1 in absolute value.
pub fn nx_plus_y(n: i32, x: i32, y: i32) -> Result<i32, Overflow> {
    mult_and_add(n, x, y, 0o7777777777)
}

/// `mult_integers(n, x)` (§105): the result must fit in 2³¹ − 1.
pub fn mult_integers(n: i32, x: i32) -> Result<i32, Overflow> {
    mult_and_add(n, x, 0, 0o17777777777)
}

/// Quotient and remainder from [`x_over_n`] and [`xn_over_d`].
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub struct DivResult {
    /// The quotient.
    pub quotient: i32,
    /// The remainder, which has the sign TeX gives it.
    pub remainder: i32,
}

/// `x_over_n` (§106): divides scaled `x` by integer `n`, truncating towards
/// zero. Division by zero is an overflow; in that case TeX returns quotient 0
/// and remainder `x`, available via [`x_over_n_unchecked`].
pub fn x_over_n(x: i32, n: i32) -> Result<DivResult, Overflow> {
    if n == 0 {
        return Err(Overflow);
    }
    Ok(x_over_n_unchecked(x, n))
}

/// [`x_over_n`] without the error, returning TeX's fallback
/// `{quotient: 0, remainder: x}` when `n == 0`.
pub fn x_over_n_unchecked(mut x: i32, mut n: i32) -> DivResult {
    if n == 0 {
        return DivResult {
            quotient: 0,
            remainder: x,
        };
    }
    let mut negative = false;
    if n < 0 {
        x = x.wrapping_neg();
        n = n.wrapping_neg();
        negative = true;
    }
    let (q, mut r) = if x >= 0 {
        (x / n, x % n)
    } else {
        // -x cannot overflow for TeX-legal values; wrap to be total.
        let nx = x.wrapping_neg();
        (-(nx / n), -(nx % n))
    };
    if negative {
        r = -r;
    }
    DivResult {
        quotient: q,
        remainder: r,
    }
}

/// `xn_over_d` (§107): computes `x*n/d` for `0 ≤ n, d ≤ 2¹⁶`, truncating
/// towards zero. TeX simulates 48-bit arithmetic with 15-bit pieces; this is
/// a direct transcription of that algorithm so that overflow is detected in
/// exactly the same cases.
pub fn xn_over_d(x: i32, n: i32, d: i32) -> Result<DivResult, Overflow> {
    const B: i64 = 0o100000; // 2^15
    let positive = x >= 0;
    let x = i64::from(x).abs();
    let (n, d) = (i64::from(n), i64::from(d));
    let t = (x % B) * n;
    let mut u = (x / B) * n + (t / B);
    let v = (u % d) * B + (t % B);
    if u / d >= B {
        return Err(Overflow);
    }
    u = B * (u / d) + (v / d);
    let (q, r) = (u as i32, (v % d) as i32);
    Ok(if positive {
        DivResult {
            quotient: q,
            remainder: r,
        }
    } else {
        DivResult {
            quotient: -q,
            remainder: -r,
        }
    })
}

/// `badness(t, s)` (§108): approximately 100(t/s)³, capped at
/// [`INF_BAD`]. `t` is the excess or shortfall, `s` the available stretch or
/// shrink.
pub fn badness(t: i32, s: i32) -> i32 {
    if t == 0 {
        return 0;
    }
    if s <= 0 {
        return INF_BAD;
    }
    let r: i32 = if t <= 7230584 {
        (t * 297) / s
    } else if s >= 1663497 {
        t / (s / 297)
    } else {
        t
    };
    if r > 1290 {
        INF_BAD
    } else {
        (r * r * r + 0o400000) / 0o1000000
    }
}

/// e-TeX `add_or_sub`: computes `x ± y`, failing if the absolute value of
/// the result would exceed `max_answer`.
pub fn add_or_sub(x: i32, mut y: i32, max_answer: i32, negative: bool) -> Result<i32, Overflow> {
    if negative {
        y = y.wrapping_neg();
    }
    if x >= 0 {
        if y <= max_answer - x {
            Ok(x + y)
        } else {
            Err(Overflow)
        }
    } else if y >= -max_answer - x {
        Ok(x + y)
    } else {
        Err(Overflow)
    }
}

/// e-TeX `quotient(n, d)`: the rounded quotient ⌊n/d + ½⌋, with the sign
/// applied symmetrically.
pub fn quotient(mut n: i32, mut d: i32) -> Result<i32, Overflow> {
    if d == 0 {
        return Err(Overflow);
    }
    let mut negative = if d > 0 {
        false
    } else {
        d = -d;
        true
    };
    if n < 0 {
        n = n.wrapping_neg();
        negative = !negative;
    }
    let mut a = n / d;
    n -= a * d;
    d = n - d;
    if d + n >= 0 {
        a += 1;
    }
    Ok(if negative { -a } else { a })
}

/// e-TeX `fract(x, n, d, max_answer)`: ⌊xn/d + ½⌋ computed without
/// overflow, failing if the result exceeds `max_answer`. This is a
/// transcription of the e-TeX algorithm.
pub fn fract(mut x: i32, mut n: i32, mut d: i32, max_answer: i32) -> Result<i32, Overflow> {
    if d == 0 {
        return Err(Overflow);
    }
    let mut a: i32 = 0;
    let mut negative = if d > 0 {
        false
    } else {
        d = -d;
        true
    };
    if x < 0 {
        x = x.wrapping_neg();
        negative = !negative;
    } else if x == 0 {
        return Ok(0);
    }
    if n < 0 {
        n = n.wrapping_neg();
        negative = !negative;
    }
    let finish = |a: i32, negative: bool| if negative { -a } else { a };
    let mut t = n / d;
    if t > max_answer / x {
        return Err(Overflow);
    }
    a = a.wrapping_add(t * x);
    n -= t * d;
    if n == 0 {
        return Ok(finish(a, negative));
    }
    t = x / d;
    if t > (max_answer - a) / n {
        return Err(Overflow);
    }
    a += t * n;
    x -= t * d;
    if x == 0 {
        return Ok(finish(a, negative));
    }
    if x < n {
        std::mem::swap(&mut x, &mut n);
    }
    // Compute f = floor(xn/d + 1/2) with 0 < n <= x < d.
    let mut f: i32 = 0;
    let mut r: i32 = (d / 2) - d;
    let h = -r;
    loop {
        if n & 1 != 0 {
            r += x;
            if r >= 0 {
                r -= d;
                f += 1;
            }
        }
        n /= 2;
        if n == 0 {
            break;
        }
        if x < h {
            x += x;
        } else {
            t = x - d;
            x += t;
            f += n;
            if x < n {
                if x == 0 {
                    break;
                }
                std::mem::swap(&mut x, &mut n);
            }
        }
    }
    if f > max_answer - a {
        return Err(Overflow);
    }
    a += f;
    Ok(finish(a, negative))
}

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

    #[test]
    fn print_scaled_examples() {
        let cases = [
            (0, "0.0"),
            (UNITY, "1.0"),
            (-UNITY, "-1.0"),
            (UNITY / 2, "0.5"),
            (1, "0.00002"),
            (MAX_DIMEN, "16383.99998"),
            (655360, "10.0"),
            (26214, "0.4"),
            (32768 + 1, "0.50002"),
            // As LuaTeX prints it.
            (i32::MIN, "--32768.0"),
        ];
        for (v, s) in cases {
            assert_eq!(Scaled(v).to_tex_string(), s, "value {v}");
        }
    }

    #[test]
    fn round_decimals_examples() {
        assert_eq!(round_decimals(&[5]), 32768);
        assert_eq!(round_decimals(&[2, 5]), 16384);
        assert_eq!(round_decimals(&[]), 0);
        assert_eq!(round_decimals(&[9, 9, 9, 9, 9, 9]), UNITY);
    }

    #[test]
    fn half_rounds_up_on_odd() {
        assert_eq!(half(3), 2);
        assert_eq!(half(-3), -1);
        assert_eq!(half(4), 2);
        assert_eq!(half(-4), -2);
    }

    #[test]
    fn badness_examples() {
        assert_eq!(badness(0, 0), 0);
        assert_eq!(badness(1, 0), INF_BAD);
        assert_eq!(badness(UNITY, UNITY), 100);
        assert_eq!(badness(UNITY, 2 * UNITY), 12);
        assert_eq!(badness(3 * UNITY, UNITY), 2698);
        assert_eq!(badness(5 * UNITY, UNITY), INF_BAD);
    }

    #[test]
    fn overflow_cases() {
        assert_eq!(nx_plus_y(2, MAX_DIMEN, 0), Err(Overflow));
        assert_eq!(x_over_n(5, 0), Err(Overflow));
        assert_eq!(quotient(5, 0), Err(Overflow));
        assert_eq!(fract(5, 1, 0, INFINITY), Err(Overflow));
        assert_eq!(add_or_sub(INFINITY, 1, INFINITY, false), Err(Overflow));
    }

    proptest! {
        /// print_scaled must produce a decimal that round_decimals maps back
        /// to the same value (TeX's documented guarantee, §103).
        #[test]
        fn print_scaled_round_trips(v in -MAX_DIMEN..=MAX_DIMEN) {
            let s = Scaled(v).to_tex_string();
            let (neg, body) = match s.strip_prefix('-') {
                Some(b) => (true, b),
                None => (false, s.as_str()),
            };
            let (int, frac) = body.split_once('.').unwrap();
            let digits: Vec<u8> = frac.bytes().map(|b| b - b'0').collect();
            let back = int.parse::<i32>().unwrap() * UNITY + round_decimals(&digits);
            prop_assert_eq!(if neg { -back } else { back }, v);
        }

        #[test]
        fn x_over_n_truncates(x in -MAX_DIMEN..=MAX_DIMEN, n in -100_000i32..100_000) {
            prop_assume!(n != 0);
            let r = x_over_n(x, n).unwrap();
            prop_assert_eq!(r.quotient, x / n);
            // TeX's remainder carries the sign of x, like Rust's `%`.
            prop_assert_eq!(r.remainder, x % n);
        }

        #[test]
        fn xn_over_d_matches_wide(x in -MAX_DIMEN..=MAX_DIMEN, n in 0i32..=0x10000, d in 1i32..=0x10000) {
            let wide = i64::from(x) * i64::from(n);
            let q = wide / i64::from(d);
            match xn_over_d(x, n, d) {
                Ok(r) => {
                    prop_assert_eq!(i64::from(r.quotient), q);
                    prop_assert_eq!(i64::from(r.remainder), wide % i64::from(d));
                }
                Err(Overflow) => prop_assert!(q.abs() >= 1 << 30),
            }
        }

        #[test]
        fn fract_matches_wide(x in -INFINITY..=INFINITY, n in -INFINITY..=INFINITY, d in -INFINITY..=INFINITY) {
            prop_assume!(d != 0);
            let (ax, an, ad) = (i128::from(x).abs(), i128::from(n).abs(), i128::from(d).abs());
            let q = (2 * ax * an + ad) / (2 * ad);
            let neg = (x < 0) ^ (n < 0) ^ (d < 0);
            match fract(x, n, d, INFINITY) {
                Ok(v) => prop_assert_eq!(i128::from(v), if neg && x != 0 { -q } else { q }),
                Err(Overflow) => prop_assert!(q > i128::from(INFINITY)),
            }
        }

        #[test]
        fn quotient_rounds_half_away(n in -INFINITY..=INFINITY, d in -INFINITY..=INFINITY) {
            prop_assume!(d != 0);
            let (an, ad) = (i64::from(n).abs(), i64::from(d).abs());
            let q = (2 * an + ad) / (2 * ad);
            let neg = (n < 0) ^ (d < 0);
            prop_assert_eq!(i64::from(quotient(n, d).unwrap()), if neg { -q } else { q });
        }
    }
}