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/Makefile.am | 1 + libdimension/cylinder.c | 34 ++-- libdimension/dimension.h | 1 + libdimension/dimension/geometry.h | 7 + libdimension/dimension/polynomial.h | 51 ++++++ libdimension/polynomial.c | 313 ++++++++++++++++++++++++++++++++++++ libdimension/sphere.c | 36 ++--- 7 files changed, 411 insertions(+), 32 deletions(-) create mode 100644 libdimension/dimension/polynomial.h create mode 100644 libdimension/polynomial.c 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 #include #include +#include #include #include #include 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 * + * * + * 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 * + * . * + *************************************************************************/ + +/* + * 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 +#include + +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 * + * * + * 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"); +} 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) */ -- cgit v1.2.3