From 7b09710392d35fb55b52031d447a542d99fc6b4b Mon Sep 17 00:00:00 2001 From: Tavian Barnes Date: Tue, 19 Aug 2014 17:10:03 -0400 Subject: Modularize the libdimension codebase. --- libdimension/polynomial.c | 441 ---------------------------------------------- 1 file changed, 441 deletions(-) delete mode 100644 libdimension/polynomial.c (limited to 'libdimension/polynomial.c') diff --git a/libdimension/polynomial.c b/libdimension/polynomial.c deleted file mode 100644 index 609364a..0000000 --- a/libdimension/polynomial.c +++ /dev/null @@ -1,441 +0,0 @@ -/************************************************************************* - * Copyright (C) 2010-2011 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 * - * . * - *************************************************************************/ - -/** - * @file - * Real root isolation algorithm based on work by Vincent, Uspensky, Collins and - * Akritas, Johnson, Krandick, and Rouillier and Zimmerman. - */ - -#include "dimension-internal.h" -#include - -/// Get the real degree of a polynomial, ignoring leading zeros. -static inline size_t -dmnsn_real_degree(const double poly[], size_t degree) -{ - for (size_t i = degree + 1; i-- > 0;) { - if (dmnsn_likely(fabs(poly[i]) >= dmnsn_epsilon)) { - return i; - } - } - - return 0; -} - -/// Divide each coefficient by the leading coefficient. -static inline void -dmnsn_polynomial_normalize(double poly[], size_t degree) -{ - for (size_t i = 0; i < degree; ++i) { - poly[i] /= poly[degree]; - } - poly[degree] = 1.0; -} - -/// Eliminate trivial zero roots from \p poly[]. -static inline void -dmnsn_eliminate_zero_roots(double **poly, size_t *degree) -{ - size_t i; - for (i = 0; i <= *degree; ++i) { - if (dmnsn_likely(fabs((*poly)[i]) >= dmnsn_epsilon)) { - break; - } - } - - *poly += i; - *degree -= i; -} - -/// Calculate a finite upper bound on the roots of a normalized polynomial. -static inline double -dmnsn_root_bound(const double poly[], size_t degree) -{ - double bound = fabs(poly[0]); - for (size_t i = 1; i < degree; ++i) { - bound = dmnsn_max(bound, fabs(poly[i])); - } - bound += 1.0; - return bound; -} - -/// Copy a polynomial. -static inline void -dmnsn_polynomial_copy(double dest[], const double src[], size_t degree) -{ - for (size_t i = 0; i <= degree; ++i) { - dest[i] = src[i]; - } -} - -/// Transform a polynomial by P'(x) = P(x + 1). -static inline void -dmnsn_polynomial_translate(double poly[], size_t degree) -{ - for (size_t i = 0; i <= degree; ++i) { - for (size_t j = degree - i; j <= degree - 1; ++j) { - poly[j] += poly[j + 1]; - } - } -} - -/// Transform a polynomial by P'(x) = P(c*x). -static inline void -dmnsn_polynomial_scale(double poly[], size_t degree, double c) -{ - double factor = c; - for (size_t i = 1; i <= degree; ++i) { - poly[i] *= factor; - factor *= c; - } -} - -/// Returns the result of Descartes' rule on x^degree * poly(1/(x + 1)). -static size_t -dmnsn_descartes_bound(const double poly[], size_t degree) -{ - // Copy the polynomial so we can be destructive - double p[degree + 1]; - dmnsn_polynomial_copy(p, poly, degree); - - // Calculate poly(1/(1/x + 1)) which avoids reversal - for (size_t i = 1; i <= degree; ++i) { - for (size_t j = i; j >= 1; --j) { - p[j] += p[j - 1]; - } - } - - // Find the number of sign changes in p[] - size_t changes = 0; - int lastsign = dmnsn_sgn(p[0]); - for (size_t i = 1; changes <= 1 && i <= degree; ++i) { - int sign = dmnsn_sgn(p[i]); - if (sign != 0 && sign != lastsign) { - ++changes; - lastsign = sign; - } - } - - return changes; -} - -/// Depth-first search of possible isolating intervals. -static size_t -dmnsn_root_bounds_recursive(double poly[], size_t degree, double *c, double *k, - double bounds[][2], size_t nbounds) -{ - size_t s = dmnsn_descartes_bound(poly, degree); - if (s >= 2) { - // Get the left child - dmnsn_polynomial_scale(poly, degree, 1.0/2.0); - *c *= 2.0; - *k /= 2.0; - double currc = *c, currk = *k; - - // Test the left child - size_t n = dmnsn_root_bounds_recursive(poly, degree, c, k, bounds, nbounds); - if (nbounds == n) { - return n; - } - bounds += n; - nbounds -= n; - - // Get the right child from the last tested polynomial - dmnsn_polynomial_translate(poly, degree); - dmnsn_polynomial_scale(poly, degree, currk/(*k)); - *c = currc + 1.0; - *k = currk; - - // Test the right child - n += dmnsn_root_bounds_recursive(poly, degree, c, k, bounds, nbounds); - return n; - } else if (s == 1) { - bounds[0][0] = (*c)*(*k); - bounds[0][1] = (*c + 1.0)*(*k); - return 1; - } else { - return 0; - } -} - -/// Find ranges that contain a single root. -static size_t -dmnsn_root_bounds(const double poly[], size_t degree, double bounds[][2], - size_t nbounds) -{ - // Copy the polynomial so we can be destructive - double p[degree + 1]; - dmnsn_polynomial_copy(p, poly, degree); - - // Scale the roots to within (0, 1] - double bound = dmnsn_root_bound(p, degree); - dmnsn_polynomial_scale(p, degree, bound); - - // Bounding intervals are of the form (c*k, (c + 1)*k) - double c = 0.0, k = 1.0; - - // Isolate the roots - size_t n = dmnsn_root_bounds_recursive(p, degree, &c, &k, bounds, nbounds); - - // Scale the roots back to within (0, bound] - for (size_t i = 0; i < n; ++i) { - bounds[i][0] *= bound; - bounds[i][1] *= bound; - } - - return n; -} - -/// Maximum number of iterations in dmnsn_bisect_root() before bailout. -#define DMNSN_BISECT_ITERATIONS 64 - -/// 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_polynomial_evaluate(poly, degree, min); - double evmax = dmnsn_polynomial_evaluate(poly, degree, max); - - // Handle equal bounds, and equal values at the bounds. - if (dmnsn_unlikely(fabs(evmax - evmin) < dmnsn_epsilon)) { - return (min + max)/2.0; - } - - double evinitial = dmnsn_min(fabs(evmin), fabs(evmax)); - double mid, evmid; - int lastsign = 0; - - for (size_t i = 0; i < DMNSN_BISECT_ITERATIONS; ++i) { - mid = (min*evmax - max*evmin)/(evmax - evmin); - evmid = dmnsn_polynomial_evaluate(poly, degree, mid); - int sign = dmnsn_sgn(evmid); - - if ((fabs(evmid) < fabs(mid)*dmnsn_epsilon - // This condition improves stability when one of the bounds is close to - // a different root than we are trying to find - && fabs(evmid) <= evinitial) - || max - min < fabs(mid)*dmnsn_epsilon) - { - break; - } - - if (mid < min) { - // This can happen due to numerical instability in the root bounding - // algorithm, so behave like the normal secant method - max = min; - evmax = evmin; - min = mid; - evmin = evmid; - } else if (mid > max) { - min = max; - evmin = evmax; - max = mid; - evmax = evmid; - } else if (sign == dmnsn_sgn(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; - } - - return mid; -} - -/// Use synthetic division to eliminate the root \p r from \p poly[]. -static inline size_t -dmnsn_eliminate_root(double poly[], size_t degree, double r) -{ - double rem = poly[degree]; - for (size_t i = degree; i-- > 0;) { - double temp = poly[i]; - poly[i] = rem; - rem = temp + r*rem; - } - return degree - 1; -} - -/// Solve a normalized linear polynomial algebraically. -static inline size_t -dmnsn_solve_linear(const double poly[2], double x[1]) -{ - x[0] = -poly[0]; - if (x[0] >= dmnsn_epsilon) - return 1; - else - return 0; -} - -/// Solve a normalized quadratic polynomial algebraically. -static inline size_t -dmnsn_solve_quadratic(const double poly[3], double x[2]) -{ - double disc = poly[1]*poly[1] - 4.0*poly[0]; - if (disc >= 0.0) { - double s = sqrt(disc); - x[0] = (-poly[1] + s)/2.0; - x[1] = (-poly[1] - s)/2.0; - - if (x[1] >= dmnsn_epsilon) - return 2; - else if (x[0] >= dmnsn_epsilon) - return 1; - else - return 0; - } else { - return 0; - } -} - -/// Solve a normalized cubic polynomial algebraically. -static inline size_t -dmnsn_solve_cubic(double poly[4], double x[3]) -{ - // Reduce to a monic trinomial (t^3 + p*t + q, t = x + b/3) - double b2 = poly[2]*poly[2]; - double p = poly[1] - b2/3.0; - double q = poly[0] - poly[2]*(9.0*poly[1] - 2.0*b2)/27.0; - - double disc = 4.0*p*p*p + 27.0*q*q; - double bdiv3 = poly[2]/3.0; - - if (disc < 0.0) { - // Three real roots -- this implies p < 0 - double msqrtp3 = -sqrt(-p/3.0); - double theta = acos(3*q/(2*p*msqrtp3))/3.0; - - // Store the roots in order from largest to smallest - x[2] = 2.0*msqrtp3*cos(theta) - bdiv3; - x[0] = -2.0*msqrtp3*cos(4.0*atan(1.0)/3.0 - theta) - bdiv3; - x[1] = -(x[0] + x[2] + poly[2]); - - if (x[2] >= dmnsn_epsilon) - return 3; - else if (x[1] >= dmnsn_epsilon) - return 2; - } else if (disc > 0.0) { - // One real root - double cbrtdiscq = cbrt(sqrt(disc/108.0) + fabs(q)/2.0); - double abst = cbrtdiscq - p/(3.0*cbrtdiscq); - - if (q >= 0) { - x[0] = -abst - bdiv3; - } else { - x[0] = abst - bdiv3; - } - } else if (fabs(p) < dmnsn_epsilon) { - // Equation is a perfect cube - x[0] = -bdiv3; - } else { - // Two real roots; one duplicate - double t1 = -(3.0*q)/(2.0*p), t2 = -2.0*t1; - x[0] = dmnsn_max(t1, t2) - bdiv3; - x[1] = dmnsn_min(t1, t2) - bdiv3; - if (x[1] >= dmnsn_epsilon) - return 2; - } - - if (x[0] >= dmnsn_epsilon) - return 1; - else - return 0; -} - -// Solve a polynomial -DMNSN_HOT size_t -dmnsn_polynomial_solve(const double poly[], size_t degree, double x[]) -{ - // Copy the polynomial so we can be destructive - double copy[degree + 1], *p = copy; - dmnsn_polynomial_copy(p, poly, degree); - - // Index into x[] - size_t i = 0; - - // Account for leading zero coefficients - degree = dmnsn_real_degree(p, degree); - // Normalize the leading coefficient to 1.0 - dmnsn_polynomial_normalize(p, degree); - // Eliminate simple zero roots - dmnsn_eliminate_zero_roots(&p, °ree); - - static const size_t max_algebraic = 3; - if (degree > max_algebraic) { - // Find isolating intervals for (degree - max_algebraic) roots of p[] - double ranges[degree - max_algebraic][2]; - size_t n = dmnsn_root_bounds(p, degree, ranges, degree - max_algebraic); - - for (size_t j = 0; j < n; ++j) { - // Bisect within the found range - double r = dmnsn_bisect_root(p, degree, ranges[j][0], ranges[j][1]); - - // Use synthetic division to eliminate the root `r' - degree = dmnsn_eliminate_root(p, degree, r); - - // Store the found root - x[i] = r; - ++i; - } - } - - switch (degree) { - case 1: - i += dmnsn_solve_linear(p, x + i); - break; - case 2: - i += dmnsn_solve_quadratic(p, x + i); - break; - case 3: - i += dmnsn_solve_cubic(p, x + i); - break; - } - - return i; -} - -// Print a polynomial -void -dmnsn_polynomial_print(FILE *file, const double poly[], size_t degree) -{ - for (size_t i = degree + 1; i-- > 0;) { - if (i < degree) { - fprintf(file, (poly[i] >= 0.0) ? " + " : " - "); - } - fprintf(file, "%.17g", fabs(poly[i])); - if (i >= 2) { - fprintf(file, "*x^%zu", i); - } else if (i == 1) { - fprintf(file, "*x"); - } - } -} -- cgit v1.2.3