src/fe_model/p2_local_model.rs

changeset 208
6be69f736c79
parent 206
ce37ff3ce507
child 209
060891c3f537
--- a/src/fe_model/p2_local_model.rs	Thu Sep 03 10:50:27 2026 -0500
+++ b/src/fe_model/p2_local_model.rs	Thu Sep 03 10:59:13 2026 -0500
@@ -8,7 +8,7 @@
 use crate::linsolve::*;
 use crate::loc::Loc;
 use crate::sets::Cube;
-use crate::sets::{NPolygon, Set, SpannedHalfspace};
+use crate::sets::Set;
 use crate::types::*;
 use numeric_literals::replace_float_literals;
 
@@ -102,13 +102,13 @@
     #[inline]
     fn p2powers(&self) -> Self::Output {
         let &Loc([x0, x1]) = self;
-        [x0 * x0, x0 * x1, x1 * x1].into()
+        [x0 * x0, 2.0 * x0 * x1, x1 * x1].into()
     }
 
     #[inline]
     fn p2powers_full(&self) -> Self::Full {
         let &Loc([x0, x1]) = self;
-        [1.0, x0, x1, x0 * x0, x0 * x1, x1 * x1].into()
+        [1.0, x0, x1, x0 * x0, 2.0 * x0 * x1, x1 * x1].into()
     }
 
     #[inline]
@@ -328,13 +328,14 @@
         let &P2LocalModel { a0, a1: Loc([a1, a2]), a2: Loc([a11, a12, a22]), .. } = self;
         let &Loc([x00, x01]) = x0;
         let d @ Loc([d0, d1]) = x1 - x0;
-        let b0 = a0 + a1 * x00 + a2 * x01 + a11 * x00 * x00 + a12 * x00 * x01 + a22 * x01 * x01;
+        let b0 =
+            a0 + a1 * x00 + a2 * x01 + a11 * x00 * x00 + 2.0 * a12 * x00 * x01 + a22 * x01 * x01;
         let b1 = a1 * d0
             + a2 * d1
             + 2.0 * a11 * d0 * x00
-            + a12 * (d0 * x01 + d1 * x00)
+            + 2.0 * a12 * (d0 * x01 + d1 * x00)
             + 2.0 * a22 * d1 * x01;
-        let b11 = a11 * d0 * d0 + a12 * d0 * d1 + a22 * d1 * d1;
+        let b11 = a11 * d0 * d0 + 2.0 * a12 * d0 * d1 + a22 * d1 * d1;
         let edge_1d_model = P2LocalModel {
             a0: b0,
             a1: Loc([b1]),
@@ -364,9 +365,9 @@
         let r = 2.0 * (a11 * a22 - a12 * a12);
         if r > 0.0 {
             // An interior solution (x[1], x[2]) has to satisfy
-            //   2a₁₁*x[1] + 2a₁₂*x[2]+a₁ =0 and 2a₂₂*x[1] + 2a₁₂*x[1]+a₂=0
+            //   2a₁₁*x[1] + 2a₁₂*x[2]+a₁ =0 and 2a₂₂*x[2] + 2a₁₂*x[1]+a₂=0
             // This gives
-            let x = [(a22 * a1 - a12 * a2) / r, (a12 * a1 - a11 * a2) / r].into();
+            let x = [(a22 * a1 - a12 * a2) / r, (a11 * a2 - a12 * a1) / r].into();
             if el.contains(&x) {
                 return (x, self.value(&x));
             }

mercurial