Add numerical polynomial solver based on Uspensky's algorithm.

7 files changed, 411 insertions, 32 deletions

diff --git a/libdimension/Makefile.am b/libdimension/Makefile.am
index 36517ee..d444682 100644
--- a/libdimension/Makefile.am
+++ b/libdimension/Makefile.am
@@ -74,6 +74,7 @@ libdimension_la_SOURCES = $(nobase_include_HEADERS)                            \
                           platform.c                                           \
                           platform.h                                           \
                           point_light.c                                        \
+                          polynomial.c                                         \
                           progress.c                                           \
                           progress-impl.h                                      \
                           prtree.c                                             \
diff --git a/libdimension/cylinder.c b/libdimension/cylinder.c
index e9db50c..ba016cf 100644
--- a/libdimension/cylinder.c
+++ b/libdimension/cylinder.c
@@ -69,19 +69,27 @@ dmnsn_cylinder_intersection_fn(const dmnsn_object *cylinder, dmnsn_line line,
 
   /* Solve (x0 + nx*t)^2 + (z0 + nz*t)^2
            == (((r2 - r1)*(y0 + ny*t) + r1 + r2)/2)^2 */
-  double a, b, c, t;
-  a = l.n.x*l.n.x + l.n.z*l.n.z - l.n.y*l.n.y*(r2 - r1)*(r2 - r1)/4.0;
-  b = 2.0*(l.n.x*l.x0.x + l.n.z*l.x0.z)
-      - l.n.y*(r2 - r1)*(l.x0.y*(r2 - r1) + r2 + r1)/2.0;
-  c = l.x0.x*l.x0.x + l.x0.z*l.x0.z
-      - (l.x0.y*(r2 - r1) + r2 + r1)*(l.x0.y*(r2 - r1) + r2 + r1)/4;
-
-  if (b*b - 4.0*a*c >= 0.0) {
-    t = (-b - sqrt(b*b - 4.0*a*c))/(2.0*a);
-    dmnsn_vector p = dmnsn_line_point(l, t);
-
-    if (t < 0.0 || p.y <= -1.0 || p.y >= 1.0) {
-      t = (-b + sqrt(b*b - 4.0*a*c))/(2.0*a);
+  double poly[3], x[2];
+  poly[2] = l.n.x*l.n.x + l.n.z*l.n.z - l.n.y*l.n.y*(r2 - r1)*(r2 - r1)/4.0;
+  poly[1] = 2.0*(l.n.x*l.x0.x + l.n.z*l.x0.z)
+            - l.n.y*(r2 - r1)*(l.x0.y*(r2 - r1) + r2 + r1)/2.0;
+  poly[0] = l.x0.x*l.x0.x + l.x0.z*l.x0.z
+            - (l.x0.y*(r2 - r1) + r2 + r1)*(l.x0.y*(r2 - r1) + r2 + r1)/4;
+
+  size_t n = dmnsn_solve_polynomial(poly, 2, x);
+
+  if (n > 0) {
+    double t = x[0];
+    dmnsn_vector p;
+    if (n == 2) {
+      t = dmnsn_min(t, x[1]);
+      p = dmnsn_line_point(l, t);
+
+      if (p.y <= -1.0 || p.y >= 1.0) {
+        t = dmnsn_max(x[0], x[1]);
+        p = dmnsn_line_point(l, t);
+      }
+    } else {
       p = dmnsn_line_point(l, t);
     }
 
diff --git a/libdimension/dimension.h b/libdimension/dimension.h
index c58eb0e..5cdc5b2 100644
--- a/libdimension/dimension.h
+++ b/libdimension/dimension.h
@@ -67,6 +67,7 @@ typedef void dmnsn_free_fn(void *ptr);
 #include <dimension/progress.h>
 #include <dimension/timer.h>
 #include <dimension/geometry.h>
+#include <dimension/polynomial.h>
 #include <dimension/color.h>
 #include <dimension/canvas.h>
 #include <dimension/gl.h>
diff --git a/libdimension/dimension/geometry.h b/libdimension/dimension/geometry.h
index 87e7ab7..4c4f8b8 100644
--- a/libdimension/dimension/geometry.h
+++ b/libdimension/dimension/geometry.h
@@ -96,6 +96,13 @@ dmnsn_degrees(double radians)
   return radians*45.0/atan(1.0);
 }
 
+DMNSN_INLINE int
+dmnsn_signbit(double n)
+{
+  /* Guarantee a 1 or 0 return, to allow testing two signs for equality */
+  return signbit(n) ? 1 : 0;
+}
+
 /* Shorthand for vector/matrix construction */
 
 DMNSN_INLINE dmnsn_vector
diff --git a/libdimension/dimension/polynomial.h b/libdimension/dimension/polynomial.h
new file mode 100644
index 0000000..52a44dd
--- /dev/null
+++ b/libdimension/dimension/polynomial.h
@@ -0,0 +1,51 @@
+/*************************************************************************
+ * Copyright (C) 2009-2010 Tavian Barnes <tavianator@gmail.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/>.                                       *
+ *************************************************************************/
+
+/*
+ * Utility functions for working with and numerically solving polynomials.
+ * Polynomials are represented as simple arrays where the ith element is the
+ * coefficient on x^i.  In general, we are only interested in positive roots.
+ */
+
+#ifndef DIMENSION_POLYNOMIAL_H
+#define DIMENSION_POLYNOMIAL_H
+
+#include <stddef.h>
+#include <stdio.h>
+
+DMNSN_INLINE double
+dmnsn_evaluate_polynomial(double poly[], size_t degree, double x)
+{
+  double ret = poly[degree];
+  ssize_t i;
+  for (i = degree - 1; i >= 0; --i) {
+    ret = ret*x + poly[i];
+  }
+  return ret;
+}
+
+/* Stores the non-negative roots of poly[] in x[], and returns the number of
+   such roots that were stored */
+size_t dmnsn_solve_polynomial(double poly[], size_t degree, double x[]);
+
+/* Helper function to print a polynomial */
+void dmnsn_print_polynomial(FILE *file, double poly[], size_t degree);
+
+#endif /* DIMENSION_POLYNOMIAL_H */
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 <tavianator@gmail.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/>.                                       *
+ *************************************************************************/
+
+#include "dimension.h"
+#include <math.h>
+
+/* 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, &degree, 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, &degree, r);
+
+      /* 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:
+    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");
+}
diff --git a/libdimension/sphere.c b/libdimension/sphere.c
index 98be4b5..490074b 100644
--- a/libdimension/sphere.c
+++ b/libdimension/sphere.c
@@ -53,29 +53,27 @@ dmnsn_sphere_intersection_fn(const dmnsn_object *sphere, dmnsn_line line,
   dmnsn_line l = dmnsn_transform_line(sphere->trans_inv, line);
 
   /* Solve (x0 + nx*t)^2 + (y0 + ny*t)^2 + (z0 + nz*t)^2 == 1 */
-  double a, b, c, t;
-  a = dmnsn_vector_dot(l.n, l.n);
-  b = 2.0*dmnsn_vector_dot(l.n, l.x0);
-  c = dmnsn_vector_dot(l.x0, l.x0) - 1.0;
+  double poly[3], x[2];
+  poly[2] = dmnsn_vector_dot(l.n, l.n);
+  poly[1] = 2.0*dmnsn_vector_dot(l.n, l.x0);
+  poly[0] = dmnsn_vector_dot(l.x0, l.x0) - 1.0;
 
-  if (b*b - 4.0*a*c >= 0.0) {
-    t = (-b - sqrt(b*b - 4.0*a*c))/(2.0*a);
-    if (t < 0.0) {
-      t = (-b + sqrt(b*b - 4.0*a*c))/(2.0*a);
-    }
+  size_t n = dmnsn_solve_polynomial(poly, 2, x);
+  if (n == 0) {
+    return false;
+  } else {
+    double t = x[0];
+    if (n == 2)
+      t = dmnsn_min(t, x[1]);
 
-    if (t >= 0.0) {
-      intersection->ray      = line;
-      intersection->t        = t;
-      intersection->normal   = dmnsn_transform_normal(sphere->trans,
+    intersection->ray      = line;
+    intersection->t        = t;
+    intersection->normal   = dmnsn_transform_normal(sphere->trans,
                                                       dmnsn_line_point(l, t));
-      intersection->texture  = sphere->texture;
-      intersection->interior = sphere->interior;
-      return true;
-    }
+    intersection->texture  = sphere->texture;
+    intersection->interior = sphere->interior;
+    return true;
   }
-
-  return false;
 }
 
 /* Return whether a point is inside a sphere (x**2 + y**2 + z**2 < 1.0) */