summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
Diffstat
-rw-r--r--bench/libdimension/polynomial.c10
-rw-r--r--libdimension/polynomial.c198
-rw-r--r--tests/libdimension/polynomial.c2
3 files changed, 115 insertions, 95 deletions
diff --git a/bench/libdimension/polynomial.c b/bench/libdimension/polynomial.c
index bcf9281..eaf7efe 100644
--- a/bench/libdimension/polynomial.c
+++ b/bench/libdimension/polynomial.c
@@ -27,29 +27,29 @@ main(void)
#define NPOLY 5
double p[NPOLY][NPOLY + 1], x[NPOLY];
- // p[0][] = x - 0.5;
+ /* p[0][] = x - 0.5; */
p[0][1] = 1.0;
p[0][0] = -0.5;
- // p[1][] = (x + 0.5)*(x - 0.5)
+ /* p[1][] = (x + 0.5)*(x - 0.5) */
p[1][2] = 1.0;
p[1][1] = 0.0;
p[1][0] = -0.25;
- // p[2][] = (x + 1)*(x - 1.2345)*(x - 100)
+ /* p[2][] = (x + 1)*(x - 1.2345)*(x - 100) */
p[2][3] = 1.0;
p[2][2] = -100.2345;
p[2][1] = 22.2155;
p[2][0] = 123.45;
- // p[3][] = (x + 1)*(x - 1.2345)*(x - 5)*(x - 100)
+ /* p[3][] = (x + 1)*(x - 1.2345)*(x - 5)*(x - 100) */
p[3][4] = 1.0;
p[3][3] = -105.2345;
p[3][2] = 523.388;
p[3][1] = 12.3725;
p[3][0] = -617.25;
- // p[4][] = (x + 1)*(x - 1.2345)*(x - 2.3456)*(x - 5)*(x - 100)
+ /* p[4][] = (x + 1)*(x - 1.2345)*(x - 2.3456)*(x - 5)*(x - 100) */
p[4][5] = 1.0;
p[4][4] = -107.5801;
p[4][3] = 770.2260432;
diff --git a/libdimension/polynomial.c b/libdimension/polynomial.c
index 14c33d0..6e2aed2 100644
--- a/libdimension/polynomial.c
+++ b/libdimension/polynomial.c
@@ -27,7 +27,7 @@ static inline size_t
dmnsn_solve_linear(const double poly[2], double x[1])
{
x[0] = -poly[0]/poly[1];
- return (x[0] > 0.0) ? 1 : 0;
+ return (x[0] >= dmnsn_epsilon) ? 1 : 0;
}
static inline size_t
@@ -38,34 +38,18 @@ dmnsn_solve_quadratic(const double poly[3], double x[2])
double s = sqrt(disc);
x[0] = (-poly[1] + s)/(2.0*poly[2]);
x[1] = (-poly[1] - s)/(2.0*poly[2]);
- return (x[0] > 0.0) ? ((x[1] > 0.0) ? 2 : 1) : 0;
+ return (x[0] >= dmnsn_epsilon) ? ((x[1] >= dmnsn_epsilon) ? 2 : 1) : 0;
} else {
return 0;
}
}
-/* Eliminate trivial zero roots from poly[] */
-static inline void
-dmnsn_eliminate_zero_roots(double poly[], size_t *degree)
-{
- size_t i, deg = *degree;
- for (i = 0; i <= deg && poly[i] == 0.0; ++i);
-
- if (i > 0) {
- for (size_t j = 0; j + i <= deg; ++j) {
- poly[j] = poly[j + i];
- }
- }
-
- *degree -= i;
-}
-
/* Get the real degree of a polynomial, ignoring leading zeros */
static inline size_t
dmnsn_real_degree(const double poly[], size_t degree)
{
for (ssize_t i = degree; i >= 0; ++i) {
- if (poly[i] != 0.0) {
+ if (fabs(poly[i]) >= dmnsn_epsilon) {
return i;
}
}
@@ -73,56 +57,20 @@ dmnsn_real_degree(const double poly[], size_t degree)
return 0;
}
-/* Use the false position method to find a root in a range that contains exactly
- one root */
-static inline double
-dmnsn_bisect_root(const double poly[], size_t degree, double min, double max)
+/* Eliminate trivial zero roots from poly[] */
+static inline void
+dmnsn_eliminate_zero_roots(double poly[], size_t *degree)
{
- double evmin = dmnsn_evaluate_polynomial(poly, degree, min);
- double evmax = dmnsn_evaluate_polynomial(poly, degree, max);
- double mid = 0.0, evmid;
- int lastsign = -1;
-
- if (max - min <= dmnsn_epsilon)
- return min;
+ size_t i, deg = *degree;
+ for (i = 0; i <= deg && fabs(poly[i]) < dmnsn_epsilon; ++i);
- do {
- mid = (min*evmax - max*evmin)/(evmax - evmin);
- evmid = dmnsn_evaluate_polynomial(poly, degree, mid);
- int sign = dmnsn_signbit(evmid);
- if (sign == dmnsn_signbit(evmax)) {
- max = mid;
- evmax = evmid;
- if (sign == lastsign) {
- /* Don't allow the algorithm to keep the same endpoint for three
- iterations in a row; this ensures superlinear convergence */
- evmin /= 2.0;
- }
- } else {
- min = mid;
- evmin = evmid;
- if (sign == lastsign) {
- evmax /= 2.0;
- }
+ if (i > 0) {
+ for (size_t j = 0; j + i <= deg; ++j) {
+ poly[j] = poly[j + i];
}
-
- lastsign = sign;
- } while (fabs(evmid) > fabs(mid)*dmnsn_epsilon);
-
- return mid;
-}
-
-/* Use synthetic division to eliminate the root `r' from poly[] */
-static inline void
-dmnsn_eliminate_root(double poly[], size_t *degree, double r)
-{
- double rem = poly[*degree];
- for (ssize_t i = *degree - 1; i >= 0; --i) {
- double temp = poly[i];
- poly[i] = rem;
- rem = temp + r*rem;
}
- --*degree;
+
+ *degree -= i;
}
/* Returns the number of sign changes between coefficients of `poly' */
@@ -132,7 +80,7 @@ dmnsn_descartes_rule(const double poly[], size_t degree)
int lastsign = 0;
size_t i;
for (i = 0; i <= degree; ++i) {
- if (poly[i] != 0.0) {
+ if (fabs(poly[i]) >= dmnsn_epsilon) {
lastsign = dmnsn_signbit(poly[i]);
break;
}
@@ -141,7 +89,7 @@ dmnsn_descartes_rule(const double poly[], size_t degree)
size_t changes = 0;
for (++i; i <= degree; ++i) {
int sign = dmnsn_signbit(poly[i]);
- if (poly[i] != 0.0 && sign != lastsign) {
+ if (fabs(poly[i]) >= dmnsn_epsilon && sign != lastsign) {
lastsign = sign;
++changes;
}
@@ -199,11 +147,29 @@ dmnsn_uspensky_bounds(const double poly[], size_t degree, double bounds[][2],
bounds[0][1] = INFINITY;
return 1;
} else {
+ /* Number of roots found so far */
size_t n = 0;
+ /* First divide poly[] by (x - 1) to test for a root at x = 1 */
+ double pdiv1[degree], rem = poly[degree];
+ for (ssize_t i = degree - 1; i >= 0; --i) {
+ pdiv1[i] = rem;
+ rem += poly[i];
+ }
+ if (fabs(rem) < dmnsn_epsilon) {
+ bounds[n][0] = 1.0;
+ bounds[n][1] = 1.0;
+ ++n;
+ if (n == max_roots)
+ return n;
+
+ --degree;
+ poly = pdiv1;
+ }
+
/* 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];
+ double a[degree + 1], b[degree + 1];
for (size_t i = 0; i <= degree; ++i) {
a[i] = poly[i];
b[i] = poly[degree - i];
@@ -214,18 +180,10 @@ dmnsn_uspensky_bounds(const double poly[], size_t degree, double bounds[][2],
}
}
- /* Test for a root at 1.0 */
- if (a[0] == 0.0) {
- bounds[n][0] = 1.0;
- bounds[n][1] = 1.0;
- ++n;
- if (n == max_roots)
- return n;
- }
-
/* 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) {
+ /* Transform the found roots of b[] into roots of poly[] */
double temp = bounds[i][0];
bounds[i][0] = 1.0/(bounds[i][1] + 1.0);
bounds[i][1] = 1.0/(temp + 1.0);
@@ -237,6 +195,7 @@ dmnsn_uspensky_bounds(const double poly[], size_t degree, double bounds[][2],
/* 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) {
+ /* Transform the found roots of a[] into roots of poly[] */
bounds[i][0] += 1.0;
bounds[i][1] += 1.0;
}
@@ -246,35 +205,96 @@ dmnsn_uspensky_bounds(const double poly[], size_t degree, double bounds[][2],
}
}
-/* Modified Uspensky's algorithm */
+/* Calculate a finite upper bound for the roots of poly[] */
+static inline double
+dmnsn_root_bound(double poly[], size_t degree)
+{
+ double bound = 0.0;
+ for (size_t i = 0; i < degree; ++i) {
+ bound = dmnsn_max(bound, fabs(poly[i]));
+ }
+ bound /= fabs(poly[degree]);
+ bound += 1.0;
+ return bound;
+}
+
+/* Use the false position method to find a root in a range that contains exactly
+ one root */
+static inline double
+dmnsn_bisect_root(const double poly[], size_t degree, double min, double max)
+{
+ double evmin = dmnsn_evaluate_polynomial(poly, degree, min);
+ double evmax = dmnsn_evaluate_polynomial(poly, degree, max);
+ double mid = 0.0, evmid;
+ int lastsign = -1;
+
+ if (max - min < dmnsn_epsilon)
+ return min;
+
+ do {
+ mid = (min*evmax - max*evmin)/(evmax - evmin);
+ evmid = dmnsn_evaluate_polynomial(poly, degree, mid);
+ int sign = dmnsn_signbit(evmid);
+ if (sign == dmnsn_signbit(evmax)) {
+ max = mid;
+ evmax = evmid;
+ if (sign == lastsign) {
+ /* Don't allow the algorithm to keep the same endpoint for three
+ iterations in a row; this ensures superlinear convergence */
+ evmin /= 2.0;
+ }
+ } else {
+ min = mid;
+ evmin = evmid;
+ if (sign == lastsign) {
+ evmax /= 2.0;
+ }
+ }
+
+ lastsign = sign;
+ } while (fabs(evmid) >= fabs(mid)*dmnsn_epsilon);
+
+ return mid;
+}
+
+/* Use synthetic division to eliminate the root `r' from poly[] */
+static inline void
+dmnsn_eliminate_root(double poly[], size_t *degree, double r)
+{
+ double rem = poly[*degree];
+ for (ssize_t i = *degree - 1; i >= 0; --i) {
+ double temp = poly[i];
+ poly[i] = rem;
+ rem = temp + r*rem;
+ }
+ --*degree;
+}
+
+/* Uspensky's algorithm */
size_t
dmnsn_solve_polynomial(const double poly[], size_t degree, double x[])
{
/* Copy the polynomial so we can be destructive */
double p[degree + 1];
- for (ssize_t i = degree; i >= 0; --i) {
+ for (size_t i = 0; i <= degree; ++i) {
p[i] = poly[i];
}
+ /* Index into x[] */
+ size_t i = 0;
+
/* Account for leading zero coefficients */
degree = dmnsn_real_degree(p, degree);
/* Eliminate simple zero roots */
dmnsn_eliminate_zero_roots(p, &degree);
- size_t i = 0; /* Index into x[] */
-
if (degree >= 3) {
- /* Find isolating intervals for degree - 2 roots of p[] */
+ /* 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;
+ /* Calculate a finite upper bound for the roots of p[] */
+ double absmax = dmnsn_root_bound(p, degree);
for (size_t j = 0; j < n; ++j) {
/* Replace large or infinite upper bounds with a finite one */
diff --git a/tests/libdimension/polynomial.c b/tests/libdimension/polynomial.c
index 230a1a7..8c6b4da 100644
--- a/tests/libdimension/polynomial.c
+++ b/tests/libdimension/polynomial.c
@@ -29,7 +29,7 @@ int
main()
{
double poly[6], x[5];
- // poly[] = (x + 1)*(x - 1.2345)*(x - 2.3456)*(x - 5)*(x - 100)
+ /* poly[] = (x + 1)*(x - 1.2345)*(x - 2.3456)*(x - 5)*(x - 100) */
poly[5] = 1.0;
poly[4] = -107.5801;
poly[3] = 770.2260432;
generated by cgit v1.2.3 (git 2.39.1) at 2024-11-23 14:39:43 +0000