From 111a5646d12cca087d500554340f03b8aceea54f Mon Sep 17 00:00:00 2001
From: Tavian Barnes <tavianator@gmail.com>
Date: Wed, 3 Nov 2010 02:44:00 -0400
Subject: Pull < DMNSN_NBINOM test outside of loop.

---
 libdimension/polynomial.c | 49 ++++++++++++++++++++++++++++-------------------
 1 file changed, 29 insertions(+), 20 deletions(-)

(limited to 'libdimension/polynomial.c')

diff --git a/libdimension/polynomial.c b/libdimension/polynomial.c
index af8d5df..5dde4e9 100644
--- a/libdimension/polynomial.c
+++ b/libdimension/polynomial.c
@@ -104,21 +104,17 @@ static const double dmnsn_pascals_triangle[DMNSN_NBINOM][DMNSN_NBINOM] = {
 static inline double
 dmnsn_binom(size_t n, size_t k)
 {
-  if (n < DMNSN_NBINOM && k < DMNSN_NBINOM) {
-    return dmnsn_pascals_triangle[n][k];
-  } else {
-    if (k > n - k) {
-      k = n - k;
-    }
-
-    double ret = 1.0;
-    for (size_t i = 0; i < k; ++i) {
-      ret *= n - i;
-      ret /= i + 1;
-    }
+  if (k > n - k) {
+    k = n - k;
+  }
 
-    return ret;
+  double ret = 1.0;
+  for (size_t i = 0; i < k; ++i) {
+    ret *= n - i;
+    ret /= i + 1;
   }
+
+  return ret;
 }
 
 /* Find all ranges that contain a single root, with Uspensky's algorithm */
@@ -157,13 +153,26 @@ dmnsn_uspensky_bounds(const double poly[], size_t degree, double bounds[][2],
     /* a[] is the expanded form of poly(x + 1), b[] is the expanded form of
        (x + 1)^degree * poly(1/(x + 1)) */
     double a[degree + 1], b[degree + 1];
-    for (size_t i = 0; i <= degree; ++i) {
-      a[i] = poly[i];
-      b[i] = poly[degree - i];
-      for (size_t j = i + 1; j <= degree; ++j) {
-        double binom = dmnsn_binom(j, i);
-        a[i] += binom*poly[j];
-        b[i] += binom*poly[degree - j];
+    if (degree < DMNSN_NBINOM) {
+      /* Use precomputed binomial coefficients if possible */
+      for (size_t i = 0; i <= degree; ++i) {
+        a[i] = poly[i];
+        b[i] = poly[degree - i];
+        for (size_t j = i + 1; j <= degree; ++j) {
+          double binom = dmnsn_pascals_triangle[j][i];
+          a[i] += binom*poly[j];
+          b[i] += binom*poly[degree - j];
+        }
+      }
+    } else {
+      for (size_t i = 0; i <= degree; ++i) {
+        a[i] = poly[i];
+        b[i] = poly[degree - i];
+        for (size_t j = i + 1; j <= degree; ++j) {
+          double binom = dmnsn_binom(j, i);
+          a[i] += binom*poly[j];
+          b[i] += binom*poly[degree - j];
+        }
       }
     }
 
-- 
cgit v1.2.3