adele_ring/
ball.rs

1//! The Archimedean (infinite) place, made into a type.
2//!
3//! A [`Ball`] is a rigorous rational enclosure `[lo, hi]` of a real number — never
4//! a point estimate. **Every** ordering / sign / magnitude decision in the tower
5//! goes through `Ball`, because the p-adic (finite) channels are constitutionally
6//! blind to size: they measure divisibility, not magnitude.
7//!
8//! When [`Ball::sign`] or [`Ball::cmp`] returns `None` the enclosure straddles the
9//! boundary and the caller is obligated to *refine* (bisect an algebraic interval,
10//! or pull more digits from a computable real) and retry. This is the only honest
11//! way to decide the sign of an exact real, and it replaces every place the old
12//! code silently reconstructed through `BigInt`.
13//!
14//! The `lo`/`hi` rationals stay *small*: they are low-precision enclosures whose
15//! `BigInt` content is bounded by the requested precision, not by the basis
16//! modulus `M`.
17
18use std::cmp::Ordering;
19
20use num_bigint::{BigInt, Sign};
21use num_rational::BigRational;
22use num_traits::{One, Signed, ToPrimitive, Zero};
23
24/// Exact rational type used for `Ball` endpoints.
25pub type Rational = BigRational;
26
27/// A rigorous rational enclosure of a real number.
28#[derive(Clone, Debug, PartialEq, Eq)]
29pub struct Ball {
30    pub lo: Rational,
31    pub hi: Rational,
32}
33
34impl Ball {
35    /// An enclosure `[lo, hi]` (endpoints are swapped if given out of order).
36    pub fn new(lo: Rational, hi: Rational) -> Self {
37        if lo <= hi {
38            Ball { lo, hi }
39        } else {
40            Ball { lo: hi, hi: lo }
41        }
42    }
43
44    /// A zero-width ball at the exact rational `r`.
45    pub fn point(r: &Rational) -> Self {
46        Ball { lo: r.clone(), hi: r.clone() }
47    }
48
49    /// A zero-width ball at an integer.
50    pub fn from_i64(n: i64) -> Self {
51        Self::point(&BigRational::from_integer(BigInt::from(n)))
52    }
53
54    /// Sign of the enclosed value, or `None` if the ball straddles zero
55    /// (the caller must refine and retry).
56    pub fn sign(&self) -> Option<Sign> {
57        if self.lo.is_positive() {
58            Some(Sign::Plus)
59        } else if self.hi.is_negative() {
60            Some(Sign::Minus)
61        } else if self.lo.is_zero() && self.hi.is_zero() {
62            Some(Sign::NoSign)
63        } else {
64            None
65        }
66    }
67
68    /// Ordering against another ball, or `None` if the two overlap (refine).
69    pub fn cmp(&self, other: &Ball) -> Option<Ordering> {
70        if self.hi < other.lo {
71            Some(Ordering::Less)
72        } else if self.lo > other.hi {
73            Some(Ordering::Greater)
74        } else if self.lo == self.hi && other.lo == other.hi && self.lo == other.lo {
75            Some(Ordering::Equal)
76        } else {
77            None
78        }
79    }
80
81    /// Whether the enclosure contains zero.
82    pub fn contains_zero(&self) -> bool {
83        !self.lo.is_positive() && !self.hi.is_negative()
84    }
85
86    /// Whether the enclosure contains the rational `v`.
87    pub fn contains(&self, v: &Rational) -> bool {
88        self.lo <= *v && *v <= self.hi
89    }
90
91    /// Width `hi - lo`.
92    pub fn width(&self) -> Rational {
93        &self.hi - &self.lo
94    }
95
96    /// Midpoint `(lo + hi) / 2`.
97    pub fn midpoint(&self) -> Rational {
98        (&self.lo + &self.hi) / BigRational::from_integer(BigInt::from(2))
99    }
100
101    /// Nearest `f64` to the midpoint (endpoints are small, so this is finite).
102    pub fn to_f64(&self) -> f64 {
103        rational_to_f64(&self.midpoint())
104    }
105
106    /// Interval addition.
107    pub fn add(&self, b: &Ball) -> Ball {
108        Ball::new(&self.lo + &b.lo, &self.hi + &b.hi)
109    }
110
111    /// Interval negation.
112    pub fn neg(&self) -> Ball {
113        Ball::new(-&self.hi, -&self.lo)
114    }
115
116    /// Interval subtraction.
117    pub fn sub(&self, b: &Ball) -> Ball {
118        Ball::new(&self.lo - &b.hi, &self.hi - &b.lo)
119    }
120
121    /// Interval multiplication (min/max over the four endpoint products).
122    pub fn mul(&self, b: &Ball) -> Ball {
123        let products = [
124            &self.lo * &b.lo,
125            &self.lo * &b.hi,
126            &self.hi * &b.lo,
127            &self.hi * &b.hi,
128        ];
129        let mut min = products[0].clone();
130        let mut max = products[0].clone();
131        for p in &products[1..] {
132            if *p < min {
133                min = p.clone();
134            }
135            if *p > max {
136                max = p.clone();
137            }
138        }
139        Ball::new(min, max)
140    }
141
142    /// Interval reciprocal, or `None` when the enclosure contains zero.
143    pub fn recip(&self) -> Option<Ball> {
144        if self.contains_zero() {
145            return None;
146        }
147        let one = BigRational::one();
148        Some(Ball::new(&one / &self.hi, &one / &self.lo))
149    }
150}
151
152/// Convert a (small) rational to `f64` via bit-shift scaling so neither endpoint
153/// overflows during the division.
154pub fn rational_to_f64(r: &Rational) -> f64 {
155    let p = r.numer();
156    let q = r.denom();
157    if q.is_zero() {
158        return f64::NAN;
159    }
160    if p.is_zero() {
161        return 0.0;
162    }
163    let negative = (p.sign() == Sign::Minus) ^ (q.sign() == Sign::Minus);
164    let mut pm = p.magnitude().clone();
165    let mut qm = q.magnitude().clone();
166    let pb = pm.bits() as i64;
167    let qb = qm.bits() as i64;
168    let shift = 60 + qb - pb;
169    if shift > 0 {
170        pm <<= shift as u64;
171    } else {
172        qm <<= (-shift) as u64;
173    }
174    let quo = &pm / &qm;
175    let mag = quo.to_f64().unwrap_or(f64::INFINITY);
176    let value = mag * 2f64.powi(-(shift as i32));
177    if negative {
178        -value
179    } else {
180        value
181    }
182}
183
184/// Build a rational from an `i64` ratio.
185pub fn rat(p: i64, q: i64) -> Rational {
186    BigRational::new(BigInt::from(p), BigInt::from(q))
187}
188
189#[cfg(test)]
190mod tests {
191    use super::*;
192
193    #[test]
194    fn sign_and_straddle() {
195        assert_eq!(Ball::new(rat(1, 2), rat(3, 2)).sign(), Some(Sign::Plus));
196        assert_eq!(Ball::new(rat(-3, 2), rat(-1, 2)).sign(), Some(Sign::Minus));
197        assert_eq!(Ball::new(rat(-1, 1), rat(1, 1)).sign(), None);
198        assert_eq!(Ball::from_i64(0).sign(), Some(Sign::NoSign));
199    }
200
201    #[test]
202    fn ordering() {
203        let a = Ball::from_i64(1);
204        let b = Ball::from_i64(2);
205        assert_eq!(a.cmp(&b), Some(Ordering::Less));
206        assert_eq!(b.cmp(&a), Some(Ordering::Greater));
207        // overlapping intervals → None
208        assert_eq!(Ball::new(rat(0, 1), rat(2, 1)).cmp(&Ball::new(rat(1, 1), rat(3, 1))), None);
209    }
210
211    #[test]
212    fn interval_arith() {
213        let a = Ball::new(rat(1, 1), rat(2, 1));
214        let b = Ball::new(rat(3, 1), rat(4, 1));
215        assert_eq!(a.add(&b), Ball::new(rat(4, 1), rat(6, 1)));
216        assert_eq!(a.mul(&b), Ball::new(rat(3, 1), rat(8, 1)));
217        assert!(a.recip().is_some());
218        assert!(Ball::new(rat(-1, 1), rat(1, 1)).recip().is_none());
219    }
220
221    #[test]
222    fn to_float() {
223        assert!((Ball::new(rat(1, 3), rat(1, 3)).to_f64() - (1.0 / 3.0)).abs() < 1e-15);
224    }
225}
adele_ring/ball.rs

adele_ring/
ball.rs
