diff --git a/crates/geom/src/homography.rs b/crates/geom/src/homography.rs
index 6a49cc5d..b27ecef2 100644
--- a/crates/geom/src/homography.rs
+++ b/crates/geom/src/homography.rs
@@ -2,7 +2,7 @@
 
 use serde::{Deserialize, Serialize};
 
-use crate::{Affine, Point};
+use crate::{Affine, Point, Real};
 
 /// Row-major 3×3 matrix mapping homogeneous points: p' = H·p.
 #[derive(Clone, Copy, Debug, PartialEq, Serialize, Deserialize)]
@@ -23,10 +23,15 @@ impl Homography {
     }
 
     pub fn apply(&self, p: Point) -> Point {
+        let (x, y) = self.apply_real(p.x, p.y);
+        Point::new(x, y)
+    }
+
+    /// [`Homography::apply`] for any [`Real`] (e.g. [`Interval`](crate::Interval) bounds over many points).
+    pub fn apply_real<T: Real>(&self, x: T, y: T) -> (T, T) {
         let m = &self.0;
-        let w = m[6] * p.x + m[7] * p.y + m[8];
-        let w = if w.abs() < 1e-300 { 1e-300 } else { w };
-        Point::new((m[0] * p.x + m[1] * p.y + m[2]) / w, (m[3] * p.x + m[4] * p.y + m[5]) / w)
+        let w = (x * m[6] + y * m[7] + m[8]).clamp_tiny(1e-300);
+        ((x * m[0] + y * m[1] + m[2]) / w, (x * m[3] + y * m[4] + m[5]) / w)
     }
 
     pub fn mul(&self, o: &Homography) -> Homography {
diff --git a/crates/geom/src/lib.rs b/crates/geom/src/lib.rs
index 327c60fc..6de121b0 100644
--- a/crates/geom/src/lib.rs
+++ b/crates/geom/src/lib.rs
@@ -7,9 +7,11 @@
 
 mod crop;
 mod homography;
+mod real;
 
 pub use crop::{CropGeometry, crop_fit_angle, max_inscribed_scale};
 pub use homography::Homography;
+pub use real::{Interval, Real};
 
 use serde::{Deserialize, Serialize};
 
@@ -213,8 +215,13 @@ impl Affine {
         Affine::translate(p.to_vec()) * Affine::rotate(a) * Affine::translate(-p.to_vec())
     }
     pub fn apply(&self, p: Point) -> Point {
+        let (x, y) = self.apply_real(p.x, p.y);
+        Point::new(x, y)
+    }
+    /// [`Affine::apply`] for any [`Real`] (e.g. [`Interval`] bounds over many points).
+    pub fn apply_real<T: Real>(&self, x: T, y: T) -> (T, T) {
         let [a, b, c, d, e, f] = self.0;
-        Point::new(a * p.x + c * p.y + e, b * p.x + d * p.y + f)
+        (x * a + y * c + e, x * b + y * d + f)
     }
     pub fn apply_vec(&self, v: Vec2) -> Vec2 {
         let [a, b, c, d, ..] = self.0;
diff --git a/crates/geom/src/real.rs b/crates/geom/src/real.rs
new file mode 100644
index 00000000..07da59f3
--- /dev/null
+++ b/crates/geom/src/real.rs
@@ -0,0 +1,229 @@
+//! Numbers the coordinate formulas are generic over: `f64` (the reference evaluation) and [`Interval`]
+//! (bounds on the `f64` results over a whole range of inputs, e.g. a block of output pixels).
+//!
+//! Formulas written once over [`Real`] serve both. For `f64` every operator is the plain IEEE operation, so a
+//! generic formula evaluated with `f64` is bit-identical to the same formula written with `f64` directly, as long
+//! as it keeps the operation order (`a * k` instead of `k * a` is fine: IEEE `+` and `×` are commutative).
+
+use std::ops::{Add, Div, Mul, Sub};
+
+/// A scalar the geometry formulas can be evaluated with. Mixed operations take a parameter (`f64`) on the right.
+pub trait Real:
+    Copy
+    + Add<Output = Self>
+    + Sub<Output = Self>
+    + Mul<Output = Self>
+    + Div<Output = Self>
+    + Add<f64, Output = Self>
+    + Sub<f64, Output = Self>
+    + Mul<f64, Output = Self>
+    + Div<f64, Output = Self>
+{
+    /// `if |v| < tiny { tiny } else { v }`: the guard of a denominator against zero ([`crate::Homography::apply`]).
+    fn clamp_tiny(self, tiny: f64) -> Self;
+}
+
+impl Real for f64 {
+    fn clamp_tiny(self, tiny: f64) -> f64 {
+        if self.abs() < tiny { tiny } else { self }
+    }
+}
+
+/// A closed interval `[lo, hi]` with outward rounding, or *undecided* (both NaN).
+///
+/// Invariant: if the `f64` operands of an operation lie in the operand intervals, the `f64` result (rounded to
+/// nearest) lies in the result interval. Each bound is computed in `f64` and then moved one ulp outward
+/// (`next_down` / `next_up`); since rounding to nearest is monotonic and errs by less than one ulp, the moved
+/// bound is on the far side of every rounded result. By induction a formula evaluated with intervals bounds the
+/// *`f64` evaluation* of the same formula (not only the exact real value) for every input in the input intervals.
+///
+/// A result that can't be bounded becomes undecided and stays so: division by an interval that contains zero,
+/// [`Real::clamp_tiny`] when some value could be clamped, and any infinite or NaN bound. Undecided intervals fail
+/// every comparison of their bounds, so callers must test the bounds with conditions that are false for NaN.
+#[derive(Clone, Copy, Debug)]
+pub struct Interval {
+    lo: f64,
+    hi: f64,
+}
+
+impl Interval {
+    const UNDECIDED: Interval = Interval { lo: f64::NAN, hi: f64::NAN };
+
+    /// `[lo, hi]` exactly (the bounds are `f64` values already); undecided unless finite with `lo <= hi`.
+    pub fn new(lo: f64, hi: f64) -> Interval {
+        if lo.is_finite() && hi.is_finite() && lo <= hi { Interval { lo, hi } } else { Interval::UNDECIDED }
+    }
+
+    /// The lower bound (NaN if undecided).
+    pub fn lo(self) -> f64 {
+        self.lo
+    }
+
+    /// The upper bound (NaN if undecided).
+    pub fn hi(self) -> f64 {
+        self.hi
+    }
+
+    /// Whether the bounds are usable (not undecided).
+    pub fn is_decided(self) -> bool {
+        self.lo <= self.hi
+    }
+
+    /// Bounds computed in round-to-nearest, moved one ulp outward.
+    fn outward(lo: f64, hi: f64) -> Interval {
+        Interval::new(lo.next_down(), hi.next_up())
+    }
+
+    /// The hull of four products or quotients, moved outward. NaN in any of them makes it undecided.
+    fn hull(v: [f64; 4]) -> Interval {
+        if v.iter().any(|v| !v.is_finite()) {
+            return Interval::UNDECIDED;
+        }
+        let lo = v.iter().copied().fold(f64::INFINITY, f64::min);
+        let hi = v.iter().copied().fold(f64::NEG_INFINITY, f64::max);
+        Interval::outward(lo, hi)
+    }
+}
+
+impl From<f64> for Interval {
+    fn from(v: f64) -> Interval {
+        Interval::new(v, v)
+    }
+}
+
+impl Add for Interval {
+    type Output = Interval;
+    fn add(self, o: Interval) -> Interval {
+        Interval::outward(self.lo + o.lo, self.hi + o.hi)
+    }
+}
+
+impl Sub for Interval {
+    type Output = Interval;
+    fn sub(self, o: Interval) -> Interval {
+        Interval::outward(self.lo - o.hi, self.hi - o.lo)
+    }
+}
+
+impl Mul for Interval {
+    type Output = Interval;
+    fn mul(self, o: Interval) -> Interval {
+        Interval::hull([self.lo * o.lo, self.lo * o.hi, self.hi * o.lo, self.hi * o.hi])
+    }
+}
+
+impl Div for Interval {
+    type Output = Interval;
+    /// Undecided when the divisor may be zero (or is undecided).
+    fn div(self, o: Interval) -> Interval {
+        if !(o.lo > 0.0 || o.hi < 0.0) {
+            return Interval::UNDECIDED;
+        }
+        Interval::hull([self.lo / o.lo, self.lo / o.hi, self.hi / o.lo, self.hi / o.hi])
+    }
+}
+
+impl Add<f64> for Interval {
+    type Output = Interval;
+    fn add(self, k: f64) -> Interval {
+        Interval::outward(self.lo + k, self.hi + k)
+    }
+}
+
+impl Sub<f64> for Interval {
+    type Output = Interval;
+    fn sub(self, k: f64) -> Interval {
+        Interval::outward(self.lo - k, self.hi - k)
+    }
+}
+
+impl Mul<f64> for Interval {
+    type Output = Interval;
+    fn mul(self, k: f64) -> Interval {
+        let (a, b) = (self.lo * k, self.hi * k);
+        if k >= 0.0 { Interval::outward(a, b) } else { Interval::outward(b, a) }
+    }
+}
+
+impl Div<f64> for Interval {
+    type Output = Interval;
+    /// Undecided when `k` is zero (or NaN).
+    fn div(self, k: f64) -> Interval {
+        let (a, b) = (self.lo / k, self.hi / k);
+        if k > 0.0 {
+            Interval::outward(a, b)
+        } else if k < 0.0 {
+            Interval::outward(b, a)
+        } else {
+            Interval::UNDECIDED
+        }
+    }
+}
+
+impl Real for Interval {
+    /// Unchanged when no value in the interval would be clamped, else undecided.
+    fn clamp_tiny(self, tiny: f64) -> Interval {
+        if self.lo >= tiny || self.hi <= -tiny { self } else { Interval::UNDECIDED }
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use proptest::prelude::*;
+
+    fn contains(i: Interval, v: f64) -> bool {
+        i.lo <= v && v <= i.hi
+    }
+
+    /// A value inside `[lo, hi]` at fraction `t`.
+    fn pick(i: Interval, t: f64) -> f64 {
+        (i.lo + (i.hi - i.lo) * t).clamp(i.lo, i.hi)
+    }
+
+    proptest! {
+        /// Every `f64` result of `a ∘ b` with `a ∈ A`, `b ∈ B` lies in `A ∘ B` (or that is undecided).
+        #[test]
+        fn operations_enclose_the_f64_results(
+            a0 in -1e6f64..1e6, aw in 0.0f64..1e3, b0 in -1e3f64..1e3, bw in 0.0f64..10.0,
+            ta in 0.0f64..=1.0, tb in 0.0f64..=1.0, k in -1e3f64..1e3,
+        ) {
+            let (ia, ib) = (Interval::new(a0, a0 + aw), Interval::new(b0, b0 + bw));
+            let (a, b) = (pick(ia, ta), pick(ib, tb));
+            prop_assert!(contains(ia + ib, a + b));
+            prop_assert!(contains(ia - ib, a - b));
+            prop_assert!(contains(ia * ib, a * b));
+            let q = ia / ib;
+            prop_assert!(!q.is_decided() || contains(q, a / b));
+            prop_assert!(contains(ia + k, a + k));
+            prop_assert!(contains(ia - k, a - k));
+            prop_assert!(contains(ia * k, a * k));
+            prop_assert!(k == 0.0 || contains(ia / k, a / k));
+            let c = ib.clamp_tiny(1.0);
+            prop_assert!(!c.is_decided() || contains(c, b.clamp_tiny(1.0)));
+        }
+    }
+
+    #[test]
+    fn unbounded_results_are_undecided() {
+        let i = Interval::new(-1.0, 2.0);
+        assert!(!(i / Interval::new(-1e-9, 1.0)).is_decided());
+        assert!(!(i / Interval::new(0.0, 1.0)).is_decided());
+        assert!((i / Interval::new(1e-9, 1.0)).is_decided());
+        assert!(!(i / 0.0).is_decided());
+        assert!(!(Interval::new(f64::MAX, f64::MAX) * 2.0).is_decided());
+        assert!(!Interval::new(2.0, 1.0).is_decided());
+        assert!(!Interval::new(f64::NAN, 1.0).is_decided());
+        // undecided propagates through everything
+        let u = Interval::UNDECIDED;
+        for r in [u + i, i - u, u * i, i / u, u + 1.0, u * 0.0, u.clamp_tiny(1e-300)] {
+            assert!(!r.is_decided());
+            assert!(r.lo.is_nan() && r.hi.is_nan(), "undecided bounds are NaN, which fails every comparison");
+        }
+        // the denominator guard: unchanged only when no value is clamped
+        assert!(Interval::new(1e-300, 1.0).clamp_tiny(1e-300).is_decided());
+        assert!(Interval::new(-1.0, -1e-300).clamp_tiny(1e-300).is_decided());
+        assert!(!Interval::new(-1.0, 1.0).clamp_tiny(1e-300).is_decided());
+        assert!(!Interval::new(5e-301, 1.0).clamp_tiny(1e-300).is_decided());
+    }
+}
diff --git a/crates/pipeline/src/geometry.rs b/crates/pipeline/src/geometry.rs
index 56b158e2..af45526b 100644
--- a/crates/pipeline/src/geometry.rs
+++ b/crates/pipeline/src/geometry.rs
@@ -282,13 +282,10 @@ fn sample_warped(base: &Rgb32f, sx: f64, sy: f64, wp: &Warp, o2t: Affine, w: usi
     let mut out = Rgb32f::new(w, h);
     par_rows(&mut out.data, w, |y, row| {
         for (x, px) in row.iter_mut().enumerate() {
-            let t = o2t.apply(Point::new(x as f64 + 0.5, y as f64 + 0.5));
-            let c = wp.to_corrected(t);
-            let g = wp.corrected_to_source(c, 1);
-            if !wp.inside(g) {
+            let Some((c, g)) = wp.frame(&o2t, x, y) else {
                 *px = BLANK;
                 continue;
-            }
+            };
             let mut v = base.sample_bilinear((g.x * sx) as f32, (g.y * sy) as f32);
             if per_channel {
                 for ch in [0usize, 2] {
diff --git a/crates/pipeline/src/optics.rs b/crates/pipeline/src/optics.rs
index 7c2cec27..ad789fc8 100644
--- a/crates/pipeline/src/optics.rs
+++ b/crates/pipeline/src/optics.rs
@@ -17,7 +17,7 @@ use std::sync::Mutex;
 
 use lightcraft_color::perceptual::{oklab_from_2020, oklab_to_2020};
 use lightcraft_develop::{DevelopSettings, EmbeddedLens, EmbeddedVignette, EmbeddedWarp};
-use lightcraft_geom::{Homography, Orientation, Point};
+use lightcraft_geom::{Affine, Homography, Interval, Orientation, Point, Real};
 use lightcraft_raster::{Plane, Rgb32f};
 
 use crate::for_rows;
@@ -122,12 +122,18 @@ impl Warp {
 
     /// Transformed px → lens-corrected px (inverse perspective).
     pub fn to_corrected(&self, p: Point) -> Point {
+        let (x, y) = self.to_corrected_real(p.x, p.y);
+        Point::new(x, y)
+    }
+
+    /// [`Warp::to_corrected`] for any [`Real`].
+    fn to_corrected_real<T: Real>(&self, x: T, y: T) -> (T, T) {
         if self.persp == Homography::IDENTITY {
-            return p;
+            return (x, y);
         }
         let (c, l) = (self.centre(), self.half_long());
-        let q = self.persp_inv.apply(Point::new((p.x - c.x) / l, (p.y - c.y) / l));
-        Point::new(c.x + q.x * l, c.y + q.y * l)
+        let (qx, qy) = self.persp_inv.apply_real((x - c.x) / l, (y - c.y) / l);
+        (qx * l + c.x, qy * l + c.y)
     }
 
     /// Lens-corrected px → transformed px (forward perspective).
@@ -142,26 +148,32 @@ impl Warp {
 
     /// Lens-corrected px → source px for colour plane `ch` (0 = R, 1 = G, 2 = B).
     pub fn corrected_to_source(&self, p: Point, ch: usize) -> Point {
+        let (x, y) = self.corrected_to_source_real(p.x, p.y, ch);
+        Point::new(x, y)
+    }
+
+    /// [`Warp::corrected_to_source`] for any [`Real`].
+    fn corrected_to_source_real<T: Real>(&self, x: T, y: T, ch: usize) -> (T, T) {
         let c = self.centre();
         let hd = (self.w.hypot(self.h) / 2.0).max(1e-9);
-        let mut q = p;
+        let (mut qx, mut qy) = (x, y);
         if self.k1 != 0.0 {
-            let (dx, dy) = (q.x - c.x, q.y - c.y);
+            let (dx, dy) = (qx - c.x, qy - c.y);
             let r2 = (dx * dx + dy * dy) / (hd * hd);
-            let f = 1.0 + self.k1 * r2;
-            q = Point::new(c.x + dx * f, c.y + dy * f);
+            let f = r2 * self.k1 + 1.0;
+            (qx, qy) = (dx * f + c.x, dy * f + c.y);
         }
         if self.lens_dist != 0.0
             && let Some(wp) = self.lens.and_then(|l| l.warp)
         {
-            let s = embedded_warp(&wp, q, ch, self.w, self.h);
-            q = Point::new(q.x + (s.x - q.x) * self.lens_dist, q.y + (s.y - q.y) * self.lens_dist);
+            let (sx, sy) = embedded_warp_real(&wp, qx, qy, ch, self.w, self.h);
+            (qx, qy) = (qx + (sx - qx) * self.lens_dist, qy + (sy - qy) * self.lens_dist);
         }
         let k = self.ca[ch];
         if k != 0.0 {
-            q = Point::new(c.x + (q.x - c.x) * (1.0 + k), c.y + (q.y - c.y) * (1.0 + k));
+            (qx, qy) = ((qx - c.x) * (1.0 + k) + c.x, (qy - c.y) * (1.0 + k) + c.y);
         }
-        q
+        (qx, qy)
     }
 
     /// Transformed px → source px for plane `ch`.
@@ -185,24 +197,86 @@ impl Warp {
         g as f32
     }
 
-    /// Whether a source position lies on the image (with a half-pixel margin).
-    pub fn inside(&self, s: Point) -> bool {
-        s.x >= -0.5 && s.y >= -0.5 && s.x <= self.w + 0.5 && s.y <= self.h + 0.5
+    /// Whether a source position lies on the image (within [`edges`](Self::edges)).
+    fn inside(&self, s: Point) -> bool {
+        let (lo, hx, hy) = self.edges();
+        s.x >= lo && s.y >= lo && s.x <= hx && s.y <= hy
+    }
+
+    /// The source area that frames output pixels, `lo..=hx × lo..=hy`: the image with a half-pixel margin.
+    fn edges(&self) -> (f64, f64, f64) {
+        (-0.5, self.w + 0.5, self.h + 0.5)
+    }
+
+    /// The framing chain of output position `(x, y)` (`o2t`: output px → transformed px): its lens-corrected
+    /// position and its green-plane source position.
+    fn framing<T: Real>(&self, o2t: &Affine, x: T, y: T) -> ((T, T), (T, T)) {
+        let (tx, ty) = o2t.apply_real(x, y);
+        let (cx, cy) = self.to_corrected_real(tx, ty);
+        ((cx, cy), self.corrected_to_source_real(cx, cy, 1))
+    }
+
+    /// The reference framing of output pixel `(x, y)`: its lens-corrected and green-plane source positions if
+    /// the latter is [`inside`](Self::inside) the image, else `None` (the pixel is blank canvas).
+    pub fn frame(&self, o2t: &Affine, x: usize, y: usize) -> Option<(Point, Point)> {
+        let ((cx, cy), (gx, gy)) = self.framing(o2t, x as f64 + 0.5, y as f64 + 0.5);
+        let g = Point::new(gx, gy);
+        self.inside(g).then(|| (Point::new(cx, cy), g))
+    }
+
+    /// Whether output pixel `(x, y)` shows the image ([`frame`](Self::frame) is `Some`).
+    pub fn covers(&self, o2t: &Affine, x: usize, y: usize) -> bool {
+        self.frame(o2t, x, y).is_some()
+    }
+
+    /// [`covers`](Self::covers) for every pixel of the output block `x0..x1 × y0..y1` at once, when that is
+    /// certain: `Some(true)` if all are covered, `Some(false)` if none is, `None` if the block may straddle the
+    /// image edge (or can't be bounded) and each pixel must be decided alone.
+    ///
+    /// Why the answer equals the per-pixel decision: the block's pixel centres `x + 0.5` are exact `f64` values in
+    /// `[x0 + 0.5, x1 − 0.5]` (likewise y). [`framing`](Self::framing) evaluated with [`Interval`]s runs the very
+    /// same operations as its `f64` evaluation in [`frame`](Self::frame), each rounded outward, so its result
+    /// bounds the `f64` source position of every pixel of the block — the rounded values, not only the exact
+    /// ones (see [`Interval`]). Divisions by a range that may contain zero, the homography's denominator clamp
+    /// and overflow make the bounds undecided (NaN), and NaN fails every test below. So when the bounds lie
+    /// within [`edges`](Self::edges) every pixel's position does, and when they lie beyond one edge every
+    /// pixel's position does: the same `f64` comparisons [`inside`](Self::inside) makes, no margin needed.
+    pub fn block_coverage(&self, o2t: &Affine, x0: usize, x1: usize, y0: usize, y1: usize) -> Option<bool> {
+        if x1 <= x0 || y1 <= y0 {
+            return None;
+        }
+        let px = Interval::new(x0 as f64 + 0.5, (x1 - 1) as f64 + 0.5);
+        let py = Interval::new(y0 as f64 + 0.5, (y1 - 1) as f64 + 0.5);
+        let (_, (sx, sy)) = self.framing(o2t, px, py);
+        let (lo, hx, hy) = self.edges();
+        if sx.lo() >= lo && sy.lo() >= lo && sx.hi() <= hx && sy.hi() <= hy {
+            Some(true)
+        } else if sx.hi() < lo || sy.hi() < lo || sx.lo() > hx || sy.lo() > hy {
+            Some(false)
+        } else {
+            None
+        }
     }
 }
 
 /// DNG `WarpRectilinear` source position for corrected position `p` (px) and plane `ch`.
 pub fn embedded_warp(wp: &EmbeddedWarp, p: Point, ch: usize, w: f64, h: f64) -> Point {
+    let (x, y) = embedded_warp_real(wp, p.x, p.y, ch, w, h);
+    Point::new(x, y)
+}
+
+/// [`embedded_warp`] for any [`Real`].
+fn embedded_warp_real<T: Real>(wp: &EmbeddedWarp, x: T, y: T, ch: usize, w: f64, h: f64) -> (T, T) {
     let (cx, cy) = (wp.center.x * w, wp.center.y * h);
     let m = (wp.radius * w.max(h)).max(1e-9);
     let [kr0, kr1, kr2, kr3, kt0, kt1] = wp.planes[ch.min(2)];
-    let dx = (p.x - cx) / m;
-    let dy = (p.y - cy) / m;
+    let dx = (x - cx) / m;
+    let dy = (y - cy) / m;
     let r2 = dx * dx + dy * dy;
-    let f = kr0 + r2 * (kr1 + r2 * (kr2 + r2 * kr3));
-    let sx = dx * f + kt0 * 2.0 * dx * dy + kt1 * (r2 + 2.0 * dx * dx);
-    let sy = dy * f + kt1 * 2.0 * dx * dy + kt0 * (r2 + 2.0 * dy * dy);
-    Point::new(cx + m * sx, cy + m * sy)
+    let f = (((r2 * kr3 + kr2) * r2 + kr1) * r2) + kr0;
+    let sx = dx * f + dx * (kt0 * 2.0) * dy + (r2 + dx * 2.0 * dx) * kt1;
+    let sy = dy * f + dx * (kt1 * 2.0) * dy + (r2 + dy * 2.0 * dy) * kt0;
+    (sx * m + cx, sy * m + cy)
 }
 
 /// DNG `FixVignetteRadial` gain at source position `p` (px).
@@ -514,6 +588,224 @@ fn min_filter(p: &Plane, r: usize, max: bool) -> Plane {
 mod tests {
     use super::*;
 
+    /// splitmix64: deterministic pseudo-random numbers for the warp tests.
+    struct Rng(u64);
+
+    impl Rng {
+        fn next(&mut self) -> u64 {
+            self.0 = self.0.wrapping_add(0x9E37_79B9_7F4A_7C15);
+            let mut z = self.0;
+            z = (z ^ (z >> 30)).wrapping_mul(0xBF58_476D_1CE4_E5B9);
+            z = (z ^ (z >> 27)).wrapping_mul(0x94D0_49BB_1331_11EB);
+            z ^ (z >> 31)
+        }
+        /// Uniform in `lo..hi`.
+        fn range(&mut self, lo: f64, hi: f64) -> f64 {
+            lo + (hi - lo) * ((self.next() >> 11) as f64 / (1u64 << 53) as f64)
+        }
+        /// `range(lo, hi)` or, one time in four, exactly 0 (to take the "no correction" branches).
+        fn maybe(&mut self, lo: f64, hi: f64) -> f64 {
+            if self.next().is_multiple_of(4) { 0.0 } else { self.range(lo, hi) }
+        }
+    }
+
+    /// A pseudo-random warp: perspective, manual distortion, CA and an embedded lens warp, each sometimes off.
+    /// Sizes are Pythagorean (exact `hypot`), so results don't depend on the platform's libm.
+    fn random_warp(r: &mut Rng) -> Warp {
+        const SIZES: [(f64, f64); 6] = [(4000.0, 3000.0), (3000.0, 4000.0), (1200.0, 900.0), (1200.0, 500.0), (56.0, 33.0), (6000.0, 4500.0)];
+        let (w, h) = SIZES[(r.next() % SIZES.len() as u64) as usize];
+        let mut wp = Warp::identity(w, h);
+        if !r.next().is_multiple_of(4) {
+            let m = [
+                1.0 + r.range(-0.2, 0.2),
+                r.range(-0.2, 0.2),
+                r.range(-0.1, 0.1),
+                r.range(-0.2, 0.2),
+                1.0 + r.range(-0.2, 0.2),
+                r.range(-0.1, 0.1),
+                r.range(-0.4, 0.4),
+                r.range(-0.4, 0.4),
+                1.0 + r.range(-0.1, 0.1),
+            ];
+            wp.persp_inv = Homography(m);
+            wp.persp = wp.persp_inv.inverse().unwrap_or(Homography::IDENTITY);
+        }
+        wp.k1 = r.maybe(-0.2, 0.2);
+        wp.ca = [r.maybe(-0.003, 0.003), r.maybe(-0.003, 0.003), r.maybe(-0.003, 0.003)];
+        if !r.next().is_multiple_of(4) {
+            let mut plane = || {
+                [
+                    1.0 + r.range(-0.05, 0.05),
+                    r.range(-0.08, 0.08),
+                    r.range(-0.05, 0.05),
+                    r.range(-0.02, 0.02),
+                    r.range(-0.003, 0.003),
+                    r.range(-0.003, 0.003),
+                ]
+            };
+            let planes = [plane(), plane(), plane()];
+            let warp = EmbeddedWarp { planes, center: Point::new(r.range(0.4, 0.6), r.range(0.4, 0.6)), radius: r.range(0.5, 0.8) };
+            wp.lens = Some(EmbeddedLens { warp: Some(warp), vignette: None });
+            wp.lens_dist = r.maybe(0.0, 2.0);
+        }
+        wp
+    }
+
+    /// FNV-1a over the bits of the f64 warp results (`Affine::apply`, `Homography::apply`, `to_corrected`,
+    /// `from_corrected`, `corrected_to_source`, `to_source`, `embedded_warp`) for pseudo-random warps and points.
+    /// The expected value was captured before these formulas became generic over `Real`: the f64 path must keep
+    /// the same operations in the same order, bit for bit. To recompute it after changing the inputs, run this
+    /// test (printing `hash`) on a scratch copy of 419261c (thumb-split, pre-generic formulas) and paste it here.
+    #[test]
+    fn warp_f64_results_are_bit_identical_to_the_pre_generic_formulas() {
+        let mut hash = 0xCBF2_9CE4_8422_2325u64;
+        let mut eat = |p: Point| {
+            for v in [p.x, p.y] {
+                for b in v.to_bits().to_le_bytes() {
+                    hash = (hash ^ b as u64).wrapping_mul(0x0000_0100_0000_01B3);
+                }
+            }
+        };
+        let mut r = Rng(42);
+        for _ in 0..2000 {
+            let wp = random_warp(&mut r);
+            let a = lightcraft_geom::Affine([
+                r.range(-2.0, 2.0),
+                r.range(-2.0, 2.0),
+                r.range(-2.0, 2.0),
+                r.range(-2.0, 2.0),
+                r.range(-500.0, 500.0),
+                r.range(-500.0, 500.0),
+            ]);
+            for _ in 0..40 {
+                let p = Point::new(r.range(-0.2, 1.2) * wp.w, r.range(-0.2, 1.2) * wp.h);
+                eat(a.apply(p));
+                eat(wp.persp_inv.apply(Point::new(p.x / wp.w - 0.5, p.y / wp.h - 0.5)));
+                eat(wp.to_corrected(p));
+                eat(wp.from_corrected(p));
+                for ch in 0..3 {
+                    eat(wp.corrected_to_source(p, ch));
+                    eat(wp.to_source(p, ch));
+                }
+                if let Some(ew) = wp.lens.and_then(|l| l.warp) {
+                    for ch in 0..4 {
+                        eat(embedded_warp(&ew, p, ch, wp.w, wp.h));
+                    }
+                }
+            }
+        }
+        // `Homography::apply`'s near-zero denominator clamp, including subnormal and signed-zero denominators and
+        // ±1e-300 exactly (not clamped: `<`, not `<=`)
+        let hz = Homography([1.0, 0.5, 0.25, -0.5, 1.0, 0.125, 1.0, 0.0, 0.0]);
+        for x in [0.0, -0.0, 1e-310, -1e-310, 1e-300, -1e-300, 1e-299, -1e-299, 0.5] {
+            eat(hz.apply(Point::new(x, 0.75)));
+        }
+        assert_eq!(hash, 0x3FA5_7473_D795_9E32, "f64 warp results changed");
+    }
+
+    /// The interval evaluation of `framing` bounds the f64 evaluation of every point of the box.
+    #[test]
+    fn interval_framing_encloses_every_f64_result() {
+        let mut r = Rng(7);
+        let (mut decided, mut total) = (0, 0);
+        for _ in 0..500 {
+            let wp = random_warp(&mut r);
+            let a = lightcraft_geom::Affine([
+                r.range(0.5, 2.0),
+                r.range(-0.2, 0.2),
+                r.range(-0.2, 0.2),
+                r.range(0.5, 2.0),
+                r.range(-300.0, 300.0),
+                r.range(-300.0, 300.0),
+            ]);
+            for _ in 0..20 {
+                let (x0, y0) = (r.range(-0.2, 1.1) * wp.w, r.range(-0.2, 1.1) * wp.h);
+                let (bw, bh) = (r.range(0.0, 40.0), r.range(0.0, 40.0));
+                let (ix, iy) = (Interval::new(x0, x0 + bw), Interval::new(y0, y0 + bh));
+                let (_, (sx, sy)) = wp.framing(&a, ix, iy);
+                total += 1;
+                if !(sx.is_decided() && sy.is_decided()) {
+                    continue;
+                }
+                decided += 1;
+                for i in 0..=6 {
+                    for j in 0..=6 {
+                        let (x, y) = if i == 6 && j == 6 {
+                            (r.range(x0, x0 + bw), r.range(y0, y0 + bh))
+                        } else {
+                            (x0 + bw * i as f64 / 5.0, y0 + bh * j as f64 / 5.0)
+                        };
+                        let (x, y) = (x.clamp(ix.lo(), ix.hi()), y.clamp(iy.lo(), iy.hi()));
+                        let (_, (fx, fy)) = wp.framing(&a, x, y);
+                        assert!(
+                            sx.lo() <= fx && fx <= sx.hi() && sy.lo() <= fy && fy <= sy.hi(),
+                            "({x}, {y}) → ({fx}, {fy}) outside {sx:?} × {sy:?}"
+                        );
+                    }
+                }
+            }
+        }
+        assert!(decided * 10 > total * 9, "{decided}/{total} boxes bounded");
+    }
+
+    #[test]
+    fn block_coverage_agrees_with_every_pixel() {
+        let lens = EmbeddedLens {
+            warp: Some(EmbeddedWarp {
+                planes: [[1.0, -0.03, 0.01, 0.0, 0.001, -0.002], [1.0, -0.028, 0.01, 0.0, 0.001, -0.002], [1.0, -0.026, 0.01, 0.0, 0.001, -0.002]],
+                center: Point::new(0.52, 0.48),
+                radius: 0.6,
+            }),
+            vignette: None,
+        };
+        let mut cases: Vec<DevelopSettings> = Vec::new();
+        let mut s = DevelopSettings::default();
+        s.optics.lens_profile = true;
+        s.optics.ca_red = 30.0;
+        s.geometry.horizontal = -15.0;
+        cases.push(s.clone());
+        s.orientation = Orientation::Rotate270;
+        cases.push(s.clone());
+        let mut s = DevelopSettings::default();
+        s.optics.distortion = -60.0;
+        s.optics.ca_blue = 50.0;
+        s.geometry.vertical = 25.0;
+        s.geometry.rotate = 7.0;
+        cases.push(s.clone());
+        s.optics.distortion = 80.0;
+        s.geometry.horizontal = 30.0;
+        s.orientation = Orientation::Rotate90;
+        cases.push(s);
+        for (i, s) in cases.iter().enumerate() {
+            let f = crate::geometry::Frame::with_lens(900, 600, s, true, Some(&lens));
+            let wp = f.warp.as_ref().expect("warp");
+            let (w, h) = f.fit(700, 700);
+            let o2t = f.out_to_oriented(w, h);
+            let (mut decided, mut tiles, mut blank) = (0, 0, 0);
+            for (tw, th) in [(32, 8), (7, 5), (64, 64)] {
+                for y0 in (0..h).step_by(th) {
+                    for x0 in (0..w).step_by(tw) {
+                        let (x1, y1) = ((x0 + tw).min(w), (y0 + th).min(h));
+                        tiles += 1;
+                        let Some(all) = wp.block_coverage(&o2t, x0, x1, y0, y1) else { continue };
+                        decided += 1;
+                        for y in y0..y1 {
+                            for x in x0..x1 {
+                                assert_eq!(wp.covers(&o2t, x, y), all, "case {i}: pixel {x},{y} in block {x0}..{x1} × {y0}..{y1}");
+                                blank += !all as usize;
+                            }
+                        }
+                    }
+                }
+            }
+            // most blocks are decided wholesale, including blank ones away from the edge
+            assert!(decided * 10 > tiles * 8, "case {i}: {decided}/{tiles} blocks decided");
+            // strong distortion leaves blank margins wide enough for whole blocks
+            assert!(blank > 0 || i < 2, "case {i}: no blank block decided");
+            assert_eq!(wp.block_coverage(&o2t, 5, 5, 0, 1), None, "empty block");
+        }
+    }
+
     /// A grid of dark lines on white with lateral CA: red magnified by `ar`, blue by `ab` (about the centre).
     fn ca_grid(w: usize, h: usize, ar: f64, ab: f64) -> Rgb32f {
         let (cx, cy) = (w as f64 / 2.0, h as f64 / 2.0);