retex-arith 0.9.9

//! METAFONT's random number generator, which LuaTeX uses for
//! `\uniformdeviate`, `\normaldeviate` and the `tex` library's randomizers.
//!
//! Everything follows *METAFONT: The Program* (§107–§152): the `fraction`
//! arithmetic it needs ([`make_fraction`], [`take_fraction`], [`m_log`],
//! [`ab_vs_cd`]) and the generator itself ([`Randoms`]: `init_randoms`,
//! `new_randoms`, `unif_rand`, `norm_rand`). A `fraction` is an integer
//! with 28 bits after the binary point, so [`FRACTION_ONE`] stands for 1.
//!
//! The arithmetic routines compute the values the book specifies with
//! 64-bit intermediates instead of its 32-bit loops; the results are the
//! same. On overflow the book sets `arith_error` and returns ±`el_gordo`;
//! these routines return the same clamped value (the generator never
//! overflows for arguments of magnitude below 2³¹).

use std::cmp::Ordering;

/// 2²⁸, the `fraction` that represents 1.0.
pub const FRACTION_ONE: i32 = 1 << 28;
/// 2²⁷, the `fraction` that represents 0.5.
pub const FRACTION_HALF: i32 = 1 << 27;
/// 2³⁰, the `fraction` that represents 4.0.
pub const FRACTION_FOUR: i32 = 1 << 30;
/// 2³¹ − 1, the largest value METAFONT likes (`el_gordo`).
pub const EL_GORDO: i32 = i32::MAX;

/// Clamps a result to ±`el_gordo`, as METAFONT does when it sets
/// `arith_error`.
fn clamp(x: i64) -> i32 {
    x.clamp(-i64::from(EL_GORDO), i64::from(EL_GORDO)) as i32
}

/// `make_fraction(p, q)` (§107): ⌊2²⁸p/q + ½⌋ for positive operands, with
/// the sign of p/q otherwise.
///
/// # Panics
///
/// If `q` is zero (METAFONT's `confusion("/")`).
pub fn make_fraction(p: i32, q: i32) -> i32 {
    assert!(q != 0, "make_fraction: division by zero");
    let negative = (p < 0) != (q < 0);
    let (p, q) = (i64::from(p).abs(), i64::from(q).abs());
    if p / q >= 8 {
        return if negative { -EL_GORDO } else { EL_GORDO };
    }
    let f = ((p << 29) + q) / (2 * q);
    clamp(if negative { -f } else { f })
}

/// `take_fraction(q, f)` (§109): ⌊qf/2²⁸ + ½⌋ for positive operands, with
/// the sign of qf otherwise.
pub fn take_fraction(q: i32, f: i32) -> i32 {
    clamp(take_fraction_wide(i64::from(q), i64::from(f)))
}

/// [`take_fraction`] for any operands that fit in 32 bits in magnitude,
/// without clamping.
fn take_fraction_wide(q: i64, f: i64) -> i64 {
    let negative = (q < 0) != (f < 0);
    let p = (q.abs() * f.abs() + i64::from(FRACTION_HALF)) >> 28;
    if negative { -p } else { p }
}

/// The values of 2²⁷·ln(1/(1 − 2⁻ᵏ)), rounded, used by [`m_log`] (§131).
const SPEC_LOG: [i32; 29] = {
    let mut t = [0; 29];
    let head = [
        93032640, 38612034, 17922280, 8662214, 4261238, 2113709, 1052693, 525315, 262400, 131136,
        65552, 32772, 16385,
    ];
    let mut k = 1;
    while k <= 13 {
        t[k] = head[k - 1];
        k += 1;
    }
    while k <= 27 {
        t[k] = 1 << (27 - k);
        k += 1;
    }
    t[28] = 1;
    t
};

/// `m_log(x)` (§132): 2⁸ times the natural logarithm of the `scaled`
/// value `x`, i.e. an integer approximation of 2²⁴·ln(x/2¹⁶), for `x > 0`.
/// METAFONT reports an error for `x ≤ 0` and returns 0, as this does.
pub fn m_log(mut x: i32) -> i32 {
    if x <= 0 {
        return 0;
    }
    let mut y: i32 = 1_302_456_956 + 4 - 100;
    let mut z: i32 = 27595 + 6_553_600;
    while x < FRACTION_FOUR {
        x += x;
        y -= 93_032_639;
        z -= 48782;
    }
    y += z / crate::UNITY;
    let mut k = 2;
    while x > FRACTION_FOUR + 4 {
        // Increase k until x can be multiplied by a factor of 2^-k.
        let mut z = (x - 1) / (1 << k) + 1;
        while x < FRACTION_FOUR + z {
            z = (z + 1) / 2;
            k += 1;
        }
        y += SPEC_LOG[k];
        x -= z;
    }
    y / 8
}

/// `ab_vs_cd(a, b, c, d)` (§117): how `a·b` compares with `c·d`.
pub fn ab_vs_cd(a: i32, b: i32, c: i32, d: i32) -> Ordering {
    (i64::from(a) * i64::from(b)).cmp(&(i64::from(c) * i64::from(d)))
}

/// The state of the generator: METAFONT's `randoms` and `j_random`
/// (§148).
#[derive(Clone, Debug, PartialEq, Eq)]
pub struct Randoms {
    /// The last 55 random fractions generated.
    randoms: [i32; 55],
    /// The number of unused entries of `randoms`.
    j: usize,
}

impl Randoms {
    /// A generator initialized with `seed` (`init_randoms`, §150).
    pub fn new(seed: i32) -> Self {
        let mut r = Randoms {
            randoms: [0; 55],
            j: 0,
        };
        r.init(seed);
        r
    }

    /// Restarts the generator from `seed` (`init_randoms`, §150).
    pub fn init(&mut self, seed: i32) {
        // |seed| does not fit for i32::MIN; LuaTeX then produces fractions
        // out of range, which nothing relies on.
        let mut j = i64::from(seed).abs();
        while j >= i64::from(FRACTION_ONE) {
            j /= 2;
        }
        let mut j = j as i32;
        let mut k: i32 = 1;
        for i in 0..55 {
            let jj = k;
            k = j - k;
            j = jj;
            if k < 0 {
                k += FRACTION_ONE;
            }
            self.randoms[(i * 21) % 55] = j;
        }
        // "Warm up" the array.
        self.new_randoms();
        self.new_randoms();
        self.new_randoms();
    }

    /// `new_randoms` (§149): the next 55 values of
    /// xₙ = (xₙ₋₅₅ − xₙ₋₂₄) mod 2²⁸.
    fn new_randoms(&mut self) {
        let r = &mut self.randoms;
        for k in 0..55 {
            let mut x = r[k] - r[if k < 24 { k + 31 } else { k - 24 }];
            if x < 0 {
                x += FRACTION_ONE;
            }
            r[k] = x;
        }
        self.j = 54;
    }

    /// The next random fraction, 0 ≤ f < [`FRACTION_ONE`] (`next_random`,
    /// §148: the table is consumed backwards).
    pub fn next_fraction(&mut self) -> i32 {
        if self.j == 0 {
            self.new_randoms();
        } else {
            self.j -= 1;
        }
        self.randoms[self.j]
    }

    /// `unif_rand(x)` (§151): a uniform deviate u with 0 ≤ u < x,
    /// 0 ≥ u > x, or u = 0 when x = 0. The bounds are inclusive at 0 and
    /// exclusive at x, and 0 has about half the probability of the other
    /// values because of rounding.
    pub fn unif_rand(&mut self, x: i32) -> i32 {
        let f = self.next_fraction();
        let ax = i64::from(x).abs();
        let y = take_fraction_wide(ax, i64::from(f));
        if y == ax {
            0
        } else if x > 0 {
            y as i32
        } else {
            (-y) as i32
        }
    }

    /// `norm_rand` (§152): a normal deviate with mean 0 and standard
    /// deviation 1, as a `scaled` value (ratio method, Algorithm 3.4.1R of
    /// *The Art of Computer Programming*).
    pub fn norm_rand(&mut self) -> i32 {
        loop {
            let (x, u) = loop {
                // 2^16·sqrt(8/e) ≈ 112428.82793.
                let x = take_fraction(112_429, self.next_fraction() - FRACTION_HALF);
                let u = self.next_fraction();
                if x.abs() < u {
                    break (x, u);
                }
            };
            let x = make_fraction(x, u);
            // 2^24·12·ln 2 ≈ 139548959.6165.
            let l = 139_548_960 - m_log(u);
            if ab_vs_cd(1024, l, x, x) != Ordering::Less {
                return x;
            }
        }
    }
}

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

    // The sequences below were printed by LuaTeX 1.22 (`tex.init_rand`
    // followed by `tex.uniform_rand` or `tex.normal_rand`, and
    // `\setrandomseed` with `\uniformdeviate`).

    fn uniform(seed: i32, x: i32, n: usize) -> Vec<i32> {
        let mut r = Randoms::new(seed);
        (0..n).map(|_| r.unif_rand(x)).collect()
    }

    fn normal(seed: i32, n: usize) -> Vec<i32> {
        let mut r = Randoms::new(seed);
        (0..n).map(|_| r.norm_rand()).collect()
    }

    #[test]
    fn uniform_seed_42() {
        assert_eq!(
            uniform(42, 1000, 30),
            [
                788, 356, 360, 48, 154, 905, 286, 377, 754, 921, 893, 58, 813, 784, 283, 776, 46,
                483, 570, 839, 427, 743, 537, 577, 15, 614, 88, 711, 425, 755
            ]
        );
        assert_eq!(
            uniform(42, -1000, 10),
            [-788, -356, -360, -48, -154, -905, -286, -377, -754, -921]
        );
        assert_eq!(
            uniform(42, 6, 30),
            [
                5, 2, 2, 0, 1, 5, 2, 2, 5, 0, 5, 0, 5, 5, 2, 5, 0, 3, 3, 5, 3, 4, 3, 3, 0, 4, 1, 4,
                3, 5
            ]
        );
        assert_eq!(uniform(42, 1, 10), [0; 10]);
        assert_eq!(uniform(42, 0, 3), [0; 3]);
        assert_eq!(
            uniform(42, EL_GORDO, 5),
            [1691999567, 764279616, 772600120, 103650880, 330031416]
        );
        assert_eq!(
            uniform(42, -EL_GORDO, 5),
            [-1691999567, -764279616, -772600120, -103650880, -330031416]
        );
    }

    #[test]
    fn uniform_other_seeds() {
        // The seed's sign does not matter, and seeds of 2^28 or more are
        // halved until they are smaller.
        assert_eq!(uniform(-42, 1000, 3), [788, 356, 360]);
        assert_eq!(uniform(43, 1000, 3), [253, 728, 550]);
        assert_eq!(uniform(0, 1000, 5), [241, 724, 352, 544, 553]);
        assert_eq!(uniform(1, 1000, 5), [706, 96, 542, 580, 544]);
        assert_eq!(uniform(2, 1000, 5), [172, 468, 733, 615, 534]);
        assert_eq!(uniform(FRACTION_ONE + 5, 1000, 3), [672, 968, 733]);
        assert_eq!(uniform(FRACTION_ONE, 1000, 1), [741]);
        assert_eq!(uniform(EL_GORDO, 1000, 5), [775, 352, 161, 508, 563]);
        assert_eq!(uniform(-EL_GORDO, 1000, 1), [775]);
        assert_eq!(uniform(720_667_377, 1_000_000, 1), [565328]);
        assert_eq!(
            uniform(123_456_789, 1_000_000_007, 60),
            [
                339280105, 694342525, 33406381, 663534847, 31468977, 920872412, 890775731,
                763103345, 352619304, 970706429, 556395658, 661018212, 990435626, 129798145,
                863619408, 619866, 169370619, 823904878, 327600645, 795750054, 24217479, 720953581,
                705383831, 74538809, 206004822, 709324257, 411170676, 944563850, 472524974,
                992473296, 501359921, 564391203, 182562610, 884105465, 115442634, 164019470,
                3255386, 306177771, 471044626, 532215923, 638790694, 964107289, 74024737,
                964409153, 666308233, 83044526, 4268754, 977702509, 471425347, 884187835,
                276770061, 37269928, 402504205, 659495034, 479703087, 780795758, 541605458,
                103255720, 732023453, 240494633
            ]
        );
    }

    #[test]
    fn uniform_raw_fractions() {
        // unif_rand(2^28) is the fraction itself; 200 values cross several
        // refills of the table.
        let got = uniform(1234, FRACTION_ONE, 200);
        let mut r = Randoms::new(1234);
        let raw: Vec<i32> = (0..200).map(|_| r.next_fraction()).collect();
        assert_eq!(got, raw);
        assert_eq!(
            got[..12],
            [
                148933938, 868208, 172117591, 200264536, 217942135, 6639995, 144635461, 178301072,
                146807680, 191625536, 108818557, 117198514
            ]
        );
        assert_eq!(
            got[190..],
            [
                254766369, 19665803, 176208551, 174348460, 140378967, 3606447, 115579941, 15397593,
                34731790, 30432757
            ]
        );
    }

    #[test]
    fn normal_sequences() {
        assert_eq!(
            normal(42, 24),
            [
                90948, -63802, 31008, 44963, -31499, -105635, 9407, -11063, 7160, -88755, -65191,
                -11236, 42221, -119113, -145102, -121987, -103122, 112162, 117308, 18788, -106915,
                4377, -88893, -24580
            ]
        );
        assert_eq!(
            normal(7, 40),
            [
                -477, 26353, -2376, 34964, 138853, -737, -63343, 19393, 16883, -40185, -48533,
                40548, -16297, 4018, -63948, 32782, -10939, 19134, -112276, 29490, -45353, 41991,
                60635, 97826, -31185, -6217, -46115, -18400, -47706, -5467, 6379, -12410, -38074,
                -37842, 44554, -97469, 27547, -36117, 42533, -25108
            ]
        );
    }

    #[test]
    fn normal_after_uniforms() {
        // \setrandomseed 42, three \uniformdeviate 1000, two \normaldeviate.
        let mut r = Randoms::new(42);
        let u: Vec<i32> = (0..3).map(|_| r.unif_rand(1000)).collect();
        assert_eq!(u, [788, 356, 360]);
        assert_eq!([r.norm_rand(), r.norm_rand()], [-18408, 112861]);
        assert_eq!(r.unif_rand(-1000), -776);
        // One uniform, then a normal, then a uniform.
        let mut r = Randoms::new(42);
        assert_eq!(
            [r.unif_rand(1000), r.norm_rand(), r.unif_rand(1000)],
            [788, -45034, 48]
        );
    }

    #[test]
    fn fraction_arithmetic() {
        assert_eq!(make_fraction(1, 2), FRACTION_HALF);
        assert_eq!(make_fraction(-1, 3), -89478485);
        assert_eq!(make_fraction(1, -3), -89478485);
        assert_eq!(make_fraction(8, 1), EL_GORDO);
        assert_eq!(make_fraction(-9, 1), -EL_GORDO);
        assert_eq!(take_fraction(1000, FRACTION_HALF), 500);
        assert_eq!(take_fraction(3, FRACTION_HALF), 2);
        assert_eq!(take_fraction(-3, FRACTION_HALF), -2);
        assert_eq!(take_fraction(EL_GORDO, 4 * FRACTION_ONE), EL_GORDO);
        assert_eq!(ab_vs_cd(2, 3, 1, 6), Ordering::Equal);
        assert_eq!(ab_vs_cd(-2, 3, 1, 6), Ordering::Less);
        // ln 1 = 0, ln 2·2^24 ≈ 11629080.
        assert_eq!(m_log(crate::UNITY), 0);
        assert_eq!(m_log(2 * crate::UNITY), 11629080);
        assert_eq!(m_log(0), 0);
    }
}