implicit_predicates/generated/
points.rs

1// Generated by implicit-predicates-gen from predicates/points.pred. Do not edit.
2// To regenerate, run from the directory that the paths are relative to:
3// implicit-predicates-gen crates/implicit-predicates/predicates --stages crates/implicit-predicates/src/generated --checks crates/implicit-predicates/tests/generated
4
5//! The point types [`Explicit`], [`Lnc`], and [`AnyPoint`], which holds a point
6//! of any of them.
7
8use core::fmt;
9
10use implicit_predicates_core::{
11    ExactCoords, ExactError, Mark, Point, Scratch, ScratchExhausted, Span, Weight, eft, expansion,
12};
13
14/// Why the constructor of a point type rejected its inputs.
15#[derive(Clone, Copy, Debug, PartialEq, Eq)]
16pub enum PointError {
17    /// An input is NaN or infinite.
18    NonFinite,
19    /// The parameters do not satisfy the conditions of the type.
20    ParamOutOfRange,
21}
22
23impl fmt::Display for PointError {
24    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
25        f.write_str(match self {
26            PointError::NonFinite => "an input of the point is not finite",
27            PointError::ParamOutOfRange => "a parameter of the point is out of range",
28        })
29    }
30}
31
32impl core::error::Error for PointError {}
33
34/// A point of type `Explicit`: the point with the coordinates `p`, given
35/// exactly. Its approximation is its coordinates and its error is zero.
36#[derive(Clone, Copy, Debug, PartialEq)]
37#[repr(transparent)]
38pub struct Explicit([f64; 3]);
39
40impl Explicit {
41    /// The point with the coordinates `p`.
42    ///
43    /// # Errors
44    ///
45    /// [`PointError::NonFinite`] if a coordinate is NaN or infinite.
46    #[inline]
47    pub fn new(p: [f64; 3]) -> Result<Self, PointError> {
48        if !(p.iter().all(|x| x.is_finite())) {
49            return Err(PointError::NonFinite);
50        }
51        Ok(Self::new_unchecked(p))
52    }
53
54    /// The point with the coordinates `p`, without checking them. A coordinate that
55    /// is not finite makes [`Point::exact`] fail with [`ExactError::OutOfRange`].
56    #[inline]
57    pub fn new_unchecked(p: [f64; 3]) -> Self {
58        Self(p)
59    }
60}
61
62impl Point for Explicit {
63    #[inline]
64    fn approx(&self) -> &[f64; 3] {
65        &self.0
66    }
67
68    #[inline]
69    fn err(&self) -> &[f64; 3] {
70        &[0.0; 3]
71    }
72
73    #[inline]
74    fn is_explicit(&self) -> bool {
75        true
76    }
77
78    #[inline]
79    fn exact(&self, scratch: &mut Scratch<'_>) -> Result<ExactCoords, ExactError> {
80        self.0.exact(scratch)
81    }
82}
83
84/// A point of type `Lnc`, computed from `p`, `q` and `t`:
85///
86/// ```text
87/// x = q.x + t * (p.x - q.x)
88/// y = q.y + t * (p.y - q.y)
89/// z = q.z + t * (p.z - q.z)
90/// where 0 <= t, t <= 1
91/// ```
92///
93/// The constructors compute the coordinates in `f64` with a bound on the error
94/// of each one, [`Point::approx`] and [`Point::err`]; [`Point::exact`] computes
95/// them exactly.
96#[derive(Clone, Copy, Debug, PartialEq)]
97pub struct Lnc {
98    p: [f64; 3],
99    q: [f64; 3],
100    t: f64,
101    approx: [f64; 3],
102    err: [f64; 3],
103}
104
105impl Lnc {
106    /// The slack factor of the error bound: the smallest `f64` `S` with `S * (1 -
107    /// 2^-53)^6 >= 1`, where 6 is the number of roundings of the computation of the
108    /// error of a coordinate on its longest chain of operations.
109    const ERR_SLACK: f64 = 1.0000000000000009;
110
111    /// The point of type `Lnc` with the parameters `p`, `q` and `t`.
112    ///
113    /// # Errors
114    ///
115    /// [`PointError::NonFinite`] if an input is NaN or infinite, and
116    /// [`PointError::ParamOutOfRange`] if the parameters do not satisfy the
117    /// conditions of the type.
118    pub fn new(p: [f64; 3], q: [f64; 3], t: f64) -> Result<Self, PointError> {
119        if !(p.iter().all(|x| x.is_finite()) && q.iter().all(|x| x.is_finite()) && t.is_finite()) {
120            return Err(PointError::NonFinite);
121        }
122        if t < 0.0 {
123            return Err(PointError::ParamOutOfRange);
124        }
125        if t > 1.0 {
126            return Err(PointError::ParamOutOfRange);
127        }
128        Ok(Self::new_unchecked(p, q, t))
129    }
130
131    /// The point of type `Lnc` with the parameters `p`, `q` and `t`, without
132    /// checking them: any `f64` is accepted. A coordinate whose approximation or
133    /// error bound is not finite gets an error of `+∞`, and [`Point::exact`]
134    /// fails with [`ExactError::OutOfRange`] when an input is not finite.
135    pub fn new_unchecked(p: [f64; 3], q: [f64; 3], t: f64) -> Self {
136        let (approx, err) = Self::evaluate([p, q], [t]);
137        Self {
138            p,
139            q,
140            t,
141            approx,
142            err,
143        }
144    }
145
146    /// The parameter `p`.
147    #[inline]
148    pub fn p(&self) -> &[f64; 3] {
149        &self.p
150    }
151
152    /// The parameter `q`.
153    #[inline]
154    pub fn q(&self) -> &[f64; 3] {
155        &self.q
156    }
157
158    /// The parameter `t`.
159    #[inline]
160    pub fn t(&self) -> f64 {
161        self.t
162    }
163
164    /// The coordinates in `f64` and the bounds on their errors, from the `F64x3`
165    /// parameters `c` and the `f64` parameters `s`, in the order of the
166    /// declaration.
167    fn evaluate(c: [[f64; 3]; 2], s: [f64; 1]) -> ([f64; 3], [f64; 3]) {
168        let v0 = c[1][0];
169        let v1 = c[1][1];
170        let v2 = c[1][2];
171        let v3 = s[0];
172        let v4 = c[0][0];
173        let v5 = c[0][1];
174        let v6 = c[0][2];
175        let v7 = v4 - v0;
176        let e7 = sum_rounding(v4, v0, v7);
177        let v8 = v5 - v1;
178        let e8 = sum_rounding(v5, v1, v8);
179        let v9 = v6 - v2;
180        let e9 = sum_rounding(v6, v2, v9);
181        let v10 = v3 * v7;
182        let e10 = covered_product(v3.abs(), e7) + product_rounding(v3, v7, v10);
183        let v11 = v3 * v8;
184        let e11 = covered_product(v3.abs(), e8) + product_rounding(v3, v8, v11);
185        let v12 = v3 * v9;
186        let e12 = covered_product(v3.abs(), e9) + product_rounding(v3, v9, v12);
187        let v13 = v0 + v10;
188        let e13 = e10 + sum_rounding(v0, v10, v13);
189        let v14 = v1 + v11;
190        let e14 = e11 + sum_rounding(v1, v11, v14);
191        let v15 = v2 + v12;
192        let e15 = e12 + sum_rounding(v2, v12, v15);
193        let approx = [v13, v14, v15];
194        let err = [
195            finish(approx[0], e13, Self::ERR_SLACK),
196            finish(approx[1], e14, Self::ERR_SLACK),
197            finish(approx[2], e15, Self::ERR_SLACK),
198        ];
199        (approx, err)
200    }
201
202    /// Writes the exact coordinates in `scratch`, from the `F64x3` parameters `c`
203    /// and the `f64` parameters `s`, in the order of the declaration. Fails with
204    /// [`ExactError::OutOfRange`] if an input is not finite, or an operand of a sum
205    /// of expansions reaches `2^1020`, or an operand of an exact product leaves the
206    /// product domain; on an error the arena is as it was.
207    fn write_exact(
208        scratch: &mut Scratch<'_>,
209        c: [[f64; 3]; 2],
210        s: [f64; 1],
211    ) -> Result<ExactCoords, ExactError> {
212        if !(c.iter().flatten().all(|x| x.is_finite()) && s.iter().all(|x| x.is_finite())) {
213            return Err(ExactError::OutOfRange);
214        }
215        let mark = scratch.mark();
216        let mut in_range = true;
217        let x = scratch.write(5, |out| {
218            let v0 = c[1][0];
219            let v3 = s[0];
220            let v4 = c[0][0];
221            if !(summable(&[v4]) && summable(&[v0])) {
222                in_range = false;
223                return 0;
224            }
225            let mut x7 = [0.0; 2];
226            let n7 = expansion::from_two_diff(v4, v0, &mut x7);
227            if !(multipliable(&x7[..n7]) && multipliable(&[v3])) {
228                in_range = false;
229                return 0;
230            }
231            let mut x10 = [0.0; 4];
232            let n10 = expansion::scale(&x7[..n7], v3, &mut x10);
233            if !(summable(&x10[..n10]) && summable(&[v0])) {
234                in_range = false;
235                return 0;
236            }
237            expansion::grow(&x10[..n10], v0, out)
238        });
239        let x = written(scratch, mark, x, in_range)?;
240        let y = scratch.write(5, |out| {
241            let v1 = c[1][1];
242            let v3 = s[0];
243            let v5 = c[0][1];
244            if !(summable(&[v5]) && summable(&[v1])) {
245                in_range = false;
246                return 0;
247            }
248            let mut x8 = [0.0; 2];
249            let n8 = expansion::from_two_diff(v5, v1, &mut x8);
250            if !(multipliable(&x8[..n8]) && multipliable(&[v3])) {
251                in_range = false;
252                return 0;
253            }
254            let mut x11 = [0.0; 4];
255            let n11 = expansion::scale(&x8[..n8], v3, &mut x11);
256            if !(summable(&x11[..n11]) && summable(&[v1])) {
257                in_range = false;
258                return 0;
259            }
260            expansion::grow(&x11[..n11], v1, out)
261        });
262        let y = written(scratch, mark, y, in_range)?;
263        let z = scratch.write(5, |out| {
264            let v2 = c[1][2];
265            let v3 = s[0];
266            let v6 = c[0][2];
267            if !(summable(&[v6]) && summable(&[v2])) {
268                in_range = false;
269                return 0;
270            }
271            let mut x9 = [0.0; 2];
272            let n9 = expansion::from_two_diff(v6, v2, &mut x9);
273            if !(multipliable(&x9[..n9]) && multipliable(&[v3])) {
274                in_range = false;
275                return 0;
276            }
277            let mut x12 = [0.0; 4];
278            let n12 = expansion::scale(&x9[..n9], v3, &mut x12);
279            if !(summable(&x12[..n12]) && summable(&[v2])) {
280                in_range = false;
281                return 0;
282            }
283            expansion::grow(&x12[..n12], v2, out)
284        });
285        let z = written(scratch, mark, z, in_range)?;
286        Ok(ExactCoords {
287            x,
288            y,
289            z,
290            w: Weight::One,
291        })
292    }
293}
294
295impl Point for Lnc {
296    #[inline]
297    fn approx(&self) -> &[f64; 3] {
298        &self.approx
299    }
300
301    #[inline]
302    fn err(&self) -> &[f64; 3] {
303        &self.err
304    }
305
306    #[inline]
307    fn is_explicit(&self) -> bool {
308        false
309    }
310
311    fn exact(&self, scratch: &mut Scratch<'_>) -> Result<ExactCoords, ExactError> {
312        Self::write_exact(scratch, [self.p, self.q], [self.t])
313    }
314}
315
316/// A point of any of the types of this module.
317#[derive(Clone, Copy, Debug, PartialEq)]
318pub enum AnyPoint {
319    /// A point of type [`Explicit`].
320    Explicit(Explicit),
321    /// A point of type [`Lnc`].
322    Lnc(Lnc),
323}
324
325impl Point for AnyPoint {
326    #[inline]
327    fn approx(&self) -> &[f64; 3] {
328        match self {
329            AnyPoint::Explicit(p) => p.approx(),
330            AnyPoint::Lnc(p) => p.approx(),
331        }
332    }
333
334    #[inline]
335    fn err(&self) -> &[f64; 3] {
336        match self {
337            AnyPoint::Explicit(p) => p.err(),
338            AnyPoint::Lnc(p) => p.err(),
339        }
340    }
341
342    #[inline]
343    fn is_explicit(&self) -> bool {
344        match self {
345            AnyPoint::Explicit(p) => p.is_explicit(),
346            AnyPoint::Lnc(p) => p.is_explicit(),
347        }
348    }
349
350    #[inline]
351    fn exact(&self, scratch: &mut Scratch<'_>) -> Result<ExactCoords, ExactError> {
352        match self {
353            AnyPoint::Explicit(p) => p.exact(scratch),
354            AnyPoint::Lnc(p) => p.exact(scratch),
355        }
356    }
357}
358
359impl From<Explicit> for AnyPoint {
360    #[inline]
361    fn from(p: Explicit) -> Self {
362        AnyPoint::Explicit(p)
363    }
364}
365
366impl From<Lnc> for AnyPoint {
367    #[inline]
368    fn from(p: Lnc) -> Self {
369        AnyPoint::Lnc(p)
370    }
371}
372
373/// `2^-53`, the unit roundoff of `f64`.
374const UNIT: f64 = 1.1102230246251565e-16;
375
376/// `2^-1074`, the smallest positive `f64`.
377const TINY: f64 = 5e-324;
378
379/// `fl(fl(a * b) + 2^-1074)` for non-negative `a` and `b`, or zero if either is
380/// zero: at least `(1 - 2^-53)^2 * a * b`, also when `a * b` underflows, since
381/// adding `2^-1074` to a value below `2^-1022` is exact.
382#[inline]
383fn covered_product(a: f64, b: f64) -> f64 {
384    if a == 0.0 || b == 0.0 {
385        0.0
386    } else {
387        a * b + TINY
388    }
389}
390
391/// The bound on the rounding error of the sum or difference `s` of `x` and `y`:
392/// the covered product of `2^-53` and `|s|`, or zero if an operand is zero,
393/// when the operation is exact.
394#[inline]
395fn sum_rounding(x: f64, y: f64, s: f64) -> f64 {
396    if x != 0.0 && y != 0.0 {
397        covered_product(UNIT, s.abs())
398    } else {
399        0.0
400    }
401}
402
403/// The bound on the rounding error of the product `p` of `x` and `y`: the
404/// covered product of `2^-53` and `|p|`, plus `2^-1074` for underflow, or zero
405/// if an operand is zero, when the operation is exact. The test is on the
406/// operands: a product that underflows is zero without being exact.
407#[inline]
408fn product_rounding(x: f64, y: f64, p: f64) -> f64 {
409    if x != 0.0 && y != 0.0 {
410        covered_product(UNIT, p.abs()) + TINY
411    } else {
412        0.0
413    }
414}
415
416/// The error bound stored for a coordinate whose value is `approx` and whose
417/// computed error is `err`: `err * slack` rounded up, which is at least the
418/// error since the slack factor covers the roundings of its computation; zero
419/// when `err` is zero; and `+∞` when `approx` or the bound is not finite.
420#[inline]
421fn finish(approx: f64, err: f64, slack: f64) -> f64 {
422    let bound = err * slack;
423    if !(approx.is_finite() && bound.is_finite()) {
424        f64::INFINITY
425    } else if err == 0.0 {
426        0.0
427    } else {
428        bound.next_up()
429    }
430}
431
432/// `2^1020`: the expansions of the core require every operand of a sum to be
433/// below it in magnitude.
434const SUM_OPERAND_LIMIT: f64 = 1.1235582092889474e307;
435
436/// Whether every component of `e` is below `SUM_OPERAND_LIMIT` in magnitude, as
437/// the components of an operand of a sum of expansions must be.
438#[inline]
439fn summable(e: &[f64]) -> bool {
440    e.iter().all(|x| x.abs() < SUM_OPERAND_LIMIT)
441}
442
443/// Whether every component of `e` lies in the product domain of the core, zero
444/// or a magnitude from `2^-480` to `2^480`, as the components of an operand of
445/// an exact product must.
446#[inline]
447fn multipliable(e: &[f64]) -> bool {
448    e.iter().all(|&x| eft::in_product_domain(x))
449}
450
451/// The span of a coordinate just written, or the error of `write_exact` once
452/// the arena is reset to `mark`: [`ExactError::OutOfRange`] when `in_range` is
453/// false.
454#[inline]
455fn written(
456    scratch: &mut Scratch<'_>,
457    mark: Mark,
458    span: Result<Span, ScratchExhausted>,
459    in_range: bool,
460) -> Result<Span, ExactError> {
461    match span {
462        Ok(span) if in_range => Ok(span),
463        Ok(_) => {
464            scratch.reset(mark);
465            Err(ExactError::OutOfRange)
466        }
467        Err(e) => {
468            scratch.reset(mark);
469            Err(e.into())
470        }
471    }
472}
implicit_predicates/generated/points.rs

implicit_predicates/generated/
points.rs
