path: root/libdimension/polynomial.c

/*************************************************************************
 * Copyright (C) 2009-2010 Tavian Barnes <tavianator@tavianator.com>     *
 *                                                                       *
 * 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                         *
 * <http://www.gnu.org/licenses/>.                                       *
 *************************************************************************/

/**
 * @file
 * Polynomials.
 */

#include "dimension-impl.h"
#include <math.h>

/** 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_normalize_polynomial(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 size_t
dmnsn_eliminate_zero_roots(double poly[], size_t degree)
{
  size_t i;
  for (i = 0; i <= degree && fabs(poly[i]) < dmnsn_epsilon; ++i);

  if (i > 0) {
    for (size_t j = 0; j + i <= degree; ++j) {
      poly[j] = poly[j + i];
    }
  }

  return degree - i;
}

/** Returns the number of sign changes between coefficients of \p poly[]. */
static inline size_t
dmnsn_descartes_rule(const double poly[], size_t degree)
{
  size_t changes = 0;
  int lastsign = dmnsn_signbit(poly[degree]);
  for (size_t i = degree; i-- > 0;) {
    int sign = dmnsn_signbit(poly[i]);
    if (fabs(poly[i]) >= dmnsn_epsilon && sign != lastsign) {
      lastsign = sign;
      ++changes;
    }
  }

  return changes;
}

/** How many levels of Pascal's triangle to precompute. */
#define DMNSN_NBINOM 11

/** Pre-computed values of Pascal's triangle. */
static const double dmnsn_pascals_triangle[DMNSN_NBINOM][DMNSN_NBINOM] = {
  { 1.0 },
  { 1.0, 1.0 },
  { 1.0, 2.0, 1.0 },
  { 1.0, 3.0, 3.0, 1.0 },
  { 1.0, 4.0, 6.0, 4.0, 1.0 },
  { 1.0, 5.0, 10.0, 10.0, 5.0, 1.0 },
  { 1.0, 6.0, 15.0, 20.0, 15.0, 6.0, 1.0 },
  { 1.0, 7.0, 21.0, 35.0, 35.0, 21.0, 7.0, 1.0 },
  { 1.0, 8.0, 28.0, 56.0, 70.0, 56.0, 28.0, 8.0, 1.0 },
  { 1.0, 9.0, 36.0, 84.0, 126.0, 126.0, 84.0, 36.0, 9.0, 1.0 },
  { 1.0, 10.0, 45.0, 120.0, 210.0, 252.0, 210.0, 120.0, 45.0, 10.0, 1.0 }
};

/** Compute the (n k) binomial coefficient. */
static inline 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 ranges that contain a single root, with Uspensky's algorithm. */
static size_t
dmnsn_uspensky_bounds(const double poly[], size_t degree, double bounds[][2],
                      size_t max_roots)
{
  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 {
    /* 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 (size_t i = degree; i-- > 0;) {
      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 + 1], b[degree + 1];
    if (dmnsn_likely(degree < DMNSN_NBINOM)) {
      /* Use precomputed binomial coefficients if possible */
      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) {
          double binom = dmnsn_pascals_triangle[j][i];
          a[i] += binom*poly[j];
          b[i] += binom*poly[degree - j];
        }
      }
    } else {
      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) {
          double binom = dmnsn_binom(j, i);
          a[i] += binom*poly[j];
          b[i] += binom*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) {
      /* 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);
    }
    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) {
      /* Transform the found roots of a[] into roots of poly[] */
      bounds[i][0] += 1.0;
      bounds[i][1] += 1.0;
    }
    n += na;

    return n;
  }
}

/** Calculate a finite upper bound for the roots of the normalized polynomial
    \p poly[]. */
static inline double
dmnsn_root_bound(const 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 += 1.0;
  return bound;
}

/** Improve a root with Newton's method. */
static inline double
dmnsn_improve_root(const double poly[], size_t degree, double x)
{
  double p;
  do {
    /* Calculate the value of the polynomial and its derivative at once */
    p = poly[degree];
    double dp = 0.0;
    for (size_t i = degree; i-- > 0;) {
      dp = dp*x + p;
      p  = p*x  + poly[i];
    }

    x -= p/dp;
  } while (fabs(p) >= fabs(x)*dmnsn_epsilon);

  return x;
}

/** 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)
{
  if (max - min < dmnsn_epsilon) {
    /* Bounds are equal, use Newton's method */
    return dmnsn_improve_root(poly, degree, min);
  }

  double evmin = dmnsn_evaluate_polynomial(poly, degree, min);
  double evmax = dmnsn_evaluate_polynomial(poly, degree, max);
  double initial_evmin = fabs(evmin), initial_evmax = fabs(evmax);
  double mid = 0.0, evmid;
  int lastsign = -1;

  do {
    mid = (min*evmax - max*evmin)/(evmax - evmin);
    evmid = dmnsn_evaluate_polynomial(poly, degree, mid);
    int sign = dmnsn_signbit(evmid);

    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_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
           /* These conditions improve stability when one of the bounds is
              close to a different root than we are trying to find */
           || fabs(evmid) > initial_evmin
           || fabs(evmid) > initial_evmax);

  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;

  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/p, t2 = -t1/2.0;
    x[0] = dmnsn_max(t1, t2);
    x[1] = dmnsn_min(t1, t2);
    if (x[1] >= dmnsn_epsilon)
      return 2;
  }

  if (x[0] >= dmnsn_epsilon)
    return 1;
  else
    return 0;
}

/* Uspensky's algorithm */
DMNSN_HOT 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 (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);
  /* Normalize the leading coefficient to 1.0 */
  dmnsn_normalize_polynomial(p, degree);
  /* Eliminate simple zero roots */
  degree = dmnsn_eliminate_zero_roots(p, degree);

  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_uspensky_bounds(p, degree, ranges, degree - max_algebraic);

    /* 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 */
      ranges[j][1] = dmnsn_min(ranges[j][1], absmax);

      /* 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_print_polynomial(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, "%.15g", fabs(poly[i]));
    if (i >= 2) {
      fprintf(file, "*x^%zu", i);
    } else if (i == 1) {
      fprintf(file, "*x");
    }
  }
}