rustymath_core/
domain.rs

1//! What a kernel computes over: a call's dtype and overflow policy (`Variant`), and the values that an argument or a result can take (`Range`). Formulas (`facts`), recipes and kernels read both.
2use crate::buffer::{DType, FloatWidth};
3
4/// The values other than NaN that the items of an argument or a result can take: each lies from `lo` to `hi`. ¯0 and +0 count as one value, because the sign of a zero isn't promised and can differ by route. Any range can also hold NaN. A rule's range is sound where it holds every result that the object's kernels give, within their accuracy.
5#[derive(Clone, Copy, Debug, PartialEq)]
6pub struct Range { pub lo: f64, pub hi: f64 }
7
8impl Range {
9    /// Every value.
10    pub const ANY: Self = Self { lo: f64::NEG_INFINITY, hi: f64::INFINITY };
11    /// Every value from 0 up.
12    pub const NONNEGATIVE: Self = Self { lo: 0.0, hi: f64::INFINITY };
13    /// No value, as an array of no items holds. It's `join`'s identity: its `lo` is ∞ and its `hi` ¯∞.
14    pub const EMPTY: Self = Self { lo: f64::INFINITY, hi: f64::NEG_INFINITY };
15    /// The range of an argument that is the one number `c`, or `ANY` where `c` is NaN or `None`, as for an argument that isn't one number.
16    #[inline]
17    pub fn of(c: Option<f64>) -> Self { c.filter(|c| !c.is_nan()).map_or(Self::ANY, |c| Self { lo: c, hi: c }) }
18    /// The values that items of `dtype` hold: 0 and 1 for Booleans, from each integer width's least integer to its greatest, and any value for other dtypes.
19    #[inline]
20    pub fn of_dtype(dtype: DType) -> Self {
21        match dtype {
22            DType::Boolean => Self { lo: 0.0, hi: 1.0 }, DType::Integer(width) => with_width!(width, T => Self { lo: T::MIN as f64, hi: T::MAX as f64 }),
23            _ => Self::ANY
24        }
25    }
26    /// The one number that the range holds, if it holds only one.
27    #[inline]
28    pub fn constant(self) -> Option<f64> { (self.lo == self.hi).then_some(self.lo) }
29    /// Whether the range holds no negative number.
30    #[inline]
31    pub fn nonnegative(self) -> bool { self.lo >= 0.0 }
32    /// Whether the range holds no value (`EMPTY`).
33    #[inline]
34    pub fn is_empty(self) -> bool { self.lo > self.hi }
35    /// Whether every value lies from `lo` to `hi`.
36    #[inline]
37    pub fn within(self, lo: f64, hi: f64) -> bool { lo <= self.lo && self.hi <= hi }
38    /// Whether no value equals `c`.
39    #[inline]
40    pub fn excludes(self, c: f64) -> bool { c < self.lo || self.hi < c }
41    /// The values that both ranges hold.
42    #[inline]
43    pub fn meet(self, other: Self) -> Self { Self { lo: self.lo.max(other.lo), hi: self.hi.min(other.hi) } }
44    /// The least range that holds both ranges.
45    #[inline]
46    pub fn join(self, other: Self) -> Self { Self { lo: self.lo.min(other.lo), hi: self.hi.max(other.hi) } }
47    /// The range of |x| for the values x of this range.
48    pub fn magnitude(self) -> Self {
49        if self.lo >= 0.0 { self } else if self.hi <= 0.0 { Self { lo: -self.hi, hi: -self.lo } } else { Self { lo: 0.0, hi: self.hi.max(-self.lo) } }
50    }
51    /// The range of x × x for the values x of this range, as `⊼x` squares them.
52    pub fn square(self) -> Self { let m = self.magnitude(); Self { lo: m.lo * m.lo, hi: m.hi * m.hi } }
53    /// The range of a fold's results, from its items' range, its seed's, if it has one, and `step`, which gives the range of a step's result from the ranges of its two arguments. Every lane holds an item. One bound holds for any fold: the least range that holds the items and the seed, where a step's results stay in it. Where they don't, as a sum's don't, each bound widens to an infinity, except that a bound on the far side of 0 stops at 0: [2, 3] widens to [0, ∞]. Where a step's results leave that range too, that bound is `ANY`.
54    ///
55    /// Where `lane` gives the most items that a lane holds and `step` is an associative function's, as `+ × ⌈ ⌊ ∧ ∨` are, a second bound holds, and the range is the meet of both. It starts at the items' range, and joins in `step` of the range on itself while the lane length that it covers doubles, up to `lane`: [2, 3] over lanes of up to 4 items gives [2, 12]. It holds because a lane of up to 2m items splits into two lanes of up to m items whose results, by associativity, one step combines, and `step` bounds that step for any arguments within its two ranges. A scan's running results are the results of shorter lanes. A seed then combines with that bound in one step, as associativity lets it combine with a lane's result, and its own range joins in for a lane of no items.
56    pub fn folded(items: Self, seed: Option<Self>, lane: Option<usize>, step: impl Fn([Self; 2]) -> Self) -> Self {
57        let start = seed.map_or(items, |s| s.join(items));
58        let (lo, hi) = (if start.lo >= 0.0 { 0.0 } else { f64::NEG_INFINITY }, if start.hi <= 0.0 { 0.0 } else { f64::INFINITY });
59        let kept = [start, Self { lo, hi }].into_iter().find(|&r| step([r, r]).inside(r)).unwrap_or(Self::ANY);
60        let Some(lane) = lane else { return kept };
61        let (mut covered, mut m) = (items, 1);
62        while m < lane {
63            covered = covered.join(step([covered, covered]));
64            m = m.saturating_mul(2);
65        }
66        seed.map_or(covered, |s| s.join(step([s, covered]))).meet(kept)
67    }
68    /// Whether `other` holds every value of this range.
69    fn inside(self, other: Self) -> bool { other.lo <= self.lo && self.hi <= other.hi }
70    /// This range with each bound moved `ULPS` units in the last place outward, for bounds that a kernel computes within its accuracy, as the transcendental functions do. A bound never moves past 0, and 0 or an infinite bound stays.
71    pub(crate) fn loose(self) -> Self { Self { lo: outward(self.lo, false, FloatWidth::F64), hi: outward(self.hi, true, FloatWidth::F64) } }
72    /// This range as a kernel's results of `dtype` hold it. The bounds are f64, and a kernel at a narrower float rounds each result to that width, within its accuracy there: each bound moves `ULPS` units in that width's last place outward, as `loose` moves it at f64. A kernel at an integer dtype gives whole results, flooring a quotient as integer `÷` does: each bound rounds outward to a whole number. Other dtypes keep the range.
73    pub fn at_width(self, dtype: DType) -> Self {
74        match dtype {
75            DType::Float(width @ (FloatWidth::F32 | FloatWidth::F16)) => Self { lo: outward(self.lo, false, width), hi: outward(self.hi, true, width) },
76            DType::Integer(_) => Self { lo: self.lo.floor(), hi: self.hi.ceil() }, _ => self
77        }
78    }
79}
80
81/// How far a computed bound can lie from a kernel's result on the same argument, in units in the last place: a kernel lies within 4 of the correctly rounded result, and so does the function that computes the bound.
82const ULPS: u64 = 8;
83
84/// `x` moved `ULPS` units in the last place of floats of `width` up, or down, but not past 0. A zero or an infinity stays. A bound between two floats of `width` starts from the nearer, which the move covers.
85fn outward(x: f64, up: bool, width: FloatWidth) -> f64 {
86    if x == 0.0 || x.is_infinite() { return x; }
87    let (bits, infinity): (u64, u64) = match width {
88        FloatWidth::F64 => (x.abs().to_bits(), f64::INFINITY.to_bits()), FloatWidth::F32 => ((x.abs() as f32).to_bits().into(), f32::INFINITY.to_bits().into()),
89        FloatWidth::F16 => (half::f16::from_f64(x.abs()).to_bits().into(), half::f16::INFINITY.to_bits().into())
90    };
91    let moved = if up == x.is_sign_positive() { (bits + ULPS).min(infinity) } else { bits.saturating_sub(ULPS) };
92    let moved = match width {
93        FloatWidth::F64 => f64::from_bits(moved), FloatWidth::F32 => f32::from_bits(moved as u32).into(),
94        FloatWidth::F16 => half::f16::from_bits(moved as u16).to_f64()
95    };
96    moved.copysign(x)
97}
98
99/// Policies captured when the call is built; never read from force-time preferences.
100#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash)]
101pub struct Variant { pub dtype: DType, pub checked: bool, pub wrapping: bool }
102
103impl Variant {
104    #[inline]
105    pub const fn checked(dtype: DType) -> Self { Self { dtype, checked: true, wrapping: false } }
106    /// Whether an integer result that doesn't fit its type fails, as checked arithmetic that doesn't wrap gives it, rather than wrapping or filling. It covers a computed result and a result narrowed to a smaller output type.
107    #[inline]
108    pub const fn fails_on_overflow(self) -> bool { self.checked && !self.wrapping }
109}
rustymath_core/domain.rs

rustymath_core/
domain.rs
