summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
-rw-r--r--libdimension/polynomial.c120
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, &degree, 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;