alkahest_cas/numeric/
pslq.rs

1//! Detect near-linear dependencies among floating scalars via an augmented lattice + LLL.
2//!
3//! Rows `eᵢ ⊕ ⌊β·xᵢ⌋` spanning ℤⁿ⁺¹ are reduced with [`crate::lattice::lattice_reduce_rows`]; short
4//! vectors correlate with approximate integer relations \(\sumᵢ aᵢ xᵢ ≈ 0\).
5
6use crate::errors::AlkahestError;
7use crate::lattice::{lattice_reduce_rows, LatticeError};
8use rug::ops::PowAssign;
9use rug::{Assign, Float, Integer};
10use std::cmp::Ordering;
11use std::fmt;
12
13/// Errors from [`guess_integer_relation`].
14#[derive(Debug, Clone, PartialEq, Eq)]
15pub enum PslqError {
16    TooFewCoordinates,
17    AllZeroMagnitudes,
18    PrecisionTooThin { bits: u32 },
19    Lattice(LatticeError),
20}
21
22impl fmt::Display for PslqError {
23    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
24        match self {
25            PslqError::TooFewCoordinates => write!(
26                f,
27                "guess_integer_relation needs at least two floating scalars"
28            ),
29            PslqError::AllZeroMagnitudes => write!(
30                f,
31                "scaled magnitudes vanished (check precision or literals)"
32            ),
33            PslqError::PrecisionTooThin { bits } => {
34                write!(f, "precision_bits ({bits}); require ≥64 MPFR bits")
35            }
36            PslqError::Lattice(e) => write!(f, "{e}"),
37        }
38    }
39}
40
41impl std::error::Error for PslqError {
42    fn source(&self) -> Option<&(dyn std::error::Error + 'static)> {
43        match self {
44            PslqError::Lattice(e) => Some(e),
45            _ => None,
46        }
47    }
48}
49
50impl AlkahestError for PslqError {
51    fn code(&self) -> &'static str {
52        match self {
53            PslqError::TooFewCoordinates => "E-PSLQ-001",
54            PslqError::AllZeroMagnitudes => "E-PSLQ-002",
55            PslqError::PrecisionTooThin { .. } => "E-PSLQ-003",
56            PslqError::Lattice(e) => e.code(),
57        }
58    }
59
60    fn remediation(&self) -> Option<&'static str> {
61        match self {
62            PslqError::TooFewCoordinates => Some("pass [x₀,…,x_{n−1}] with n ≥ 2"),
63            PslqError::AllZeroMagnitudes => {
64                Some("use higher-precision inputs (strings or MPFR literals)")
65            }
66            PslqError::PrecisionTooThin { .. } => {
67                Some("raise precision_bits — ≈664 bits ≈ 200 decimal digits")
68            }
69            PslqError::Lattice(e) => e.remediation(),
70        }
71    }
72}
73
74fn lin_residual(bits: u32, coeffs: &[Integer], xs: &[Float]) -> Float {
75    let mut acc = Float::with_val(bits, 0);
76    for (c, xv) in coeffs.iter().zip(xs.iter()) {
77        let mut term = Float::with_val(bits, c);
78        term *= Float::with_val(bits, xv);
79        acc += term;
80    }
81    acc.abs_mut();
82    acc
83}
84
85/// Search for integers `(a₀,…,a_{n−1})` with \(|\sum_i a_i x_i|\) below a precision-derived threshold.
86///
87/// * `max_abs_coeff` — optional filter rejecting candidates with any `|a_i|` above the bound.
88pub fn guess_integer_relation(
89    xs: &[Float],
90    precision_bits: u32,
91    max_abs_coeff: Option<u128>,
92) -> Result<Option<Vec<Integer>>, PslqError> {
93    let n = xs.len();
94    if n < 2 {
95        return Err(PslqError::TooFewCoordinates);
96    }
97    if precision_bits < 64 {
98        return Err(PslqError::PrecisionTooThin {
99            bits: precision_bits,
100        });
101    }
102    let bits = precision_bits.min(16_384);
103
104    let mut normed: Vec<Float> = xs.iter().map(|xv| Float::with_val(bits, xv)).collect();
105    let mut ymax = Float::with_val(bits, 0);
106    for v in &normed {
107        let mut cp = Float::with_val(bits, v);
108        cp.abs_mut();
109        if cp.partial_cmp(&ymax) == Some(Ordering::Greater) {
110            ymax.assign(&cp);
111        }
112    }
113    let zero = Float::with_val(bits, 0);
114    if ymax.partial_cmp(&zero) == Some(Ordering::Equal) {
115        return Err(PslqError::AllZeroMagnitudes);
116    }
117
118    for v in &mut normed {
119        let mut quot = Float::with_val(bits, &*v);
120        quot /= &ymax;
121        v.assign(&quot);
122    }
123
124    let shift_amt = (bits * 3 / 4).min(1536);
125    let mut scale = Integer::from(1u32);
126    scale <<= shift_amt;
127
128    let mut augmented: Vec<Vec<Integer>> = Vec::with_capacity(n);
129    for i in 0..n {
130        let mut row = vec![Integer::from(0); n + 1];
131        row[i] = Integer::from(1);
132        let mut xf = Float::with_val(bits, &normed[i]);
133        xf *= Float::with_val(bits, &scale);
134        let tail = xf.to_integer().ok_or(PslqError::AllZeroMagnitudes)?;
135        row[n].assign(&tail);
136        augmented.push(row);
137    }
138
139    let reduced = lattice_reduce_rows(&augmented).map_err(PslqError::Lattice)?;
140
141    let mut tol = Float::with_val(bits, 2);
142    let exp_lim = ((-((bits as f64) * 0.75).floor()) as i32).min(-1);
143    tol.pow_assign(exp_lim);
144    tol *= Float::with_val(bits, (n.max(1)) as i32);
145
146    let mut best: Option<(Vec<Integer>, Float)> = None;
147    for row in &reduced {
148        let coeffs: Vec<Integer> = row.iter().take(n).cloned().collect();
149        if coeffs.iter().all(Integer::is_zero) {
150            continue;
151        }
152        if let Some(limit) = max_abs_coeff {
153            let lim = Integer::from(limit);
154            let mut ok = true;
155            for z in &coeffs {
156                let mut a = z.clone();
157                a.abs_mut();
158                if a.cmp(&lim) == Ordering::Greater {
159                    ok = false;
160                    break;
161                }
162            }
163            if !ok {
164                continue;
165            }
166        }
167        let resid = lin_residual(bits, &coeffs, &normed);
168        let take = match &best {
169            None => true,
170            Some((_, r0)) => resid.partial_cmp(r0) == Some(Ordering::Less),
171        };
172        if take {
173            best = Some((coeffs, resid));
174        }
175    }
176
177    Ok(best.and_then(|(v, r)| {
178        if r.partial_cmp(&tol) != Some(Ordering::Greater) {
179            Some(v)
180        } else {
181            None
182        }
183    }))
184}
185
186#[cfg(test)]
187mod tests {
188    use super::*;
189
190    #[test]
191    fn relation_on_1_2_3() {
192        let bits = 256u32;
193        let xs = vec![
194            Float::with_val(bits, 1),
195            Float::with_val(bits, 2),
196            Float::with_val(bits, 3),
197        ];
198        let rel = guess_integer_relation(&xs, bits, Some(10_000))
199            .unwrap()
200            .unwrap();
201        let r = lin_residual(bits, &rel, &xs);
202        let mut tol = Float::with_val(bits, 2);
203        tol.pow_assign(-((bits as f64 * 0.75).floor() as i32));
204        tol *= Float::with_val(bits, 3);
205        assert!(
206            r.partial_cmp(&tol) != Some(Ordering::Greater),
207            "residual {r:?} tol {tol:?}"
208        );
209    }
210}
alkahest_cas/numeric/pslq.rs

alkahest_cas/numeric/
pslq.rs
