summaryrefslogtreecommitdiffstats
path: root/libdimension
diff options
context:
space:
mode:
Diffstat (limited to 'libdimension')
-rw-r--r--libdimension/compiler.h4
-rw-r--r--libdimension/geometry.c4
-rw-r--r--libdimension/polynomial.c31
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, &degree);
+ 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, &degree, r);
+ degree = dmnsn_eliminate_root(p, degree, r);
/* Store the found root */
x[i] = r;