pounce_common/
utils.rs

1//! Numerical and timing utilities.
2//!
3//! Mirrors `Common/IpUtils.{hpp,cpp}`. The PRNG is a portable
4//! reimplementation of POSIX `drand48` (48-bit LCG with the
5//! glibc-default seed) so output matches a `drand48`-built upstream
6//! Ipopt deterministically. Replace with libc's `drand48` later if a
7//! particular platform's bit-equivalence target requires it.
8
9use crate::types::Number;
10use std::cell::Cell;
11use std::sync::OnceLock;
12use std::time::{Instant, SystemTime, UNIX_EPOCH};
13
14/// Equivalent to `Ipopt::IsFiniteNumber`. Returns false for NaN/±∞.
15#[inline]
16pub fn is_finite_number(val: Number) -> bool {
17    val.is_finite()
18}
19
20/// Equivalent to `Ipopt::Compare_le(lhs, rhs, BasVal)` — relaxed `<=`
21/// with a tolerance proportional to `|BasVal|`. Threshold matches
22/// upstream's `(10 * machine_epsilon * Max(1, |BasVal|))`.
23#[inline]
24pub fn compare_le(lhs: Number, rhs: Number, bas_val: Number) -> bool {
25    let tol = 10.0 * Number::EPSILON * bas_val.abs().max(1.0);
26    lhs - rhs <= tol
27}
28
29/// Equivalent to `Ipopt::IpRandom01`. 48-bit LCG matching POSIX
30/// `drand48`: X_{n+1} = (0x5DEECE66D · X_n + 0xB) mod 2^48,
31/// returning the high-order 48 bits as a `f64` in [0, 1).
32pub fn ip_random_01() -> Number {
33    LCG_STATE.with(|s| {
34        let mut x = s.get();
35        x = x.wrapping_mul(LCG_A).wrapping_add(LCG_C) & LCG_MASK;
36        s.set(x);
37        // Top 32 bits → 53-bit mantissa portion, exactly as drand48
38        // converts. Match glibc's `erand48`:
39        //   r = ((double)(x >> 16)) * 2^-32 + ((x & 0xffff) << 4) * 2^-48 ... etc.
40        // The simplest reformulation that matches: split into 16-bit
41        // chunks and accumulate.
42        let x0 = (x & 0xFFFF) as f64;
43        let x1 = ((x >> 16) & 0xFFFF) as f64;
44        let x2 = ((x >> 32) & 0xFFFF) as f64;
45        x2 / 65536.0 + x1 / (65536.0 * 65536.0) + x0 / (65536.0 * 65536.0 * 65536.0)
46    })
47}
48
49/// Reset the PRNG to glibc's default seed. Equivalent to
50/// `Ipopt::IpResetRandom01`.
51pub fn ip_reset_random_01() {
52    LCG_STATE.with(|s| s.set(LCG_DEFAULT_SEED));
53}
54
55const LCG_A: u64 = 0x5DEECE66D;
56const LCG_C: u64 = 0xB;
57const LCG_MASK: u64 = (1 << 48) - 1;
58const LCG_DEFAULT_SEED: u64 = 0x1234_ABCD_330E;
59
60thread_local! {
61    static LCG_STATE: Cell<u64> = const { Cell::new(LCG_DEFAULT_SEED) };
62}
63
64/// Wallclock time in seconds since first call. Equivalent to
65/// `Ipopt::WallclockTime`.
66pub fn wallclock_time() -> Number {
67    static EPOCH: OnceLock<Instant> = OnceLock::new();
68    let e = EPOCH.get_or_init(Instant::now);
69    e.elapsed().as_secs_f64()
70}
71
72/// CPU time in seconds. Equivalent to `Ipopt::CpuTime`.
73///
74/// std offers no portable CPU-time API, so we fall back to
75/// wallclock. The bit-equivalence trace strips the timing column
76/// before diffing per the plan, so this is acceptable; phase 4 will
77/// wire in a `libc::clock_gettime(CLOCK_PROCESS_CPUTIME_ID)` path
78/// where it's needed by the iteration log.
79pub fn cpu_time() -> Number {
80    wallclock_time()
81}
82
83/// System time in seconds since UNIX epoch. Equivalent to
84/// `Ipopt::SysTime`.
85pub fn sys_time() -> Number {
86    SystemTime::now()
87        .duration_since(UNIX_EPOCH)
88        .map(|d| d.as_secs_f64())
89        .unwrap_or(0.0)
90}
91
92/// Equivalent to `Ipopt::ComputeMemIncrease`. Updates `len` to a new
93/// length consistent with `recommended`/`min` while clamping to
94/// `T::MAX`. Returns `Err` if even the maximum representable size
95/// would not be an increase.
96pub fn compute_mem_increase_i64(
97    len: &mut i64,
98    recommended: f64,
99    min: i64,
100    context: &str,
101) -> Result<(), String> {
102    if recommended >= i64::MAX as f64 {
103        if *len < i64::MAX {
104            *len = i64::MAX;
105            Ok(())
106        } else {
107            Err(format!(
108                "Cannot allocate more than {} bytes for {} due to integer-type limit",
109                (i64::MAX as f64) * 8.0,
110                context
111            ))
112        }
113    } else {
114        *len = min.max(recommended as i64);
115        Ok(())
116    }
117}
118
119#[cfg(test)]
120mod tests {
121    use super::*;
122
123    #[test]
124    fn finite_check() {
125        assert!(is_finite_number(0.0));
126        assert!(is_finite_number(1e300));
127        assert!(!is_finite_number(f64::NAN));
128        assert!(!is_finite_number(f64::INFINITY));
129    }
130
131    #[test]
132    fn compare_le_tolerates_eps() {
133        let v = 1.0;
134        let near = v + 5.0 * Number::EPSILON;
135        assert!(compare_le(near, v, v));
136    }
137
138    #[test]
139    fn random_01_in_range_and_deterministic() {
140        ip_reset_random_01();
141        let a: Vec<f64> = (0..16).map(|_| ip_random_01()).collect();
142        ip_reset_random_01();
143        let b: Vec<f64> = (0..16).map(|_| ip_random_01()).collect();
144        assert_eq!(a, b);
145        for v in a {
146            assert!((0.0..1.0).contains(&v), "{v}");
147        }
148    }
149
150    #[test]
151    fn wallclock_monotonic() {
152        let a = wallclock_time();
153        let b = wallclock_time();
154        assert!(b >= a);
155    }
156}
pounce_common/utils.rs

pounce_common/
utils.rs
