From 7366a6c53deb8ddf94694fe88a1edab552b0d64f Mon Sep 17 00:00:00 2001
From: phkahler <14852918+phkahler@users.noreply.github.com>
Date: Sat, 2 May 2020 14:56:50 -0400
Subject: [PATCH] Bernstein polynomials with no branching. (#591)

---
 src/srf/ratpoly.cpp | 91 ++++++++++-----------------------------------
 src/srf/surface.h   |  4 --
 2 files changed, 19 insertions(+), 76 deletions(-)

diff --git a/src/srf/ratpoly.cpp b/src/srf/ratpoly.cpp
index 70c66b55..f4180cef 100644
--- a/src/srf/ratpoly.cpp
+++ b/src/srf/ratpoly.cpp
@@ -13,84 +13,31 @@
 // and convergence should be fast by now.
 #define RATPOLY_EPS (LENGTH_EPS/(1e2))
 
-double SolveSpace::Bernstein(int k, int deg, double t)
+static double Bernstein(int k, int deg, double t)
 {
-    if(k > deg || k < 0) return 0;
+// indexed by [degree][k][exponent]
+    static const double bernstein_coeff[4][4][4] = {
+    { { 1.0,0.0,0.0,0.0 },  { 1.0,0.0,0.0,0.0 }, { 1.0,0.0,0.0,0.0 }, { 1.0,0.0,0.0,0.0 } },
+    { { 1.0,-1.0,0.0,0.0 }, { 0.0,1.0,0.0,0.0 }, { 0.0,0.0,0.0,0.0 }, { 0.0,0.0,0.0,0.0 } },
+    { { 1.0,-2.0,1.0,0.0 }, { 0.0,2.0,-2.0,0.0 },{ 0.0,0.0,1.0,0.0 }, { 0.0,0.0,0.0,0.0 } },
+    { { 1.0,-3.0,3.0,-1.0 },{ 0.0,3.0,-6.0,3.0 },{ 0.0,0.0,3.0,-3.0}, { 0.0,0.0,0.0,1.0 } } };
 
-    switch(deg) {
-        case 0:
-            return 1;
-
-        case 1:
-            if(k == 0) {
-                return (1 - t);
-            } else if(k == 1) {
-                return t;
-            }
-            break;
-
-        case 2:
-            if(k == 0) {
-                return (1 - t)*(1 - t);
-            } else if(k == 1) {
-                return 2*(1 - t)*t;
-            } else if(k == 2) {
-                return t*t;
-            }
-            break;
-
-        case 3:
-            if(k == 0) {
-                return (1 - t)*(1 - t)*(1 - t);
-            } else if(k == 1) {
-                return 3*(1 - t)*(1 - t)*t;
-            } else if(k == 2) {
-                return 3*(1 - t)*t*t;
-            } else if(k == 3) {
-                return t*t*t;
-            }
-            break;
-    }
-    ssassert(false, "Unexpected degree of spline");
+    const double *c;
+    c = bernstein_coeff[deg][k];
+    return (((c[3]*t+c[2])*t)+c[1])*t+c[0];
 }
 
-double SolveSpace::BernsteinDerivative(int k, int deg, double t)
+static double BernsteinDerivative(int k, int deg, double t)
 {
-    switch(deg) {
-        case 0:
-            return 0;
+    static const double bernstein_derivative_coeff[4][4][3] = {
+    { { 0.0,0.0,0.0 },  { 0.0,0.0,0.0 }, { 0.0,0.0,0.0 }, { 0.0,0.0,0.0 } },
+    { { -1.0,0.0,0.0 }, { 1.0,0.0,0.0 }, { 0.0,0.0,0.0 }, { 0.0,0.0,0.0 } },
+    { { -2.0,2.0,0.0 }, { 2.0,-4.0,0.0 },{ 0.0,2.0,0.0 }, { 0.0,0.0,0.0 } },
+    { { -3.0,6.0,-3.0 },{ 3.0,-12.0,9.0 },{ 0.0,6.0,-9.0}, { 0.0,0.0,3.0 } } };
 
-        case 1:
-            if(k == 0) {
-                return -1;
-            } else if(k == 1) {
-                return 1;
-            }
-            break;
-
-        case 2:
-            if(k == 0) {
-                return -2 + 2*t;
-            } else if(k == 1) {
-                return 2 - 4*t;
-            } else if(k == 2) {
-                return 2*t;
-            }
-            break;
-
-        case 3:
-            if(k == 0) {
-                return -3 + 6*t - 3*t*t;
-            } else if(k == 1) {
-                return 3 - 12*t + 9*t*t;
-            } else if(k == 2) {
-                return 6*t - 9*t*t;
-            } else if(k == 3) {
-                return 3*t*t;
-            }
-            break;
-    }
-    ssassert(false, "Unexpected degree of spline");
+    const double *c;
+    c = bernstein_derivative_coeff[deg][k];
+    return ((c[2]*t)+c[1])*t+c[0];
 }
 
 Vector SBezier::PointAt(double t) const {
diff --git a/src/srf/surface.h b/src/srf/surface.h
index 99418c53..6fabf66a 100644
--- a/src/srf/surface.h
+++ b/src/srf/surface.h
@@ -10,10 +10,6 @@
 #ifndef SOLVESPACE_SURFACE_H
 #define SOLVESPACE_SURFACE_H
 
-// Utility functions, Bernstein polynomials of order 1-3 and their derivatives.
-double Bernstein(int k, int deg, double t);
-double BernsteinDerivative(int k, int deg, double t);
-
 class SBezierList;
 class SSurface;
 class SCurvePt;