diff options
| -rw-r--r-- | libdimension/polynomial.c | 120 |
1 files changed, 63 insertions, 57 deletions
diff --git a/libdimension/polynomial.c b/libdimension/polynomial.c index 32dd922..e968832 100644 --- a/libdimension/polynomial.c +++ b/libdimension/polynomial.c @@ -91,7 +91,7 @@ dmnsn_improve_root(double poly[], size_t degree, double x) dx = p/dp; x -= dx; - } while (fabs(dx) > dmnsn_epsilon); + } while (fabs(dx/x) > dmnsn_epsilon); return x; } @@ -174,51 +174,63 @@ dmnsn_binom(size_t n, size_t k) return ret; } -/* Find a range that contains a single root, with Uspensky's algorithm */ -static bool -dmnsn_uspensky_bound(double poly[], size_t degree, double *min, double *max) +/* Find all ranges that contain a single root, with Uspensky's algorithm */ +static size_t +dmnsn_uspensky_bounds(double poly[], size_t degree, double bounds[][2], + size_t max_roots) { - size_t n = dmnsn_descartes_rule(poly, degree); - if (n == 0) { - return false; - } else if (n == 1) { - return true; + size_t signchanges = dmnsn_descartes_rule(poly, degree); + if (signchanges == 0) { + return 0; + } else if (signchanges == 1) { + bounds[0][0] = +0.0; + bounds[0][1] = INFINITY; + return 1; } else { - double a[degree]; - /* a is the expanded form of poly(x + 1) */ + size_t n = 0; + + /* 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], b[degree]; 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) { - a[i] += dmnsn_binom(j, i)*poly[j]; + double binom = dmnsn_binom(j, i); + a[i] += binom*poly[j]; + b[i] += binom*poly[degree - j]; } } + /* Test for a root at 1.0 */ if (a[0] == 0.0) { - *max = *min; - return true; - } else if (dmnsn_uspensky_bound(a, degree, min, max)) { - *min += 1.0; - *max += 1.0; - return true; + bounds[n][0] = 1.0; + bounds[n][1] = 1.0; + ++n; + if (n == max_roots) + return n; } - double b[degree]; - /* b is the expanded form of (x + 1)^degree * poly(1/(x + 1)) */ - for (size_t i = 0; i <= degree; ++i) { - b[i] = poly[degree - i]; - for (size_t j = i + 1; j <= degree; ++j) { - b[i] += dmnsn_binom(j, i)*poly[degree - j]; - } + /* Recursively test for roots in b[] */ + size_t nb = dmnsn_uspensky_bounds(b, degree, bounds + n, max_roots - n); + for (size_t i = n; i < n + nb; ++i) { + double temp = bounds[i][0]; + bounds[i][0] = 1.0/(bounds[i][1] + 1.0); + bounds[i][1] = 1.0/(temp + 1.0); } - - if (dmnsn_uspensky_bound(b, degree, min, max)) { - double temp = *min; - *min = 1.0/(*max + 1.0); - *max = 1.0/(temp + 1.0); - return true; + n += nb; + if (n == max_roots) + return n; + + /* Recursively test for roots in a[] */ + size_t na = dmnsn_uspensky_bounds(a, degree, bounds + n, max_roots - n); + for (size_t i = n; i < n + na; ++i) { + bounds[i][0] += 1.0; + bounds[i][1] += 1.0; } + n += na; - return false; + return n; } } @@ -239,22 +251,25 @@ dmnsn_solve_polynomial(double poly[], size_t degree, double x[]) /* Account for leading zero coefficients */ degree = dmnsn_real_degree(p, degree); - while (degree > 2) { - /* Get a bound on the range of positive roots */ - double min = +0.0, max = INFINITY; - if (dmnsn_uspensky_bound(p, degree, &min, &max)) { - if (isinf(max)) { - /* Replace an infinite upper bound with a finite one due to Cauchy */ - max = 0.0; - for (size_t j = 0; j < degree; ++j) { - max = dmnsn_max(max, fabs(p[j])); - } - max /= fabs(p[degree]); - max += 1.0; - } + if (degree >= 3) { + /* Find isolating intervals for degree - 2 roots of p[] */ + double ranges[degree - 2][2]; + size_t n = dmnsn_uspensky_bounds(p, degree, ranges, degree - 2); + + /* Calculate an finite upper bound for the roots of p[] */ + double absmax = 0.0; + for (size_t j = 0; j < degree; ++j) { + absmax = dmnsn_max(absmax, fabs(p[j])); + } + absmax /= fabs(p[degree]); + absmax += 1.0; + + for (size_t j = 0; j < n; ++j) { + /* Replace large or infinite upper bounds with a finite one */ + ranges[j][1] = dmnsn_min(ranges[j][1], absmax); /* Bisect within the found range */ - double r = dmnsn_bisect_root(p, degree, min, max); + double r = dmnsn_bisect_root(p, degree, ranges[j][0], ranges[j][1]); r = dmnsn_improve_root(p, degree, r); /* Use synthetic division to eliminate the root `r' */ @@ -263,22 +278,13 @@ dmnsn_solve_polynomial(double poly[], size_t degree, double x[]) /* Store the found root */ x[i] = r; ++i; - } else { - break; } - - i += dmnsn_zero_roots(p, °ree, x + i); - degree = dmnsn_real_degree(p, degree); } - switch (degree) { - case 1: + if (degree == 1) { i += dmnsn_solve_linear(p, x + i); - break; - - case 2: + } else if (degree == 2) { i += dmnsn_solve_quadratic(p, x + i); - break; } return i; |