diff options
| -rw-r--r-- | libdimension/compiler.h | 4 | ||||
| -rw-r--r-- | libdimension/geometry.c | 4 | ||||
| -rw-r--r-- | libdimension/polynomial.c | 31 |
3 files changed, 19 insertions, 20 deletions
diff --git a/libdimension/compiler.h b/libdimension/compiler.h index 1be1c28..b647aea 100644 --- a/libdimension/compiler.h +++ b/libdimension/compiler.h @@ -37,8 +37,8 @@ #define dmnsn_likely(test) __builtin_expect(!!(test), true) #define dmnsn_unlikely(test) __builtin_expect(!!(test), false) #else - #define dmnsn_likely(test) (test) - #define dmnsn_unlikely(test) (test) + #define dmnsn_likely(test) (!!(test)) + #define dmnsn_unlikely(test) (!!(test)) #endif #ifdef __GNUC__ diff --git a/libdimension/geometry.c b/libdimension/geometry.c index 7e2fd77..827e78a 100644 --- a/libdimension/geometry.c +++ b/libdimension/geometry.c @@ -23,7 +23,7 @@ * Geometrical function implementations. */ -#include "dimension.h" +#include "dimension-impl.h" #include <math.h> /* Identity matrix */ @@ -152,7 +152,7 @@ dmnsn_matrix_inverse(dmnsn_matrix A) dmnsn_matrix2 P, Q, R, S, Pi, RPi, PiQ, RPiQ, PP, QQ, RR, SS; double Pdet = A.n[0][0]*A.n[1][1] - A.n[0][1]*A.n[1][0]; - if (fabs(Pdet) < dmnsn_epsilon) { + if (dmnsn_unlikely(fabs(Pdet) < dmnsn_epsilon)) { /* If P is close to singular, try a more generic algorithm; this is very * unlikely, but not impossible, eg. * [ 1 1 0 0 ] diff --git a/libdimension/polynomial.c b/libdimension/polynomial.c index 7a166ac..ed4679b 100644 --- a/libdimension/polynomial.c +++ b/libdimension/polynomial.c @@ -31,7 +31,7 @@ static inline size_t dmnsn_real_degree(const double poly[], size_t degree) { for (size_t i = degree + 1; i-- > 0;) { - if (fabs(poly[i]) >= dmnsn_epsilon) { + if (dmnsn_likely(fabs(poly[i]) >= dmnsn_epsilon)) { return i; } } @@ -50,19 +50,19 @@ dmnsn_normalize_polynomial(double poly[], size_t degree) } /** Eliminate trivial zero roots from \p poly[]. */ -static inline void -dmnsn_eliminate_zero_roots(double poly[], size_t *degree) +static inline size_t +dmnsn_eliminate_zero_roots(double poly[], size_t degree) { - size_t i, deg = *degree; - for (i = 0; i <= deg && fabs(poly[i]) < dmnsn_epsilon; ++i); + size_t i; + for (i = 0; i <= degree && fabs(poly[i]) < dmnsn_epsilon; ++i); if (i > 0) { - for (size_t j = 0; j + i <= deg; ++j) { + for (size_t j = 0; j + i <= degree; ++j) { poly[j] = poly[j + i]; } } - *degree -= i; + return degree - i; } /** Returns the number of sign changes between coefficients of \p poly[]. */ @@ -153,7 +153,7 @@ dmnsn_uspensky_bounds(const double poly[], size_t degree, double bounds[][2], /* 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 (degree < DMNSN_NBINOM) { + if (dmnsn_likely(degree < DMNSN_NBINOM)) { /* Use precomputed binomial coefficients if possible */ for (size_t i = 0; i <= degree; ++i) { a[i] = poly[i]; @@ -294,17 +294,16 @@ dmnsn_bisect_root(const double poly[], size_t degree, double min, double max) } /** Use synthetic division to eliminate the root \p r from \p poly[]. */ -static inline void -dmnsn_eliminate_root(double poly[], size_t *degree, double r) +static inline size_t +dmnsn_eliminate_root(double poly[], size_t degree, double r) { - size_t deg = *degree; - double rem = poly[deg]; - for (size_t i = deg; i-- > 0;) { + double rem = poly[degree]; + for (size_t i = degree; i-- > 0;) { double temp = poly[i]; poly[i] = rem; rem = temp + r*rem; } - --*degree; + return degree - 1; } /** Solve a normalized linear polynomial algebraically. */ @@ -411,7 +410,7 @@ dmnsn_solve_polynomial(const double poly[], size_t degree, double x[]) /* Normalize the leading coefficient to 1.0 */ dmnsn_normalize_polynomial(p, degree); /* Eliminate simple zero roots */ - dmnsn_eliminate_zero_roots(p, °ree); + degree = dmnsn_eliminate_zero_roots(p, degree); static const size_t max_algebraic = 3; if (degree > max_algebraic) { @@ -430,7 +429,7 @@ dmnsn_solve_polynomial(const double poly[], size_t degree, double x[]) double r = dmnsn_bisect_root(p, degree, ranges[j][0], ranges[j][1]); /* Use synthetic division to eliminate the root `r' */ - dmnsn_eliminate_root(p, °ree, r); + degree = dmnsn_eliminate_root(p, degree, r); /* Store the found root */ x[i] = r; |