# sajty/wfmath forked from worldforge/wfmath

Merge branch 'epsilon_cleanup'

2 parents 399dad4 + bc0c32c commit b95c87e463c0bd7a2a8803f57a1efb2adff988b1 alriddoch committed Feb 28, 2012
2 TODO
 @@ -1,3 +1,5 @@ +There is a "double dummy" in ball_funcs.h + Inline in the foo.h any method in foo_funcs.h that does not require another header, and has no loop
 @@ -137,7 +137,7 @@ inline void Quaternion::fromAtlas(const AtlasInType& a) CoordType norm = std::sqrt(m_w * m_w + m_vec.sqrMag()); - if (norm <= WFMATH_EPSILON) { + if (norm <= numeric_constants::epsilon()) { m_valid = false; m_vec.setValid(false); return;
 @@ -84,7 +84,7 @@ class AxisBox AxisBox& operator=(const AxisBox& a) {m_low = a.m_low; m_high = a.m_high; return *this;} - bool isEqualTo(const AxisBox& b, CoordType epsilon = WFMATH_EPSILON) const; + bool isEqualTo(const AxisBox& b, CoordType epsilon = numeric_constants::epsilon()) const; bool operator==(const AxisBox& a) const {return isEqualTo(a);} bool operator!=(const AxisBox& a) const {return !isEqualTo(a);}
 @@ -83,7 +83,7 @@ class Ball Ball& operator=(const Ball& b) {m_radius = b.m_radius; m_center = b.m_center; return *this;} - bool isEqualTo(const Ball& b, CoordType epsilon = WFMATH_EPSILON) const; + bool isEqualTo(const Ball& b, CoordType epsilon = numeric_constants::epsilon()) const; bool operator==(const Ball& b) const {return isEqualTo(b);} bool operator!=(const Ball& b) const {return !isEqualTo(b);}
 @@ -75,7 +75,7 @@ Ball BoundingSphere(const container, std::allocator > double dummy; #endif assert("Check that bounding sphere is good to library accuracy" && - m.accuracy(dummy) < WFMATH_EPSILON); + m.accuracy(dummy) < numeric_constants::epsilon()); w = m.center(); Point center;
 @@ -57,6 +57,9 @@ class Quaternion; // Constants +/// Determines how close to machine precision the library tries to come. +#define WFMATH_PRECISION_FUDGE_FACTOR 30 + template struct numeric_constants { @@ -72,6 +75,8 @@ struct numeric_constants static FloatType sqrt3(); /// The natural logarithm of 2 static FloatType log2(); + /// This is the attempted precision of the library. + static FloatType epsilon(); }; template<> @@ -95,6 +100,10 @@ struct numeric_constants static float log2() { return 0.69314718055994530941723212145817656807550013436025F; } + static float epsilon() { + return (WFMATH_PRECISION_FUDGE_FACTOR * + std::numeric_limits::epsilon()); + } }; template<> @@ -118,10 +127,12 @@ struct numeric_constants static double log2() { return 0.69314718055994530941723212145817656807550013436025; } + static double epsilon() { + return (WFMATH_PRECISION_FUDGE_FACTOR * + std::numeric_limits::epsilon()); + } }; -/// Determines how close to machine precision the library tries to come. -#define WFMATH_PRECISION_FUDGE_FACTOR 30 /// How long we can let RotMatrix and Quaternion go before fixing normalization #define WFMATH_MAX_NORM_AGE ((WFMATH_PRECISION_FUDGE_FACTOR * 2) / 3) @@ -135,7 +146,7 @@ typedef float CoordType; double _ScaleEpsilon(double x1, double x2, double epsilon); float _ScaleEpsilon(float x1, float x2, float epsilon); CoordType _ScaleEpsilon(const CoordType* x1, const CoordType* x2, - int length, CoordType epsilon = WFMATH_EPSILON); + int length, CoordType epsilon = numeric_constants::epsilon()); /// Test for equality up to precision epsilon /** @@ -146,12 +157,12 @@ CoordType _ScaleEpsilon(const CoordType* x1, const CoordType* x2, * compare equal, but Equal(0.00010000, 0.00010002, 1.0e-3) will. **/ template -inline bool Equal(const C& c1, const C& c2, CoordType epsilon = WFMATH_EPSILON) +inline bool Equal(const C& c1, const C& c2, CoordType epsilon = numeric_constants::epsilon()) {return c1.isEqualTo(c2, epsilon);} -bool Equal(double x1, double x2, double epsilon = WFMATH_EPSILON); +bool Equal(double x1, double x2, double epsilon = numeric_constants::epsilon()); // Avoid template and expensive casts from float to doubles. -bool Equal(float x1, float x2, float epsilon = WFMATH_EPSILON); +bool Equal(float x1, float x2, float epsilon = numeric_constants::epsilon()); // These let us avoid including for the sake of // std::max() and std::min().
 @@ -57,7 +57,7 @@ void test_general(const C& c) } // We lose precision in string conversion - assert(Equal(c3, c, FloatMax(WFMATH_EPSILON, 1e-5F))); + assert(Equal(c3, c, FloatMax(numeric_constants::epsilon(), 1e-5F))); } } // namespace WFMath
 @@ -31,7 +31,7 @@ static const unsigned ul_max_digits = (unsigned) (8 * sizeof(unsigned long) // number of bits * log_10_of_2 // base 10 vs. base 2 digits + 1 // log(1) == 0, have to add one for leading digit - + WFMATH_EPSILON); // err on the safe side of roundoff + + WFMath::numeric_constants::epsilon()); // err on the safe side of roundoff #else // _MSC_VER static const unsigned ul_max_digits = 10; #endif // _MSC_VER
 @@ -280,7 +280,7 @@ bool Intersect<3>(const RotBox<3>& r, const AxisBox<3>& b, bool proper) assert(false); } - if(axis.sqrMag() < WFMATH_EPSILON * WFMATH_EPSILON) { + if(axis.sqrMag() < numeric_constants::epsilon() * numeric_constants::epsilon()) { // Parallel edges, this reduces to the 2d case above. We've already // checked the bounding box intersections, so we know they intersect. // We don't need to scale WFMATH_EPSILON, det(m_orient) = 1
 @@ -113,7 +113,7 @@ template inline bool Intersect(const Ball& b, const Point& p, bool proper) { return _LessEq(SquaredDistance(b.m_center, p), b.m_radius * b.m_radius - * (1 + WFMATH_EPSILON), proper); + * (1 + numeric_constants::epsilon()), proper); } template @@ -152,7 +152,7 @@ inline bool Contains(const Ball& b, const AxisBox& a, bool proper) sqr_dist += furthest * furthest; } - return _LessEq(sqr_dist, b.m_radius * b.m_radius * (1 + WFMATH_EPSILON), proper); + return _LessEq(sqr_dist, b.m_radius * b.m_radius * (1 + numeric_constants::epsilon()), proper); } template
 @@ -65,7 +65,7 @@ class Line Line& operator=(const Line& a); /// generic: check if two classes are equal, up to a given tolerance - bool isEqualTo(const Line& s, float epsilon = WFMATH_EPSILON) const; + bool isEqualTo(const Line& s, float epsilon = numeric_constants::epsilon()) const; /// generic: check if two classes are equal, up to tolerance WFMATH_EPSILON bool operator==(const Line& s) const {return isEqualTo(s);} /// generic: check if two classes are not equal, up to tolerance WFMATH_EPSILON @@ -89,7 +89,7 @@ class Line Point getCenter() const {return Barycenter(m_points);} // Add before i'th corner, zero is beginning, numCorners() is end - bool addCorner(size_t i, const Point& p, float = WFMATH_EPSILON) + bool addCorner(size_t i, const Point& p, float = numeric_constants::epsilon()) {m_points.insert(m_points.begin() + i, p); return true;} // Remove the i'th corner
 @@ -120,7 +120,7 @@ class Point Point& operator= (const Point& rhs); - bool isEqualTo(const Point &p, CoordType epsilon = WFMATH_EPSILON) const; + bool isEqualTo(const Point &p, CoordType epsilon = numeric_constants::epsilon()) const; bool operator== (const Point& rhs) const {return isEqualTo(rhs);} bool operator!= (const Point& rhs) const {return !isEqualTo(rhs);}
 @@ -190,7 +190,7 @@ Point Barycenter(const container, std::allocator > >& } // Make sure the weights don't add up to zero - if (max_weight <= 0 || std::fabs(tot_weight) <= max_weight * WFMATH_EPSILON) { + if (max_weight <= 0 || std::fabs(tot_weight) <= max_weight * numeric_constants::epsilon()) { return out; }
 @@ -70,7 +70,7 @@ class Polygon<2> Polygon& operator=(const Polygon& p) {m_points = p.m_points; return *this;} - bool isEqualTo(const Polygon& p, CoordType epsilon = WFMATH_EPSILON) const; + bool isEqualTo(const Polygon& p, CoordType epsilon = numeric_constants::epsilon()) const; bool operator==(const Polygon& p) const {return isEqualTo(p);} bool operator!=(const Polygon& p) const {return !isEqualTo(p);} @@ -88,14 +88,14 @@ class Polygon<2> // interface, and the epsilon argument is ignored // Add before i'th corner, zero is beginning, numCorners() is end - bool addCorner(size_t i, const Point<2>& p, double = WFMATH_EPSILON) + bool addCorner(size_t i, const Point<2>& p, CoordType = numeric_constants::epsilon()) {m_points.insert(m_points.begin() + i, p); return true;} // Remove the i'th corner void removeCorner(size_t i) {m_points.erase(m_points.begin() + i);} // Move the i'th corner to p - bool moveCorner(size_t i, const Point<2>& p, double = WFMATH_EPSILON) + bool moveCorner(size_t i, const Point<2>& p, CoordType = numeric_constants::epsilon()) {m_points[i] = p; return true;} // Remove all points @@ -230,7 +230,7 @@ class _Poly2Orient // Try to convert a point from dim dimensions into 2D, expanding the // basis if necessary. Returns true on success. On failure, the // basis is unchanged. - bool expand(const Point& pd, Point<2>& p2, CoordType epsilon = WFMATH_EPSILON); + bool expand(const Point& pd, Point<2>& p2, CoordType epsilon = numeric_constants::epsilon()); // Reduce the basis to the minimum necessary to span the points in // poly (with the exception of skip). Returns _Poly2Reorient, which needs @@ -321,7 +321,7 @@ class Polygon Polygon& operator=(const Polygon& p) {m_orient = p.m_orient; m_poly = p.m_poly; return *this;} - bool isEqualTo(const Polygon& p2, CoordType epsilon = WFMATH_EPSILON) const; + bool isEqualTo(const Polygon& p2, CoordType epsilon = numeric_constants::epsilon()) const; bool operator==(const Polygon& p) const {return isEqualTo(p);} bool operator!=(const Polygon& p) const {return !isEqualTo(p);} @@ -339,7 +339,7 @@ class Polygon // Add before i'th corner, zero is beginning, numCorners() is end // Only succeeds if p lies in a plane with all current points - bool addCorner(size_t i, const Point& p, CoordType epsilon = WFMATH_EPSILON); + bool addCorner(size_t i, const Point& p, CoordType epsilon = numeric_constants::epsilon()); // Remove the i'th corner void removeCorner(size_t i); @@ -348,7 +348,7 @@ class Polygon // lies in the same plane as all the other points. Note that, // under certain circumstances, this plane may not contain the // original location of the point. - bool moveCorner(size_t i, const Point& p, CoordType epsilon = WFMATH_EPSILON); + bool moveCorner(size_t i, const Point& p, CoordType epsilon = numeric_constants::epsilon()); // Remove all points void clear() {m_poly.clear(); m_orient = _Poly2Orient();}
 @@ -162,7 +162,7 @@ _Poly2Reorient _Poly2Orient::reduce(const Polygon<2>& poly, size_t skip) } int exponent; (void) std::frexp(size, &exponent); - epsilon = std::ldexp(WFMATH_EPSILON, exponent); + epsilon = std::ldexp(numeric_constants::epsilon(), exponent); i = 0; if(skip == 0)
 @@ -89,7 +89,7 @@ bool _Poly2Orient<3>::checkIntersectPlane(const AxisBox<3>& b, Point<2>& p2, for(int i = 0; i < 3; ++i) { float foo = std::fabs(normal[i]); - if(std::fabs(normal[i]) < normal_mag * WFMATH_EPSILON) + if(std::fabs(normal[i]) < normal_mag * numeric_constants::epsilon()) axis_direction[i] = AXIS_FLAT; else if(normal[i] > 0) { axis_direction[i] = AXIS_UP; @@ -109,7 +109,7 @@ bool _Poly2Orient<3>::checkIntersectPlane(const AxisBox<3>& b, Point<2>& p2, CoordType perp_size = Dot(normal, high_corner - low_corner) / normal_mag; assert(perp_size >= 0); - if(perp_size < normal_mag * WFMATH_EPSILON) { + if(perp_size < normal_mag * numeric_constants::epsilon()) { // We have a very flat box, lying parallel to the plane return !proper && checkContained(Midpoint(high_corner, low_corner), p2); }
 @@ -68,7 +68,7 @@ inline bool _Poly2Orient::checkContained(const Point& pd, Point<2> & p for(int i = 0; i < dim; ++i) sqrsum += pd[i] * pd[i]; - return off.sqrMag() < WFMATH_EPSILON * sqrsum; + return off.sqrMag() < numeric_constants::epsilon() * sqrsum; } template<> @@ -187,7 +187,7 @@ int _Intersect(const _Poly2Orient &o1, const _Poly2Orient &o2, // Don't need to scale, the m_axes are unit vectors sqrmag1 = basis1.sqrMag(); - if(sqrmag1 > WFMATH_EPSILON * WFMATH_EPSILON) + if(sqrmag1 > numeric_constants::epsilon() * numeric_constants::epsilon()) basis_size = 1; if(o1.m_axes[1].isValid()) { @@ -200,7 +200,7 @@ int _Intersect(const _Poly2Orient &o1, const _Poly2Orient &o2, basis2 -= basis1 * (Dot(basis1, basis2) / sqrmag1); sqrmag2 = basis2.sqrMag(); - if(sqrmag2 > WFMATH_EPSILON * WFMATH_EPSILON) { + if(sqrmag2 > numeric_constants::epsilon() * numeric_constants::epsilon()) { if(basis_size++ == 0) { basis1 = basis2; sqrmag1 = sqrmag2; @@ -329,7 +329,7 @@ int _Intersect(const _Poly2Orient &o1, const _Poly2Orient &o2, data.off[1] = Dot(o2.m_axes[1], off); off_copy -= o1.m_axes[1] * data.off[1]; - if(off_copy.sqrMag() > off_sqr_mag * WFMATH_EPSILON) + if(off_copy.sqrMag() > off_sqr_mag * numeric_constants::epsilon()) return -1; // The planes are different } else
 @@ -214,7 +214,7 @@ Quaternion& Quaternion::rotation(const Vector<3>& axis, CoordType angle) CoordType axis_mag = axis.mag(); CoordType half_angle = angle / 2; - if (axis_mag < WFMATH_EPSILON) { + if (axis_mag < numeric_constants::epsilon()) { m_valid = false; return *this; } @@ -233,7 +233,7 @@ Quaternion& Quaternion::rotation(const Vector<3>& axis) CoordType axis_mag = axis.mag(); CoordType half_angle = axis_mag / 2; - if (axis_mag < WFMATH_EPSILON) { + if (axis_mag < numeric_constants::epsilon()) { m_valid = false; return *this; } @@ -252,13 +252,13 @@ Quaternion& Quaternion::rotation(const Vector<3>& from, const Vector<3>& to) CoordType mag_prod = std::sqrt(from.sqrMag() * to.sqrMag()); CoordType ctheta_plus_1 = Dot(from, to) / mag_prod + 1; - if (mag_prod < WFMATH_EPSILON) { + if (mag_prod < numeric_constants::epsilon()) { m_valid = false; return *this; } // antiparallel vectors - if(ctheta_plus_1 < WFMATH_EPSILON) // same check as used in the RotMatrix function + if(ctheta_plus_1 < numeric_constants::epsilon()) // same check as used in the RotMatrix function throw ColinearVectors<3>(from, to); // cosine of half the angle
 @@ -85,7 +85,7 @@ class Quaternion // This regards q and -1*q as equal, since they give the // same RotMatrix<3> - bool isEqualTo(const Quaternion &q, CoordType epsilon = WFMATH_EPSILON) const; + bool isEqualTo(const Quaternion &q, CoordType epsilon = numeric_constants::epsilon()) const; bool operator== (const Quaternion& rhs) const {return isEqualTo(rhs);} bool operator!= (const Quaternion& rhs) const {return !isEqualTo(rhs);}
 @@ -75,7 +75,7 @@ class RotBox RotBox& operator=(const RotBox& s); - bool isEqualTo(const RotBox& b, CoordType epsilon = WFMATH_EPSILON) const; + bool isEqualTo(const RotBox& b, CoordType epsilon = numeric_constants::epsilon()) const; bool operator==(const RotBox& b) const {return isEqualTo(b);} bool operator!=(const RotBox& b) const {return !isEqualTo(b);}
 @@ -156,7 +156,7 @@ RotMatrix<3>& RotMatrix<3>::rotation (const Vector<3>& axis) } bool _MatrixSetValsImpl(const int size, CoordType* vals, bool& flip, - CoordType* buf1, CoordType* buf2, double precision) + CoordType* buf1, CoordType* buf2, CoordType precision) { precision = std::fabs(precision); @@ -167,7 +167,7 @@ bool _MatrixSetValsImpl(const int size, CoordType* vals, bool& flip, // precision to WFMATH_EPSILON while(true) { - double try_prec = 0; + CoordType try_prec = 0; for(int i = 0; i < size; ++i) { for(int j = 0; j < size; ++j) { @@ -187,7 +187,7 @@ bool _MatrixSetValsImpl(const int size, CoordType* vals, bool& flip, if(try_prec > precision) return false; - if(try_prec <= WFMATH_EPSILON) + if(try_prec <= numeric_constants::epsilon()) break; // Precision needs improvement, use linear approximation scheme.