From 0075e0c37d9c33ed00e4308e6444b61b204327ba Mon Sep 17 00:00:00 2001 From: Tavian Barnes Date: Tue, 26 Oct 2010 15:42:13 -0400 Subject: Add numerical polynomial solver based on Uspensky's algorithm. --- libdimension/polynomial.c | 313 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 313 insertions(+) create mode 100644 libdimension/polynomial.c (limited to 'libdimension/polynomial.c') diff --git a/libdimension/polynomial.c b/libdimension/polynomial.c new file mode 100644 index 0000000..41313bb --- /dev/null +++ b/libdimension/polynomial.c @@ -0,0 +1,313 @@ +/************************************************************************* + * Copyright (C) 2009-2010 Tavian Barnes * + * * + * This file is part of The Dimension Library. * + * * + * The Dimension Library is free software; you can redistribute it and/ * + * or modify it under the terms of the GNU Lesser General Public License * + * as published by the Free Software Foundation; either version 3 of the * + * License, or (at your option) any later version. * + * * + * The Dimension Library is distributed in the hope that it will be * + * useful, but WITHOUT ANY WARRANTY; without even the implied warranty * + * of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * + * Lesser General Public License for more details. * + * * + * You should have received a copy of the GNU Lesser General Public * + * License along with this program. If not, see * + * . * + *************************************************************************/ + +#include "dimension.h" +#include + +/* Basic solving methods */ + +static inline size_t +dmnsn_solve_linear(double poly[2], double x[1]) +{ + x[0] = -poly[0]/poly[1]; + return (x[0] >= 0.0) ? 1 : 0; +} + +static inline size_t +dmnsn_solve_quadratic(double poly[3], double x[2]) +{ + double disc = poly[1]*poly[1] - 4.0*poly[0]*poly[2]; + if (disc >= 0.0) { + 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; + } else { + return 0; + } +} + +static inline size_t +dmnsn_zero_roots(double poly[], size_t *degree, double x[]) +{ + size_t i = 0; + while (i <= *degree && poly[i] == 0.0) { + x[i] = 0.0; + ++i; + } + + if (i > 0) { + for (size_t j = 0; j + i <= *degree; ++j) { + poly[j] = poly[j + i]; + } + } + + *degree -= i; + return i; +} + +/* Get the real degree of a polynomial, ignoring leading zeros */ +static inline size_t +dmnsn_real_degree(double poly[], size_t degree) +{ + for (ssize_t i = degree; i >= 0; ++i) { + if (poly[i] != 0.0) { + return i; + } + } + + return 0; +} + +/* Improve a root with Newton's method */ +static inline double +dmnsn_improve_root(double poly[], size_t degree, double x) +{ + double error; + do { + /* Calculate the value of the polynomial and its derrivative at once */ + double p = poly[degree], dp = 0.0; + for (ssize_t i = degree - 1; i >= 0; --i) { + dp = dp*x + p; + p = p*x + poly[i]; + } + + double dx = p/dp; + error = fabs(dx/x); + x -= dx; + } while (error > dmnsn_epsilon); + + return x; +} + +/* Use the method of bisection to find a root in a range that contains exactly + one root, counting multiplicity */ +static inline double +dmnsn_bisect_root(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); + + while (max - min > dmnsn_epsilon) { + double mid = (min + max)/2.0; + double evmid = dmnsn_evaluate_polynomial(poly, degree, mid); + if (dmnsn_signbit(evmid) == dmnsn_signbit(evmax)) { + max = mid; + evmax = evmid; + } else { + min = mid; + evmin = evmid; + } + } + + return (min + max)/2.0; +} + +/* 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; +} + +/* Returns the number of sign changes between coefficients of `poly' */ +static inline size_t +dmnsn_descartes_rule(double poly[], size_t degree) +{ + int lastsign = 0; + size_t i; + for (i = 0; i <= degree; ++i) { + if (poly[i] != 0.0) { + lastsign = dmnsn_signbit(poly[i]); + break; + } + } + + size_t changes = 0; + for (++i; i <= degree; ++i) { + int sign = dmnsn_signbit(poly[i]); + if (poly[i] != 0.0 && sign != lastsign) { + lastsign = sign; + ++changes; + } + } + + return changes; +} + +/* Get the (n k) binomial coefficient */ +static double +dmnsn_binom(size_t n, size_t k) +{ + 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; + } + + 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) +{ + size_t n = dmnsn_descartes_rule(poly, degree); + if (n == 0) { + return false; + } else if (n == 1) { + return true; + } else { + double a[degree]; + /* a is the expanded form of poly(x + 1) */ + for (size_t i = 0; i <= degree; ++i) { + a[i] = poly[i]; + for (size_t j = i + 1; j <= degree; ++j) { + a[i] += dmnsn_binom(j, i)*poly[j]; + } + } + + 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; + } + + 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]; + } + } + + 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; + } + + return false; + } +} + +/* Modified Uspensky's algorithm */ +size_t +dmnsn_solve_polynomial(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) { + p[i] = poly[i]; + } + + size_t i = 0; /* Index into x[] */ + + /* Eliminate simple zero roots */ + i += dmnsn_zero_roots(p, °ree, x + i); + /* 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; + } + + /* Bisect within the found range */ + double r = dmnsn_bisect_root(p, degree, min, max); + r = dmnsn_improve_root(p, degree, r); + + /* Use synthetic division to eliminate the root `r' */ + dmnsn_eliminate_root(p, °ree, r); + + /* 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: + i += dmnsn_solve_linear(p, x + i); + break; + + case 2: + i += dmnsn_solve_quadratic(p, x + i); + break; + } + + return i; +} + +void +dmnsn_print_polynomial(FILE *file, double poly[], size_t degree) +{ + fprintf(file, "%g*x^%zu", poly[degree], degree); + for (size_t i = degree - 1; i > 1; --i) { + if (poly[i] > 0.0) { + fprintf(file, " + %g*x^%zu", poly[i], i); + } else if (poly[i] < 0.0) { + fprintf(file, " - %g*x^%zu", -poly[i], i); + } + } + + if (poly[1] > 0.0) { + fprintf(file, " + %g*x", poly[1]); + } else if (poly[1] < 0.0) { + fprintf(file, " - %g*x", -poly[1]); + } + + if (poly[0] > 0.0) { + fprintf(file, " + %g", poly[0]); + } else if (poly[0] < 0.0) { + fprintf(file, " - %g", -poly[0]); + } + + fprintf(file, "\n"); +} -- cgit v1.2.3