summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
-rw-r--r--libdimension/polynomial.c49
1 files changed, 29 insertions, 20 deletions
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];
+ }
}
}