axiolid_exact/
fixed.rs

1//! Fixed-point intervals at a chosen precision, with certified `sin` and
2//! `cos`: the tier for values no finite arithmetic holds exactly.
3//!
4//! A [`FixedInterval`] is `[lo, hi] * 2^-bits` with big-integer bounds.
5//! Every operation rounds its bounds outward, so the true value stays
6//! inside. Unlike [`Dyadic`] the result is not exact, but
7//! its width shrinks as `bits` grows: a sign the interval cannot decide
8//! at one precision is asked again at a higher one, and any nonzero value
9//! is decided at some precision. Deciding that a value is exactly zero is
10//! the caller's business (by an identity, as for harmonics of dyadic
11//! angles).
12//!
13//! `sin` and `cos` of a dyadic angle are summed from their Taylor series
14//! after halving the angle below `1/16`, with the Lagrange remainder added
15//! to the bounds, and doubled back up in interval arithmetic.
16
17use num_bigint::BigInt;
18use num_integer::Integer;
19
20use axiolid_guarantees::Sign;
21
22use crate::dyadic::Dyadic;
23use crate::interval::Interval;
24
25/// A real number between `lo * 2^-bits` and `hi * 2^-bits`.
26#[derive(Debug, Clone, PartialEq, Eq)]
27pub struct FixedInterval {
28    lo: BigInt,
29    hi: BigInt,
30    bits: u32,
31}
32
33/// `x / 2^k`, rounded up.
34fn ceil_shr(x: &BigInt, k: u32) -> BigInt {
35    -((-x) >> k)
36}
37
38impl FixedInterval {
39    /// Exactly `value` (an integer), at `bits` fractional bits.
40    #[must_use]
41    pub fn integer(value: i64, bits: u32) -> Self {
42        let v = BigInt::from(value) << bits;
43        Self {
44            lo: v.clone(),
45            hi: v,
46            bits,
47        }
48    }
49
50    /// The narrowest interval at `bits` holding `value`.
51    #[must_use]
52    pub fn from_dyadic(value: &Dyadic, bits: u32) -> Self {
53        let shift = value.exponent() + i64::from(bits);
54        let m = value.mantissa();
55        if shift >= 0 {
56            let v = m << (shift as u64);
57            return Self {
58                lo: v.clone(),
59                hi: v,
60                bits,
61            };
62        }
63        let k = u32::try_from(-shift).unwrap_or(u32::MAX);
64        Self {
65            lo: m >> k,
66            hi: ceil_shr(m, k),
67            bits,
68        }
69    }
70
71    /// The narrowest interval at `bits` holding a finite `value`.
72    #[must_use]
73    pub fn from_f64(value: f64, bits: u32) -> Option<Self> {
74        Dyadic::try_from_f64(value).map(|d| Self::from_dyadic(&d, bits))
75    }
76
77    /// The fractional bits of the bounds.
78    #[must_use]
79    pub fn bits(&self) -> u32 {
80        self.bits
81    }
82
83    fn same(&self, other: &Self) {
84        debug_assert_eq!(self.bits, other.bits, "fixed intervals of one precision");
85    }
86
87    /// The sum.
88    #[must_use]
89    pub fn add(&self, other: &Self) -> Self {
90        self.same(other);
91        Self {
92            lo: &self.lo + &other.lo,
93            hi: &self.hi + &other.hi,
94            bits: self.bits,
95        }
96    }
97
98    /// The difference.
99    #[must_use]
100    pub fn sub(&self, other: &Self) -> Self {
101        self.add(&other.neg())
102    }
103
104    /// The negation.
105    #[must_use]
106    pub fn neg(&self) -> Self {
107        Self {
108            lo: -&self.hi,
109            hi: -&self.lo,
110            bits: self.bits,
111        }
112    }
113
114    /// The product, rounded outward.
115    #[must_use]
116    pub fn mul(&self, other: &Self) -> Self {
117        self.same(other);
118        let products = [
119            &self.lo * &other.lo,
120            &self.lo * &other.hi,
121            &self.hi * &other.lo,
122            &self.hi * &other.hi,
123        ];
124        let lo = products.iter().min().expect("four products");
125        let hi = products.iter().max().expect("four products");
126        Self {
127            lo: lo >> self.bits,
128            hi: ceil_shr(hi, self.bits),
129            bits: self.bits,
130        }
131    }
132
133    /// The product with an integer, exactly.
134    #[must_use]
135    pub fn mul_int(&self, factor: i128) -> Self {
136        let factor = BigInt::from(factor);
137        let (a, b) = (&self.lo * &factor, &self.hi * &factor);
138        if a <= b {
139            Self {
140                lo: a,
141                hi: b,
142                bits: self.bits,
143            }
144        } else {
145            Self {
146                lo: b,
147                hi: a,
148                bits: self.bits,
149            }
150        }
151    }
152
153    /// The quotient by a positive integer, rounded outward.
154    #[must_use]
155    pub fn div_int(&self, divisor: u64) -> Self {
156        let d = BigInt::from(divisor.max(1));
157        Self {
158            lo: self.lo.div_floor(&d),
159            hi: -((-&self.hi).div_floor(&d)),
160            bits: self.bits,
161        }
162    }
163
164    /// The same interval at `bits` fractional bits, rounded outward.
165    #[must_use]
166    pub fn with_bits(&self, bits: u32) -> Self {
167        if bits >= self.bits {
168            let k = bits - self.bits;
169            return Self {
170                lo: &self.lo << k,
171                hi: &self.hi << k,
172                bits,
173            };
174        }
175        let k = self.bits - bits;
176        Self {
177            lo: &self.lo >> k,
178            hi: ceil_shr(&self.hi, k),
179            bits,
180        }
181    }
182
183    /// The sign every value in the interval shares: `Zero` only for the
184    /// single point zero, `None` where the interval straddles zero.
185    #[must_use]
186    pub fn sign(&self) -> Option<Sign> {
187        let zero = BigInt::from(0);
188        if self.lo > zero {
189            Some(Sign::Positive)
190        } else if self.hi < zero {
191            Some(Sign::Negative)
192        } else if self.lo == zero && self.hi == zero {
193            Some(Sign::Zero)
194        } else {
195            None
196        }
197    }
198
199    /// A sound `f64` interval holding this one.
200    #[must_use]
201    pub fn enclosure(&self) -> Interval {
202        let exponent = -i64::from(self.bits);
203        let lo = Dyadic::from_parts(self.lo.clone(), exponent).enclosure();
204        let hi = Dyadic::from_parts(self.hi.clone(), exponent).enclosure();
205        Interval::from_bounds(lo.lo(), hi.hi())
206    }
207
208    /// Bounds on the magnitude of every value in the interval, as sound
209    /// `f64`s: `(0, _)` where the interval holds zero.
210    #[must_use]
211    pub fn magnitude(&self) -> (f64, f64) {
212        let e = self.enclosure();
213        let upper = e.lo().abs().max(e.hi().abs());
214        let lower = if e.lo() > 0.0 {
215            e.lo()
216        } else if e.hi() < 0.0 {
217            -e.hi()
218        } else {
219            0.0
220        };
221        (lower, upper)
222    }
223
224    /// The interval cut to `[-1, 1]`: sound for a value known to lie
225    /// there, such as a sine.
226    fn clamp_unit(mut self) -> Self {
227        let one = BigInt::from(1) << self.bits;
228        if self.hi > one {
229            self.hi = one.clone();
230        }
231        if self.lo < -&one {
232            self.lo = -one;
233        }
234        self
235    }
236
237    /// Certified `sin(x)` and `cos(x)` at `bits` fractional bits; `None`
238    /// for an angle beyond `2^40` in magnitude.
239    #[must_use]
240    pub fn sin_cos(x: &Dyadic, bits: u32) -> Option<(Self, Self)> {
241        let magnitude = x.to_f64().abs();
242        if magnitude > 2f64.powi(40) {
243            return None;
244        }
245        // Halvings that bring the angle below 1/16.
246        let mut halvings = 0u32;
247        while magnitude * 16.0 > 2f64.powi(halvings as i32) {
248            halvings += 1;
249        }
250        // Each doubling may double the error: guard bits for them and for
251        // the series' own rounding.
252        let w = bits + halvings + 32;
253        let y = if x.mantissa() == &BigInt::from(0) {
254            Self::integer(0, w)
255        } else {
256            let halved =
257                Dyadic::from_parts(x.mantissa().clone(), x.exponent() - i64::from(halvings));
258            Self::from_dyadic(&halved, w)
259        };
260        let y2 = y.mul(&y);
261        let one_ulp = BigInt::from(1);
262        let series = |mut term: Self, first: u64| {
263            // Terms `y^(n) / n!` with alternating signs; `first` is the
264            // power of the first term (0 for cos, 1 for sin).
265            let mut sum = term.clone();
266            let mut n = first;
267            loop {
268                term = term.mul(&y2).div_int((n + 1) * (n + 2)).neg();
269                n += 2;
270                let (_, size) = term.magnitude_ulps();
271                if size <= one_ulp {
272                    // The Lagrange remainder is at most this next term's
273                    // magnitude: widen by it instead of adding it.
274                    let pad = size + &one_ulp;
275                    sum.lo -= &pad;
276                    sum.hi += &pad;
277                    return sum;
278                }
279                sum = sum.add(&term);
280            }
281        };
282        let mut s = series(y.clone(), 1);
283        let mut c = series(Self::integer(1, w), 0);
284        for _ in 0..halvings {
285            let s2 = s.mul(&c).mul_int(2);
286            let c2 = Self::integer(1, w).sub(&s.mul(&s).mul_int(2));
287            s = s2.clamp_unit();
288            c = c2.clamp_unit();
289        }
290        Some((s.with_bits(bits), c.with_bits(bits)))
291    }
292
293    /// `(lower, upper)` bounds on the magnitude in units of `2^-bits`.
294    fn magnitude_ulps(&self) -> (BigInt, BigInt) {
295        let zero = BigInt::from(0);
296        let (a, b) = (self.lo.clone(), self.hi.clone());
297        let upper = if -&a > b { -&a } else { b.clone() };
298        let lower = if a > zero {
299            a
300        } else if b < zero {
301            -b
302        } else {
303            zero
304        };
305        (lower, upper)
306    }
307}
308
309#[cfg(test)]
310mod tests {
311    use super::*;
312
313    fn d(x: f64) -> Dyadic {
314        Dyadic::try_from_f64(x).unwrap()
315    }
316
317    #[test]
318    fn floors_and_ceilings_round_outward() {
319        let x = FixedInterval::from_f64(-0.1, 8).unwrap();
320        // -0.1 * 256 = -25.6
321        assert_eq!(
322            (x.lo.clone(), x.hi.clone()),
323            (BigInt::from(-26), BigInt::from(-25))
324        );
325        let p = x.mul(&x);
326        // 0.01 * 256 = 2.56; products 625..676 over 256.
327        assert!(p.lo <= BigInt::from(2) && p.hi >= BigInt::from(3));
328        let q = FixedInterval::integer(7, 4).div_int(3);
329        // 7/3 * 16 = 37.33
330        assert_eq!((q.lo, q.hi), (BigInt::from(37), BigInt::from(38)));
331    }
332
333    #[test]
334    fn sine_and_cosine_of_one_hold_their_digits() {
335        // floor(sin(1) * 2^200) and floor(cos(1) * 2^200), from exact
336        // rational Taylor sums.
337        let sin1: BigInt = "1352191738627887730370568382426388711294672994500390224088559"
338            .parse()
339            .unwrap();
340        let cos1: BigInt = "868232330700371202471720823065555340568207867417093696894962"
341            .parse()
342            .unwrap();
343        let (s, c) = FixedInterval::sin_cos(&d(1.0), 200).unwrap();
344        assert!(s.lo <= sin1 && sin1 < s.hi, "{s:?}");
345        assert!(c.lo <= cos1 && cos1 < c.hi, "{c:?}");
346        assert!(&s.hi - &s.lo < BigInt::from(8), "{s:?}");
347        assert!(&c.hi - &c.lo < BigInt::from(8), "{c:?}");
348    }
349
350    #[test]
351    fn large_and_tiny_angles_agree_with_f64() {
352        for x in [0.0, 1e-300, -3e-9, 0.75, -2.5, 7.0, 100.25, -12345.678] {
353            let (s, c) = FixedInterval::sin_cos(&d(x), 120).unwrap();
354            let (es, ec) = (s.enclosure(), c.enclosure());
355            assert!(
356                es.lo() <= x.sin() + 1e-15 && x.sin() - 1e-15 <= es.hi(),
357                "{x}"
358            );
359            assert!(
360                ec.lo() <= x.cos() + 1e-15 && x.cos() - 1e-15 <= ec.hi(),
361                "{x}"
362            );
363            // sin^2 + cos^2 = 1, to a few units of the last place.
364            let one = s.mul(&s).add(&c.mul(&c));
365            let unit = BigInt::from(1) << 120u32;
366            assert!(one.lo <= unit && unit <= one.hi, "{x}: {one:?}");
367            assert!(&one.hi - &one.lo < BigInt::from(64), "{x}: {one:?}");
368        }
369        assert!(FixedInterval::sin_cos(&d(1e13), 64).is_none());
370    }
371
372    #[test]
373    fn signs_are_certain_or_withheld() {
374        assert_eq!(FixedInterval::integer(0, 10).sign(), Some(Sign::Zero));
375        let (s, _) = FixedInterval::sin_cos(&d(1e-40), 64).unwrap();
376        assert_eq!(s.sign(), None);
377        let (s, _) = FixedInterval::sin_cos(&d(1e-40), 200).unwrap();
378        assert_eq!(s.sign(), Some(Sign::Positive));
379    }
380}
axiolid_exact/fixed.rs

axiolid_exact/
fixed.rs
