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 
53 namespace WFMath {
54 
55 bool 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 
63 bool 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 
71 double _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 
82 float _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 
93 CoordType _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 }
WFMath
Generic library namespace.
Definition: shape.h:41
WFMath::CoordType
double CoordType
Basic floating point type.
Definition: const.h:140