wfmath 1.0.3
A math library for the Worldforge system.
const.cpp
1// const.cpp (basic comparison functions)
2//
3// The WorldForge Project
4// Copyright (C) 2000, 2001, 2002 The WorldForge Project
5//
6// This program is free software; you can redistribute it and/or modify
7// it under the terms of the GNU General Public License as published by
8// the Free Software Foundation; either version 2 of the License, or
9// (at your option) any later version.
10//
11// This program is distributed in the hope that it will be useful,
12// but WITHOUT ANY WARRANTY; without even the implied warranty of
13// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14// GNU General Public License for more details.
15//
16// You should have received a copy of the GNU General Public License
17// along with this program; if not, write to the Free Software
18// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
19//
20// For information about WorldForge and its authors, please contact
21// the Worldforge Web Site at http://www.worldforge.org.
22//
23
24// Author: Ron Steinke
25
26#ifdef HAVE_CONFIG_H
27#include "config.h"
28#endif
29
30#include <wfmath/const.h>
31
32#include <cmath>
33
34#include <cassert>
35
36// Implementation of IsFloatEqual() is borrowed from Jesse Jones (thanks!).
37// The comments also belong to him. The names have been changed
38// to protect the innocent....
39//
40// Comparing floating point numbers for equality is tricky because
41// the limited precision of the hardware introduces small errors
42// so that two numbers that should compare equal don't. So what
43// we do is consider the numbers equal if their difference is less
44// than some epsilon value. But choosing this epsilon is also tricky
45// because if the numbers are large epsilon should also be large,
46// and if the numbers are very small the epsilon must be even smaller.
47// To get around this we'll use a technique from Knuth's "Seminumerical
48// Algorithms" section 4.2.2 and scale epsilon by the exponent of
49// one of the specified numbers.
50//
51//---------------------------------------------------------------
52
53namespace WFMath {
54
55bool Equal(double x1, double x2, double epsilon)
56{
57 // If the difference between the numbers is smaller than the
58 // scaled epsilon we'll consider the numbers to be equal.
59
60 return std::fabs(x1 - x2) <= _ScaleEpsilon(x1, x2, epsilon);
61}
62
63bool Equal(float x1, float x2, float epsilon)
64{
65 // If the difference between the numbers is smaller than the
66 // scaled epsilon we'll consider the numbers to be equal.
67
68 return std::fabs(x1 - x2) <= _ScaleEpsilon(x1, x2, epsilon);
69}
70
71double _ScaleEpsilon(double x1, double x2, double epsilon)
72{
73 // Get the exponent of the smaller of the two numbers (using the
74 // smaller of the two gives us a tighter epsilon value).
75 int exponent;
76 (void) std::frexp(std::fabs(x1) < std::fabs(x2) ? x1 : x2, &exponent);
77
78 // Scale epsilon by the exponent.
79 return std::ldexp(epsilon, exponent);
80}
81
82float _ScaleEpsilon(float x1, float x2, float epsilon)
83{
84 // Get the exponent of the smaller of the two numbers (using the
85 // smaller of the two gives us a tighter epsilon value).
86 int exponent;
87 (void) std::frexp(std::fabs(x1) < std::fabs(x2) ? x1 : x2, &exponent);
88
89 // Scale epsilon by the exponent.
90 return std::ldexp(epsilon, exponent);
91}
92
93CoordType _ScaleEpsilon(const CoordType* x1, const CoordType* x2,
94 int length, CoordType epsilon)
95{
96 assert(length > 0);
97
98 CoordType max1 = 0, max2 = 0;
99
100 for(int i = 0; i < length; ++i) {
101 auto val1 = std::fabs(x1[i]), val2 = std::fabs(x2[i]);
102 if(val1 > max1)
103 max1 = val1;
104 if(val2 > max2)
105 max2 = val2;
106 }
107
108 return _ScaleEpsilon(max1, max2, epsilon);
109}
110
111}
Generic library namespace.
Definition: shape.h:41
double CoordType
Basic floating point type.
Definition: const.h:140