diff options
| -rw-r--r-- | libdimension/polynomial.c | 49 |
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]; + } } } |