From e365c09141f52255495a2cff8dac1ca1c404547a Mon Sep 17 00:00:00 2001 From: dotCooCoo Date: Mon, 2 Nov 2026 07:10:04 -0800 Subject: [PATCH 1/8] Geometry: a projective map's horizon and invertibility checks do not depend on its scale A homography means the same map at any nonzero uniform scale of its entries, but `apply`, `is_invertible` and `inverse` compared w and the determinant with absolute thresholds. A map stored with small entries refused every point and could not be inverted, and a map with large entries drew points past its horizon. - `apply` refuses a point when w is at most 1e-9 times the largest magnitude in the w row. - A map counts as invertible when its determinant is more than 1e-15 times the largest of the six products the determinant sums, taken on a copy rescaled by a power of two. An inverse whose determinant overflows or underflows is taken on that copy. --- crates/geom/src/projective.rs | 116 +++++++++++++++++++++++----------- docs/target-app-parity.md | 1 + 2 files changed, 80 insertions(+), 37 deletions(-) diff --git a/crates/geom/src/projective.rs b/crates/geom/src/projective.rs index 80758d3ba..a5ae446a5 100644 --- a/crates/geom/src/projective.rs +++ b/crates/geom/src/projective.rs @@ -50,14 +50,16 @@ impl Homography { [m[0][0], m[0][1], m[0][2], m[1][0], m[1][1], m[1][2], m[2][0], m[2][1], m[2][2]] } - /// Map a point; `None` when it lands on or beyond the horizon (w ≤ 0), or when it does not map - /// to a finite page point. + /// Map a point; `None` when it lands on or beyond the horizon (w at most 1e-9 times the largest + /// magnitude in the map's `w` row, so the answer is the same at any scale of the map), or when + /// it does not map to a finite page point. pub fn apply(&self, p: Point) -> Option { let m = &self.0; let x = m[0][0] * p.x + m[0][1] * p.y + m[0][2]; let y = m[1][0] * p.x + m[1][1] * p.y + m[1][2]; let w = m[2][0] * p.x + m[2][1] * p.y + m[2][2]; - if w <= 1e-9 || !x.is_finite() || !y.is_finite() { + let k = m[2].iter().fold(0.0_f64, |k, v| k.max(v.abs())); + if w <= 1e-9 * k || !x.is_finite() || !y.is_finite() { return None; } let q = Point::new(x / w, y / w); @@ -70,16 +72,34 @@ impl Homography { + m[0][2] * (m[1][0] * m[2][1] - m[1][1] * m[2][0]) } - /// Whether the map is invertible: the magnitude of its determinant is above 1e-15, or, when the + /// Whether the map is invertible: [`Self::scaled_det`] gives a determinant, or, when the /// determinant is not finite, [`Self::inverse`] gives one. fn is_invertible(&self) -> bool { - let d = self.det(); - if d.is_finite() { - return d.abs() > 1e-15; + if self.det().is_finite() { + return self.scaled_det().is_some(); } self.inverse().is_some() } + /// [`Self::scaled`] and its determinant; `None` when the magnitude of that determinant is at + /// most 1e-15 times that of the largest of the six products it sums (a test that gives the + /// same answer at any scale of the map). + fn scaled_det(&self) -> Option<(Self, f64)> { + let s = self.scaled()?; + let m = &s.0; + let terms = [ + m[0][0] * m[1][1] * m[2][2], + m[0][0] * m[1][2] * m[2][1], + m[0][1] * m[1][0] * m[2][2], + m[0][1] * m[1][2] * m[2][0], + m[0][2] * m[1][0] * m[2][1], + m[0][2] * m[1][1] * m[2][0], + ]; + let largest = terms.iter().fold(0.0_f64, |k, v| k.max(v.abs())); + let d = s.det(); + (d.abs() > 1e-15 * largest).then_some((s, d)) + } + /// The same map with every entry divided by one power of two, chosen so that the largest of the /// six products its determinant sums is at least 1/8 and below 4. `None` when every one of /// those products is 0 or an entry is not finite. @@ -100,36 +120,14 @@ impl Homography { Some(Self(self.0.map(|row| row.map(|v| ldexp(v, -e))))) } - /// Inverse; `None` when the map is not invertible or its inverse has an entry that is not - /// finite. When the determinant is not finite, this is the inverse of [`Self::scaled`] (the - /// exact inverse multiplied by one power of two, which is the same map), and that copy counts - /// as not invertible when the magnitude of its determinant is at most 1e-15 times that of the largest - /// of the six products the determinant sums. + /// Inverse; `None` when [`Self::scaled_det`] counts the map as not invertible or its inverse + /// has an entry that is not finite. When the determinant is not finite or is below the normal + /// range, this is the inverse of [`Self::scaled`] (the exact inverse multiplied by one power of + /// two, which is the same map). pub fn inverse(&self) -> Option { + let (s, sd) = self.scaled_det()?; let d = self.det(); - let (h, d) = if d.is_finite() { - if d.abs() < 1e-15 { - return None; - } - (*self, d) - } else { - let s = self.scaled()?; - let m = &s.0; - let terms = [ - m[0][0] * m[1][1] * m[2][2], - m[0][0] * m[1][2] * m[2][1], - m[0][1] * m[1][0] * m[2][2], - m[0][1] * m[1][2] * m[2][0], - m[0][2] * m[1][0] * m[2][1], - m[0][2] * m[1][1] * m[2][0], - ]; - let largest = terms.iter().fold(0.0_f64, |k, v| k.max(v.abs())); - let d = s.det(); - if d.abs() <= 1e-15 * largest { - return None; - } - (s, d) - }; + let (h, d) = if d.is_finite() && d.abs() >= f64::MIN_POSITIVE { (*self, d) } else { (s, sd) }; let m = &h.0; let c = |r0: usize, c0: usize, r1: usize, c1: usize| m[r0][c0] * m[r1][c1] - m[r0][c1] * m[r1][c0]; let adj = [ @@ -397,12 +395,56 @@ mod tests { let want = t * far.apply(p).unwrap(); assert!(c.apply(t * p).unwrap().distance(want) <= 1e-9 * want.to_vec2().hypot(), "{c:?}"); } - // A determinant that does not overflow is judged as it is. - assert!(Homography::from_array([1e-6, 0.0, 0.0, 0.0, 1e-6, 0.0, 0.0, 0.0, 1e-6]).is_none()); + // A determinant that does not overflow is judged against the size of the map's entries: + // the identity with every entry 1e-6 is the identity, and two equal rows are not invertible. + let small = [1e-6, 0.0, 0.0, 0.0, 1e-6, 0.0, 0.0, 0.0, 1e-6]; + let q = Homography::from_array(small).unwrap().apply(Point::new(2.0, 3.0)).unwrap(); + assert!(q.distance(Point::new(2.0, 3.0)) < 1e-12); + assert!(Homography::from_array([1e-6, 2e-6, 0.0, 1e-6, 2e-6, 0.0, 0.0, 0.0, 1e-6]).is_none()); assert!(Homography::from_array([f64::NAN; 9]).is_none()); assert!(Homography::from_array([0.0; 9]).is_none()); } + #[test] + fn the_checks_do_not_depend_on_the_scale_of_the_map() { + let times = |h: Homography, f: f64| Homography(h.0.map(|row| row.map(|v| v * f))); + // The same map with its w row's largest entry 1. + let unit = |h: Homography| { + let k = h.0[2].iter().fold(0.0_f64, |k, v| k.max(v.abs())); + Homography(h.0.map(|row| row.map(|v| v / k))) + }; + // The horizon of h() is 0.001x + 0.002y + 1 = 0: two points in front of it, one with w of + // about 5e-10 (too close to it to draw), and two beyond it. + let pts = + [Point::new(30.0, -12.0), Point::new(-999.99, 0.0), Point::new(-999.9999995, 0.0), Point::new(-2000.0, 0.0), Point::new(0.0, -600.0)]; + let at_unit: Vec<_> = pts.iter().map(|&p| h().apply(p)).collect(); + assert!(at_unit[0].is_some() && at_unit[1].is_some()); + assert!(at_unit[2..].iter().all(Option::is_none)); + let inv = h().inverse().unwrap(); + // Two equal rows (up to a factor of 2), with entries that are exact in binary. + let flat = Homography([[1.0, 2.0, 3.0], [2.0, 4.0, 6.0], [0.25, 0.5, 1.0]]); + assert!(!flat.is_invertible()); + assert_eq!(flat.inverse(), None); + for f in [1e-200, 1e-12, 1e12, 1e200] { + let s = times(h(), f); + for (&p, want) in pts.iter().zip(&at_unit) { + match (s.apply(p), want) { + (Some(q), Some(w)) => assert!(q.distance(*w) <= 1e-6 * w.to_vec2().hypot(), "{f} {p:?}"), + (None, None) => {} + (got, _) => panic!("{f} {p:?}: {got:?}"), + } + } + assert!(s.is_invertible(), "{f}"); + let si = unit(s.inverse().unwrap_or_else(|| panic!("{f}"))); + for (a, b) in si.to_array().iter().zip(unit(inv).to_array()) { + assert!((a - b).abs() <= 1e-9 * b.abs().max(1.0), "{f}: {si:?}"); + } + let sf = times(flat, f); + assert!(!sf.is_invertible(), "{f}"); + assert_eq!(sf.inverse(), None, "{f}"); + } + } + #[test] fn a_product_that_overflows_is_none() { // Each factor is finite and invertible, but the product's first entry is 1e400. diff --git a/docs/target-app-parity.md b/docs/target-app-parity.md index 7e3241865..68ef6787b 100644 --- a/docs/target-app-parity.md +++ b/docs/target-app-parity.md @@ -1045,6 +1045,7 @@ Moved from the ROADMAP's "Shipped so far" on 2026-10-10. Update it when a featur Move in Perspective (its copy, perpendicular and Snap to Grid choices), Transform in Perspective (its copy), Nudge in Perspective (its ×10 and copy), Draw in Perspective (Snap to Grid), Define Perspective Grid (its reproject choice), Show Perspective Grid, Width Point (Adjust Adjoining Width Points), Liquify (its pen pressure, Simplify and Affect choices) and Puppet Warp (its pins' rest shape) refuse a value other than true, false or null, such as "true", and leave the document and its undo history as they were. Moves and transforms in perspective and Define Perspective Grid with reproject return an error, and leave the document and its undo history as they were, when the projection they would give type or a symbol instance in perspective has a number too large to store. Object › Perspective › Edit Text returns an error, and does not enter isolation mode, when the type cannot be shown flat where it is drawn: when part of it crosses the horizon, or when that view's projection has a number too large to store. Regressions: `a_product_that_overflows_is_none` in `crates/geom/src/projective.rs`, and `moving_type_whose_projection_would_overflow_is_refused` and `edit_text_whose_flat_view_would_overflow_is_refused` in `crates/engine/src/tests/persp_text.rs`. A stored projection with numbers so large that its determinant overflows is judged on a copy divided by a power of two: an invertible one is read back when the document opens as the same map at a smaller scale (divided by the largest number in its bottom row), so points that overflowed at the stored scale draw, and one that cannot be inverted, such as one whose bottom row is 0 at that scale, is not. Regression: `a_determinant_that_overflows_is_taken_at_a_smaller_scale` in `crates/geom/src/projective.rs`. + A projection gives the same answers when every one of its numbers is multiplied by the same factor: a point is past its horizon when its w is at most 1e-9 times the largest number in the bottom row, and the projection is invertible when its determinant is more than 1e-15 times the largest of the six products the determinant sums, so a stored projection with very small numbers draws and inverts, and one that cannot be inverted is refused at every scale. Regression: `the_checks_do_not_depend_on_the_scale_of_the_map` in `crates/geom/src/projective.rs`. Create Object Mosaic (its columns and rows), Transform Each (its reference point and Random seed) and Rearrange All Artboards (its columns) refuse a fraction, a negative number or a string, such as 2.5, -1 or "4"; 4.0 counts as 4. Move, Rotate, Scale, Reflect, Shear, Transform Again, Nudge and the Transform panel's values (their copy, Scale Strokes & Effects, Scale Corners, Transform Patterns, Transform Objects and area type options), Align, Distribute and Distribute Spacing (Transform Patterns), Object Properties (visible, locked, Isolate Blending and Opacity & Mask Define Knockout Shape) and Live Shape Properties (Invert Pie, Make Sides Equal and relative rounding) refuse a value other than true, false or null, such as "false", and leave the document as it was. Object Properties (its opacity), Live Shape Properties (corner radius and polygon sides), the Transform panel's reference point and Align, Distribute and Distribute Spacing to an artboard (the artboard's number) refuse a value other than a number or null, such as "50", and leave the document as it was. From bffb4e6e87ce8363d5c4755ca2c25e7dc8035afc Mon Sep 17 00:00:00 2001 From: dotCooCoo Date: Mon, 2 Nov 2026 07:21:32 -0800 Subject: [PATCH 2/8] Geometry: a projective map with entries of very different sizes stays invertible The scale-free determinant test no longer builds a rescaled copy of the map. It forms each of the six determinant products as a mantissa product and an exponent sum and compares their sum with the largest of them, so a map such as diag(1e307, 1e-157, 1e-157), whose rescaled copy overflows, is still invertible and loads from a document. The rescaled copy is built only to invert a map whose determinant is not finite or is below the normal range. --- crates/geom/src/projective.rs | 71 +++++++++++++++++++++++------------ 1 file changed, 46 insertions(+), 25 deletions(-) diff --git a/crates/geom/src/projective.rs b/crates/geom/src/projective.rs index a5ae446a5..73805baf1 100644 --- a/crates/geom/src/projective.rs +++ b/crates/geom/src/projective.rs @@ -72,32 +72,38 @@ impl Homography { + m[0][2] * (m[1][0] * m[2][1] - m[1][1] * m[2][0]) } - /// Whether the map is invertible: [`Self::scaled_det`] gives a determinant, or, when the - /// determinant is not finite, [`Self::inverse`] gives one. + /// Whether the map is invertible: [`Self::has_significant_det`], or, when the determinant is + /// not finite, [`Self::inverse`] gives one. fn is_invertible(&self) -> bool { if self.det().is_finite() { - return self.scaled_det().is_some(); + return self.has_significant_det(); } self.inverse().is_some() } - /// [`Self::scaled`] and its determinant; `None` when the magnitude of that determinant is at - /// most 1e-15 times that of the largest of the six products it sums (a test that gives the - /// same answer at any scale of the map). - fn scaled_det(&self) -> Option<(Self, f64)> { - let s = self.scaled()?; - let m = &s.0; - let terms = [ - m[0][0] * m[1][1] * m[2][2], - m[0][0] * m[1][2] * m[2][1], - m[0][1] * m[1][0] * m[2][2], - m[0][1] * m[1][2] * m[2][0], - m[0][2] * m[1][0] * m[2][1], - m[0][2] * m[1][1] * m[2][0], - ]; + /// Whether every entry is finite and the determinant's magnitude is more than 1e-15 times that + /// of the largest of the six products it sums (a test that gives the same answer at any scale + /// of the map). Each product is formed as a mantissa product and an exponent sum, and all six + /// are divided by the power of two of the largest before they are summed, so the largest is at + /// least 1/8 and below 1 and none of them overflows. + fn has_significant_det(&self) -> bool { + if !self.is_finite() { + return false; + } + let f = self.0.map(|row| row.map(frexp)); + // The columns each of the six products takes from rows 0, 1 and 2, and its sign. + let cols = [(0, 1, 2, 1.0), (0, 2, 1, -1.0), (1, 0, 2, -1.0), (1, 2, 0, 1.0), (2, 0, 1, 1.0), (2, 1, 0, -1.0)]; + let terms = cols.map(|(a, b, c, sign)| { + let (a, b, c) = (f[0][a], f[1][b], f[2][c]); + (sign * a.0 * b.0 * c.0, a.1 + b.1 + c.1) + }); + let Some(e) = terms.iter().filter(|(m, _)| *m != 0.0).map(|(_, e)| *e).max() else { + return false; + }; + let terms = terms.map(|(m, x)| ldexp(m, x - e)); let largest = terms.iter().fold(0.0_f64, |k, v| k.max(v.abs())); - let d = s.det(); - (d.abs() > 1e-15 * largest).then_some((s, d)) + let d: f64 = terms.iter().sum(); + d.abs() > 1e-15 * largest } /// The same map with every entry divided by one power of two, chosen so that the largest of the @@ -120,14 +126,17 @@ impl Homography { Some(Self(self.0.map(|row| row.map(|v| ldexp(v, -e))))) } - /// Inverse; `None` when [`Self::scaled_det`] counts the map as not invertible or its inverse - /// has an entry that is not finite. When the determinant is not finite or is below the normal - /// range, this is the inverse of [`Self::scaled`] (the exact inverse multiplied by one power of - /// two, which is the same map). + /// Inverse; `None` when [`Self::has_significant_det`] counts the map as not invertible or its + /// inverse has an entry that is not finite. When the determinant is not finite or is below the + /// normal range, this is the inverse of [`Self::scaled`] (the exact inverse multiplied by one + /// power of two, which is the same map). pub fn inverse(&self) -> Option { - let (s, sd) = self.scaled_det()?; + if !self.has_significant_det() { + return None; + } let d = self.det(); - let (h, d) = if d.is_finite() && d.abs() >= f64::MIN_POSITIVE { (*self, d) } else { (s, sd) }; + let h = if d.is_finite() && d.abs() >= f64::MIN_POSITIVE { *self } else { self.scaled()? }; + let d = h.det(); let m = &h.0; let c = |r0: usize, c0: usize, r1: usize, c1: usize| m[r0][c0] * m[r1][c1] - m[r0][c1] * m[r1][c0]; let adj = [ @@ -443,6 +452,18 @@ mod tests { assert!(!sf.is_invertible(), "{f}"); assert_eq!(sf.inverse(), None, "{f}"); } + // Entries of very different sizes: the determinant is about 1e-7 and every entry of the + // inverse is finite, though one power of two that brings the products near 1 takes the + // first entry past the largest finite number. + let uneven = Homography([[1e307, 0.0, 0.0], [0.0, 1e-157, 0.0], [0.0, 0.0, 1e-157]]); + assert!(uneven.det().is_finite()); + assert_eq!(Homography::from_array(uneven.to_array()), Some(uneven)); + let inv = uneven.inverse().unwrap(); + for p in [Point::new(1e-300, 2.0), Point::new(-3e-299, -0.5)] { + let q = uneven.apply(p).unwrap(); + let back = inv.apply(q).unwrap(); + assert!((back.x - p.x).abs() <= 1e-9 * p.x.abs() && (back.y - p.y).abs() <= 1e-9 * p.y.abs(), "{p:?} {back:?}"); + } } #[test] From b22912b25a237de9f915794a5d297be4b8e14927 Mon Sep 17 00:00:00 2001 From: dotCooCoo Date: Mon, 2 Nov 2026 07:37:37 -0800 Subject: [PATCH 3/8] Geometry: the inverse of a projective map with a tiny determinant keeps its w row When a cofactor of a map with a normal but very small determinant underflows to 0, the inverse formed from the map itself can have a w row of 0, so every point mapped through it is refused. The inverse is now kept only when it is itself invertible, and is otherwise formed from the map divided by a power of two that brings its determinant's products near 1. --- crates/geom/src/projective.rs | 43 ++++++++++++++++++++++++++--------- 1 file changed, 32 insertions(+), 11 deletions(-) diff --git a/crates/geom/src/projective.rs b/crates/geom/src/projective.rs index 73805baf1..6fe0cfdce 100644 --- a/crates/geom/src/projective.rs +++ b/crates/geom/src/projective.rs @@ -126,18 +126,31 @@ impl Homography { Some(Self(self.0.map(|row| row.map(|v| ldexp(v, -e))))) } - /// Inverse; `None` when [`Self::has_significant_det`] counts the map as not invertible or its - /// inverse has an entry that is not finite. When the determinant is not finite or is below the - /// normal range, this is the inverse of [`Self::scaled`] (the exact inverse multiplied by one - /// power of two, which is the same map). + /// Inverse; `None` when [`Self::has_significant_det`] counts the map as not invertible. The + /// inverse is formed from the map itself when its determinant is finite and in the normal + /// range, and is kept when every entry is finite and it passes [`Self::has_significant_det`]. + /// Otherwise it is the inverse of [`Self::scaled`] (the exact inverse multiplied by one power + /// of two, which is the same map), and `None` when that one fails the same checks. pub fn inverse(&self) -> Option { if !self.has_significant_det() { return None; } let d = self.det(); - let h = if d.is_finite() && d.abs() >= f64::MIN_POSITIVE { *self } else { self.scaled()? }; - let d = h.det(); - let m = &h.0; + if d.is_finite() + && d.abs() >= f64::MIN_POSITIVE + && let Some(inv) = self.adjugate_inverse() + { + return Some(inv); + } + self.scaled()?.adjugate_inverse() + } + + /// The adjugate divided by the determinant; `None` when an entry is not finite or the result + /// fails [`Self::has_significant_det`] (a cofactor that underflows to 0 can leave a whole row + /// at 0). + fn adjugate_inverse(&self) -> Option { + let d = self.det(); + let m = &self.0; let c = |r0: usize, c0: usize, r1: usize, c1: usize| m[r0][c0] * m[r1][c1] - m[r0][c1] * m[r1][c0]; let adj = [ [c(1, 1, 2, 2), -c(0, 1, 2, 2), c(0, 1, 1, 2)], @@ -150,11 +163,9 @@ impl Homography { out[r][k] = v / d; } } - if out.iter().flatten().any(|v| !v.is_finite()) { - return None; - } + let out = Self(out); // H⁻¹·(x, y, 1) = (u, v, 1)/w_H, so points in front of the camera keep w > 0. - Some(Self(out)) + out.has_significant_det().then_some(out) } /// `self ∘ other`: `other` first. @@ -464,6 +475,16 @@ mod tests { let back = inv.apply(q).unwrap(); assert!((back.x - p.x).abs() <= 1e-9 * p.x.abs() && (back.y - p.y).abs() <= 1e-9 * p.y.abs(), "{p:?} {back:?}"); } + // A normal determinant (1e-200) whose cofactor 1e-200 · 1e-200 is below the smallest + // positive f64: the inverse still has a w row that is not 0. + let tiny = Homography([[1e-200, 0.0, 0.0], [0.0, 1e-200, 0.0], [0.0, 0.0, 1e200]]); + assert_eq!(Homography::from_array(tiny.to_array()), Some(tiny)); + let inv = tiny.inverse().unwrap(); + for p in [Point::new(1e300, -2e299), Point::new(-4e299, 5e300)] { + let q = tiny.apply(p).unwrap(); + let back = inv.apply(q).unwrap(); + assert!((back.x - p.x).abs() <= 1e-9 * p.x.abs() && (back.y - p.y).abs() <= 1e-9 * p.y.abs(), "{p:?} {back:?}"); + } } #[test] From d1fb8e540cfc176685eaa025db24d2041d866fcc Mon Sep 17 00:00:00 2001 From: dotCooCoo Date: Mon, 2 Nov 2026 07:45:10 -0800 Subject: [PATCH 4/8] Geometry: a projective map with entries far apart in size keeps its inverse The inverse holds every cofactor and the determinant as a mantissa and a power of two, and forms each entry from their quotient. A map such as diag(1e-300, 1e300, 1e-300), which counts as invertible, now gets its inverse diag(1e300, 1e-300, 1e300) instead of none. --- crates/geom/src/projective.rs | 68 +++++++++++++++++++++++++---------- 1 file changed, 50 insertions(+), 18 deletions(-) diff --git a/crates/geom/src/projective.rs b/crates/geom/src/projective.rs index 6fe0cfdce..baf057065 100644 --- a/crates/geom/src/projective.rs +++ b/crates/geom/src/projective.rs @@ -127,17 +127,15 @@ impl Homography { } /// Inverse; `None` when [`Self::has_significant_det`] counts the map as not invertible. The - /// inverse is formed from the map itself when its determinant is finite and in the normal - /// range, and is kept when every entry is finite and it passes [`Self::has_significant_det`]. - /// Otherwise it is the inverse of [`Self::scaled`] (the exact inverse multiplied by one power - /// of two, which is the same map), and `None` when that one fails the same checks. + /// inverse is formed from the map itself when its determinant is finite, and is kept when every + /// entry is finite and it passes [`Self::has_significant_det`]. Otherwise it is the inverse of + /// [`Self::scaled`] (the exact inverse multiplied by one power of two, which is the same map), + /// and `None` when that one fails the same checks. pub fn inverse(&self) -> Option { if !self.has_significant_det() { return None; } - let d = self.det(); - if d.is_finite() - && d.abs() >= f64::MIN_POSITIVE + if self.det().is_finite() && let Some(inv) = self.adjugate_inverse() { return Some(inv); @@ -145,22 +143,35 @@ impl Homography { self.scaled()?.adjugate_inverse() } - /// The adjugate divided by the determinant; `None` when an entry is not finite or the result - /// fails [`Self::has_significant_det`] (a cofactor that underflows to 0 can leave a whole row - /// at 0). + /// The adjugate divided by the determinant; `None` when an entry is not finite, the + /// determinant is 0, or the result fails [`Self::has_significant_det`]. Every cofactor and the + /// determinant are held as a mantissa and a power of two, so an entry of the result is lost + /// only when that entry itself is beyond the range of `f64`. fn adjugate_inverse(&self) -> Option { - let d = self.det(); - let m = &self.0; - let c = |r0: usize, c0: usize, r1: usize, c1: usize| m[r0][c0] * m[r1][c1] - m[r0][c1] * m[r1][c0]; + if !self.is_finite() { + return None; + } + let f = self.0.map(|row| row.map(frexp)); + let mul = |a: (f64, i32), b: (f64, i32)| (a.0 * b.0, a.1 + b.1); + let c = |r0: usize, c0: usize, r1: usize, c1: usize| { + let (a, b) = (mul(f[r0][c0], f[r1][c1]), mul(f[r0][c1], f[r1][c0])); + split_sum(&[a, (-b.0, b.1)]) + }; + let neg = |(m, e): (f64, i32)| (-m, e); let adj = [ - [c(1, 1, 2, 2), -c(0, 1, 2, 2), c(0, 1, 1, 2)], - [-c(1, 0, 2, 2), c(0, 0, 2, 2), -c(0, 0, 1, 2)], - [c(1, 0, 2, 1), -c(0, 0, 2, 1), c(0, 0, 1, 1)], + [c(1, 1, 2, 2), neg(c(0, 1, 2, 2)), c(0, 1, 1, 2)], + [neg(c(1, 0, 2, 2)), c(0, 0, 2, 2), neg(c(0, 0, 1, 2))], + [c(1, 0, 2, 1), neg(c(0, 0, 2, 1)), c(0, 0, 1, 1)], ]; + // The determinant expanded along the first row, as `det` sums it. + let (dm, de) = split_sum(&[mul(f[0][0], adj[0][0]), mul(f[0][1], adj[1][0]), mul(f[0][2], adj[2][0])]); + if dm == 0.0 { + return None; + } let mut out = [[0.0; 3]; 3]; for (r, row) in adj.iter().enumerate() { - for (k, v) in row.iter().enumerate() { - out[r][k] = v / d; + for (k, (m, e)) in row.iter().enumerate() { + out[r][k] = ldexp(m / dm, e - de); } } let out = Self(out); @@ -269,6 +280,17 @@ fn frexp(v: f64) -> (f64, i32) { (f64::from_bits((bits & !(0x7ff << 52)) | (1022 << 52)), e + bias) } +/// The sum of `m·2^e` over `terms`, in order, as a mantissa `0.5 ≤ |m| < 1` (or 0) and a power of +/// two. Every term is first divided by the power of two of the largest nonzero one, so none of +/// them overflows. +fn split_sum(terms: &[(f64, i32)]) -> (f64, i32) { + let Some(e) = terms.iter().filter(|(m, _)| *m != 0.0).map(|(_, e)| *e).max() else { + return (0.0, 0); + }; + let (m, x) = frexp(terms.iter().map(|&(m, x)| ldexp(m, x - e)).sum()); + (m, e + x) +} + /// `m·2^e`, multiplied in steps of powers of two that are each within the normal range. fn ldexp(mut m: f64, mut e: i32) -> f64 { while e > 1023 && m.is_finite() { @@ -485,6 +507,16 @@ mod tests { let back = inv.apply(q).unwrap(); assert!((back.x - p.x).abs() <= 1e-9 * p.x.abs() && (back.y - p.y).abs() <= 1e-9 * p.y.abs(), "{p:?} {back:?}"); } + // A normal determinant (1e-300) whose cofactor 1e-300 · 1e-300 is below the smallest + // positive f64, and where one power of two that brings the products near 1 takes 1e300 past + // the largest finite number: the inverse is still diag(1e300, 1e-300, 1e300). + let spread = Homography([[1e-300, 0.0, 0.0], [0.0, 1e300, 0.0], [0.0, 0.0, 1e-300]]); + assert_eq!(Homography::from_array(spread.to_array()), Some(spread)); + let inv = spread.inverse().unwrap(); + let want = [1e300, 0.0, 0.0, 0.0, 1e-300, 0.0, 0.0, 0.0, 1e300]; + for (v, w) in inv.to_array().iter().zip(want) { + assert!((v - w).abs() <= 1e-12 * w.abs(), "{inv:?}"); + } } #[test] From 4494d40c333474ad3b7ff33ccd697405e236561e Mon Sep 17 00:00:00 2001 From: dotCooCoo Date: Mon, 2 Nov 2026 07:55:19 -0800 Subject: [PATCH 5/8] Geometry: the inverse of an invertible projective map is kept however small its relative determinant adjugate_inverse no longer applies the relative determinant test to its result. That test is not symmetric between a map and its inverse: a map whose determinant is 1e-8 of its largest product can have a finite, accurate inverse whose own ratio is near 1e-16, so inverse() was None while is_invertible() was true. The result is now kept when every entry is finite and its determinant, summed from the six products in split form, is not 0. --- crates/geom/src/projective.rs | 59 ++++++++++++++++++++++++++++------- 1 file changed, 47 insertions(+), 12 deletions(-) diff --git a/crates/geom/src/projective.rs b/crates/geom/src/projective.rs index baf057065..2a9ff9093 100644 --- a/crates/geom/src/projective.rs +++ b/crates/geom/src/projective.rs @@ -90,13 +90,7 @@ impl Homography { if !self.is_finite() { return false; } - let f = self.0.map(|row| row.map(frexp)); - // The columns each of the six products takes from rows 0, 1 and 2, and its sign. - let cols = [(0, 1, 2, 1.0), (0, 2, 1, -1.0), (1, 0, 2, -1.0), (1, 2, 0, 1.0), (2, 0, 1, 1.0), (2, 1, 0, -1.0)]; - let terms = cols.map(|(a, b, c, sign)| { - let (a, b, c) = (f[0][a], f[1][b], f[2][c]); - (sign * a.0 * b.0 * c.0, a.1 + b.1 + c.1) - }); + let terms = self.det_terms(); let Some(e) = terms.iter().filter(|(m, _)| *m != 0.0).map(|(_, e)| *e).max() else { return false; }; @@ -106,6 +100,18 @@ impl Homography { d.abs() > 1e-15 * largest } + /// The six signed products the determinant sums, each as a mantissa product and an exponent + /// sum. + fn det_terms(&self) -> [(f64, i32); 6] { + let f = self.0.map(|row| row.map(frexp)); + // The columns each of the six products takes from rows 0, 1 and 2, and its sign. + let cols = [(0, 1, 2, 1.0), (0, 2, 1, -1.0), (1, 0, 2, -1.0), (1, 2, 0, 1.0), (2, 0, 1, 1.0), (2, 1, 0, -1.0)]; + cols.map(|(a, b, c, sign)| { + let (a, b, c) = (f[0][a], f[1][b], f[2][c]); + (sign * a.0 * b.0 * c.0, a.1 + b.1 + c.1) + }) + } + /// The same map with every entry divided by one power of two, chosen so that the largest of the /// six products its determinant sums is at least 1/8 and below 4. `None` when every one of /// those products is 0 or an entry is not finite. @@ -128,7 +134,7 @@ impl Homography { /// Inverse; `None` when [`Self::has_significant_det`] counts the map as not invertible. The /// inverse is formed from the map itself when its determinant is finite, and is kept when every - /// entry is finite and it passes [`Self::has_significant_det`]. Otherwise it is the inverse of + /// entry is finite and its own determinant is not 0. Otherwise it is the inverse of /// [`Self::scaled`] (the exact inverse multiplied by one power of two, which is the same map), /// and `None` when that one fails the same checks. pub fn inverse(&self) -> Option { @@ -144,9 +150,10 @@ impl Homography { } /// The adjugate divided by the determinant; `None` when an entry is not finite, the - /// determinant is 0, or the result fails [`Self::has_significant_det`]. Every cofactor and the - /// determinant are held as a mantissa and a power of two, so an entry of the result is lost - /// only when that entry itself is beyond the range of `f64`. + /// determinant is 0, or the result has an entry that is not finite or a determinant of 0 (its + /// six products, each held as a mantissa and a power of two, sum to 0, as when a row underflows + /// to 0). Every cofactor and the determinant are held as a mantissa and a power of two, so an + /// entry of the result is lost only when that entry itself is beyond the range of `f64`. fn adjugate_inverse(&self) -> Option { if !self.is_finite() { return None; @@ -176,7 +183,7 @@ impl Homography { } let out = Self(out); // H⁻¹·(x, y, 1) = (u, v, 1)/w_H, so points in front of the camera keep w > 0. - out.has_significant_det().then_some(out) + (out.is_finite() && split_sum(&out.det_terms()).0 != 0.0).then_some(out) } /// `self ∘ other`: `other` first. @@ -519,6 +526,34 @@ mod tests { } } + #[test] + fn an_invertible_map_has_an_inverse_at_any_scale() { + let times = |h: Homography, f: f64| Homography(h.0.map(|row| row.map(|v| v * f))); + // Determinant 1e-17 against a largest product of 1e-9 (and 1e-8 at 1e3 times the map): its + // inverse is finite, though that inverse's own determinant is near 1e-16 of its largest + // product. + let skew = Homography([[-1000.0, -1e-8, -1e-5], [0.0, 1e-4, -1e7], [1e-8, 0.0, 1e-8]]); + let exact = [1e5, 10.0, 1e16 + 1e8, -1e16, -1e12 + 1e4, -1e27, -1e5, -10.0, -1e16]; + for f in [1.0, 1e3] { + let s = times(skew, f); + assert_eq!(Homography::from_array(s.to_array()), Some(s), "{f}"); + let inv = s.inverse().unwrap_or_else(|| panic!("{f}")); + for (v, w) in inv.to_array().iter().zip(exact) { + assert!((v * f - w).abs() <= 1e-9 * w.abs(), "{f}: {inv:?}"); + } + // The map times its inverse is the identity, to within the rounding of the products + // each entry sums (they cancel by up to 1e23 in the second row). + let id = s.then_after(&inv); + for (r, row) in id.0.iter().enumerate() { + for (k, v) in row.iter().enumerate() { + let size: f64 = (0..3).map(|i| (s.0[r][i] * inv.0[i][k]).abs()).sum(); + let want = if r == k { 1.0 } else { 0.0 }; + assert!((v - want).abs() <= 1e-12 * size, "{f}: {id:?}"); + } + } + } + } + #[test] fn a_product_that_overflows_is_none() { // Each factor is finite and invertible, but the product's first entry is 1e400. From 0f48ba21e9d8c260e756c61fc7ceec421b903c1b Mon Sep 17 00:00:00 2001 From: dotCooCoo Date: Mon, 2 Nov 2026 08:15:21 -0800 Subject: [PATCH 6/8] Geometry: a projective map at a tiny scale maps points and derivatives as at unit scale apply and affine_at multiply every entry by the power of two that brings the largest to at least 1/2 when all entries are below 1/2, so products that would be subnormal keep their precision. The identity times 5e-324 mapped (0.5, 0.5) to (0, 0) before. Maps with an entry of 1/2 or more are unchanged. --- crates/geom/src/projective.rs | 36 ++++++++++++++++++++++++++++++++--- 1 file changed, 33 insertions(+), 3 deletions(-) diff --git a/crates/geom/src/projective.rs b/crates/geom/src/projective.rs index 2a9ff9093..9d4800626 100644 --- a/crates/geom/src/projective.rs +++ b/crates/geom/src/projective.rs @@ -52,9 +52,9 @@ impl Homography { /// Map a point; `None` when it lands on or beyond the horizon (w at most 1e-9 times the largest /// magnitude in the map's `w` row, so the answer is the same at any scale of the map), or when - /// it does not map to a finite page point. + /// it does not map to a finite page point. The point is mapped through [`Self::upscaled`]. pub fn apply(&self, p: Point) -> Option { - let m = &self.0; + let m = self.upscaled(); let x = m[0][0] * p.x + m[0][1] * p.y + m[0][2]; let y = m[1][0] * p.x + m[1][1] * p.y + m[1][2]; let w = m[2][0] * p.x + m[2][1] * p.y + m[2][2]; @@ -66,6 +66,15 @@ impl Homography { q.is_finite().then_some(q) } + /// The entries, multiplied by the power of two that brings the largest magnitude to at least + /// 1/2 and below 1 when every entry is below 1/2 in magnitude (an exact step that keeps a map + /// at a tiny scale from losing precision to subnormal products), and unchanged otherwise. + fn upscaled(&self) -> [[f64; 3]; 3] { + let largest = self.0.iter().flatten().fold(0.0_f64, |k, v| k.max(v.abs())); + let e = frexp(largest).1.min(0); + self.0.map(|row| row.map(|v| ldexp(v, -e))) + } + pub fn det(&self) -> f64 { let m = &self.0; m[0][0] * (m[1][1] * m[2][2] - m[1][2] * m[2][1]) - m[0][1] * (m[1][0] * m[2][2] - m[1][2] * m[2][0]) @@ -252,7 +261,7 @@ impl Homography { /// The affine map that matches `self` to first order at `p` (its derivative there), `None` /// beyond the horizon or when the derivative is not finite. pub fn affine_at(&self, p: Point) -> Option { - let m = &self.0; + let m = self.upscaled(); let w = m[2][0] * p.x + m[2][1] * p.y + m[2][2]; let q = self.apply(p)?; // d(x/w) = (dx − (x/w)·dw) / w, and the same for y. @@ -554,6 +563,27 @@ mod tests { } } + #[test] + fn a_map_at_a_tiny_scale_maps_points_as_at_unit_scale() { + for f in [5e-324, 1e-300] { + let h = Homography(Homography::IDENTITY.0.map(|row| row.map(|v| v * f))); + assert_eq!(Homography::from_array(h.to_array()), Some(h), "{f}"); + for p in [Point::new(0.5, 0.5), Point::new(1e-20, 1e-20)] { + let q = h.apply(p).unwrap_or_else(|| panic!("{f}: {p:?}")); + assert!((q.x - p.x).abs() <= 1e-12 * p.x && (q.y - p.y).abs() <= 1e-12 * p.y, "{f}: {p:?} -> {q:?}"); + } + } + // A perspective map whose entries are exact at 2^-1060 times their size, so the scaled + // copy holds the same map. + let unit = Homography([[1.0, 0.0, 2.0], [0.0, 1.0, 3.0], [0.5, 0.25, 1.0]]); + let tiny = Homography(unit.0.map(|row| row.map(|v| v * pow2(-1022) * pow2(-38)))); + let p = Point::new(0.3, 0.7); + let (a, b) = (tiny.affine_at(p).unwrap(), unit.affine_at(p).unwrap()); + for (u, v) in a.as_coeffs().iter().zip(b.as_coeffs()) { + assert!((u - v).abs() <= 1e-12 * v.abs().max(1.0), "{a:?} {b:?}"); + } + } + #[test] fn a_product_that_overflows_is_none() { // Each factor is finite and invertible, but the product's first entry is 1e400. From b5885e35fa7461b5036ff7ac9fc64dfc05dd88c1 Mon Sep 17 00:00:00 2001 From: dotCooCoo Date: Mon, 2 Nov 2026 08:24:29 -0800 Subject: [PATCH 7/8] Geometry: a projective map at a small scale maps a far point without overflowing The power of two that brings a small map's entries up before a point is mapped is now reduced by the power of two of the point's larger coordinate, so no product of an entry and a coordinate reaches 1. A map that mapped a point before the scaling step does so again, and affine_at uses the same entries. --- crates/geom/src/projective.rs | 28 +++++++++++++++++++++------- 1 file changed, 21 insertions(+), 7 deletions(-) diff --git a/crates/geom/src/projective.rs b/crates/geom/src/projective.rs index 9d4800626..070048f38 100644 --- a/crates/geom/src/projective.rs +++ b/crates/geom/src/projective.rs @@ -54,7 +54,7 @@ impl Homography { /// magnitude in the map's `w` row, so the answer is the same at any scale of the map), or when /// it does not map to a finite page point. The point is mapped through [`Self::upscaled`]. pub fn apply(&self, p: Point) -> Option { - let m = self.upscaled(); + let m = self.upscaled(p); let x = m[0][0] * p.x + m[0][1] * p.y + m[0][2]; let y = m[1][0] * p.x + m[1][1] * p.y + m[1][2]; let w = m[2][0] * p.x + m[2][1] * p.y + m[2][2]; @@ -66,12 +66,15 @@ impl Homography { q.is_finite().then_some(q) } - /// The entries, multiplied by the power of two that brings the largest magnitude to at least - /// 1/2 and below 1 when every entry is below 1/2 in magnitude (an exact step that keeps a map - /// at a tiny scale from losing precision to subnormal products), and unchanged otherwise. - fn upscaled(&self) -> [[f64; 3]; 3] { + /// The entries for mapping `p`, multiplied by a power of two of at least 1 (an exact step that + /// keeps a map at a tiny scale from losing precision to subnormal products). When every entry + /// is below 1/2 in magnitude, the power brings the largest to at least 1/2 and below 1, divided + /// by the power of two of `p`'s larger coordinate magnitude when that is 1 or more, so that no + /// product of an entry and a coordinate reaches 1. Otherwise the entries are unchanged. + fn upscaled(&self, p: Point) -> [[f64; 3]; 3] { let largest = self.0.iter().flatten().fold(0.0_f64, |k, v| k.max(v.abs())); - let e = frexp(largest).1.min(0); + let e_p = frexp(p.x.abs().max(p.y.abs())).1.max(0); + let e = (frexp(largest).1.min(0) + e_p).min(0); self.0.map(|row| row.map(|v| ldexp(v, -e))) } @@ -261,7 +264,7 @@ impl Homography { /// The affine map that matches `self` to first order at `p` (its derivative there), `None` /// beyond the horizon or when the derivative is not finite. pub fn affine_at(&self, p: Point) -> Option { - let m = self.upscaled(); + let m = self.upscaled(p); let w = m[2][0] * p.x + m[2][1] * p.y + m[2][2]; let q = self.apply(p)?; // d(x/w) = (dx − (x/w)·dw) / w, and the same for y. @@ -584,6 +587,17 @@ mod tests { } } + #[test] + fn a_map_at_a_small_scale_maps_a_far_point() { + // Every product of an entry and a coordinate is finite, but the entries brought to 1/2 or + // more would overflow against this point. + let h = Homography([[1e-6, 1e-6, 0.0], [0.0, 1e-6, 0.0], [1e-6, 0.0, 1e-6]]); + let p = Point::new(1.75e308, 1.75e308); + let q = h.apply(p).unwrap_or_else(|| panic!("{p:?}")); + assert!((q.x - 2.0).abs() <= 1e-12 && (q.y - 1.0).abs() <= 1e-12, "{q:?}"); + assert!(h.affine_at(p).is_some()); + } + #[test] fn a_product_that_overflows_is_none() { // Each factor is finite and invertible, but the product's first entry is 1e400. From 2d0fa32d303ee77d8ed2c55709c5bf2efa44d054 Mon Sep 17 00:00:00 2001 From: dotCooCoo Date: Mon, 2 Nov 2026 08:53:47 -0800 Subject: [PATCH 8/8] Docs: the projection checks compare absolute values The horizon test compares a point's signed w with the largest absolute value in the bottom row, and the invertibility test compares the absolute value of the determinant with the largest absolute value among its six products, so a negative determinant counts as invertible. A stored projection that overflows is divided by the largest absolute value in its bottom row. --- docs/target-app-parity.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/target-app-parity.md b/docs/target-app-parity.md index 68ef6787b..043058ec0 100644 --- a/docs/target-app-parity.md +++ b/docs/target-app-parity.md @@ -1044,8 +1044,8 @@ Moved from the ROADMAP's "Shipped so far" on 2026-10-10. Update it when a featur All Artwork Above and Other Layers (Object › Lock and Object › Hide), Transform Each, Rearrange All Artboards, the width point and width profile commands, Liquify, Puppet Warp, Define Perspective Grid, Attach to Active Plane, Release with Perspective, moves and transforms in perspective, Move Plane and Move Plane to Match Object record no undo step when they change nothing. Move in Perspective (its copy, perpendicular and Snap to Grid choices), Transform in Perspective (its copy), Nudge in Perspective (its ×10 and copy), Draw in Perspective (Snap to Grid), Define Perspective Grid (its reproject choice), Show Perspective Grid, Width Point (Adjust Adjoining Width Points), Liquify (its pen pressure, Simplify and Affect choices) and Puppet Warp (its pins' rest shape) refuse a value other than true, false or null, such as "true", and leave the document and its undo history as they were. Moves and transforms in perspective and Define Perspective Grid with reproject return an error, and leave the document and its undo history as they were, when the projection they would give type or a symbol instance in perspective has a number too large to store. Object › Perspective › Edit Text returns an error, and does not enter isolation mode, when the type cannot be shown flat where it is drawn: when part of it crosses the horizon, or when that view's projection has a number too large to store. Regressions: `a_product_that_overflows_is_none` in `crates/geom/src/projective.rs`, and `moving_type_whose_projection_would_overflow_is_refused` and `edit_text_whose_flat_view_would_overflow_is_refused` in `crates/engine/src/tests/persp_text.rs`. - A stored projection with numbers so large that its determinant overflows is judged on a copy divided by a power of two: an invertible one is read back when the document opens as the same map at a smaller scale (divided by the largest number in its bottom row), so points that overflowed at the stored scale draw, and one that cannot be inverted, such as one whose bottom row is 0 at that scale, is not. Regression: `a_determinant_that_overflows_is_taken_at_a_smaller_scale` in `crates/geom/src/projective.rs`. - A projection gives the same answers when every one of its numbers is multiplied by the same factor: a point is past its horizon when its w is at most 1e-9 times the largest number in the bottom row, and the projection is invertible when its determinant is more than 1e-15 times the largest of the six products the determinant sums, so a stored projection with very small numbers draws and inverts, and one that cannot be inverted is refused at every scale. Regression: `the_checks_do_not_depend_on_the_scale_of_the_map` in `crates/geom/src/projective.rs`. + A stored projection with numbers so large that its determinant overflows is judged on a copy divided by a power of two: an invertible one is read back when the document opens as the same map at a smaller scale (divided by the largest absolute value in its bottom row), so points that overflowed at the stored scale draw, and one that cannot be inverted, such as one whose bottom row is 0 at that scale, is not. Regression: `a_determinant_that_overflows_is_taken_at_a_smaller_scale` in `crates/geom/src/projective.rs`. + A projection gives the same answers when every one of its numbers is multiplied by the same factor. A point is on or past its horizon when its w is at most 1e-9 times the largest absolute value in the bottom row, so a point with a negative w is past it. The projection is invertible when the absolute value of its determinant is more than 1e-15 times the largest absolute value among the six products the determinant sums, so a negative determinant counts the same as a positive one. A stored projection with very small numbers draws and inverts, and one that cannot be inverted is refused at every scale. Regression: `the_checks_do_not_depend_on_the_scale_of_the_map` in `crates/geom/src/projective.rs`. Create Object Mosaic (its columns and rows), Transform Each (its reference point and Random seed) and Rearrange All Artboards (its columns) refuse a fraction, a negative number or a string, such as 2.5, -1 or "4"; 4.0 counts as 4. Move, Rotate, Scale, Reflect, Shear, Transform Again, Nudge and the Transform panel's values (their copy, Scale Strokes & Effects, Scale Corners, Transform Patterns, Transform Objects and area type options), Align, Distribute and Distribute Spacing (Transform Patterns), Object Properties (visible, locked, Isolate Blending and Opacity & Mask Define Knockout Shape) and Live Shape Properties (Invert Pie, Make Sides Equal and relative rounding) refuse a value other than true, false or null, such as "false", and leave the document as it was. Object Properties (its opacity), Live Shape Properties (corner radius and polygon sides), the Transform panel's reference point and Align, Distribute and Distribute Spacing to an artboard (the artboard's number) refuse a value other than a number or null, such as "50", and leave the document as it was.