dune-common 3.0-git
float_cmp.cc
Go to the documentation of this file.
1// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2// vi: set et ts=4 sw=2 sts=2:
3#include "float_cmp.hh"
4
5#include <vector>
6#include <limits>
7#include <algorithm>
8#include <cstdlib>
10
11namespace Dune {
12
13
14 namespace FloatCmp {
15 // traits
17
21 template<class T> struct EpsilonType {
23 typedef T Type;
24 };
26
31 template<class T, typename A>
32 struct EpsilonType<std::vector<T, A> > {
35 };
37
42 template<class T, int n>
43 struct EpsilonType<FieldVector<T, n> > {
46 };
47
48 // default epsilon
49 template<class T>
51 static typename EpsilonType<T>::Type value()
52 { return std::numeric_limits<typename EpsilonType<T>::Type>::epsilon()*8; }
53 };
54 template<class T>
56 static typename EpsilonType<T>::Type value()
57 { return std::numeric_limits<typename EpsilonType<T>::Type>::epsilon()*8; }
58 };
59 template<class T>
61 static typename EpsilonType<T>::Type value()
62 { return std::max(std::numeric_limits<typename EpsilonType<T>::Type>::epsilon(), 1e-6); }
63 };
64
65 namespace Detail {
66 // basic comparison
67 template<class T, CmpStyle style = defaultCmpStyle>
68 struct eq_t;
69 template<class T>
70 struct eq_t<T, relativeWeak> {
71 static bool eq(const T &first,
72 const T &second,
74 { return std::abs(first - second) <= epsilon*std::max(std::abs(first), std::abs(second)); }
75 };
76 template<class T>
77 struct eq_t<T, relativeStrong> {
78 static bool eq(const T &first,
79 const T &second,
81 { return std::abs(first - second) <= epsilon*std::min(std::abs(first), std::abs(second)); }
82 };
83 template<class T>
84 struct eq_t<T, absolute> {
85 static bool eq(const T &first,
86 const T &second,
88 { return std::abs(first-second) <= epsilon; }
89 };
90 template<class T, CmpStyle cstyle>
91 struct eq_t<std::vector<T>, cstyle> {
92 static bool eq(const std::vector<T> &first,
93 const std::vector<T> &second,
95 unsigned int size = first.size();
96 if(size != second.size()) return false;
97 for(unsigned int i = 0; i < size; ++i)
98 if(!eq_t<T, cstyle>(first[i], second[i], epsilon))
99 return false;
100 return true;
101 }
102 };
103 template<class T, int n, CmpStyle cstyle>
104 struct eq_t<Dune::FieldVector<T, n>, cstyle> {
105 static bool eq(const Dune::FieldVector<T, n> &first,
106 const Dune::FieldVector<T, n> &second,
107 typename EpsilonType<T>::Type epsilon = DefaultEpsilon<T>::value()) {
108 for(int i = 0; i < n; ++i)
109 if(!eq_t<T, cstyle>(first[i], second[i], epsilon))
110 return false;
111 return true;
112 }
113 };
114 } // namespace Detail
115
116 // operations in functional style
117 template <class T, CmpStyle style>
118 bool eq(const T &first,
119 const T &second,
120 typename EpsilonType<T>::Type epsilon)
121 {
122 return Detail::eq_t<T, style>::eq(first, second, epsilon);
123 }
124 template <class T, CmpStyle style>
125 bool ne(const T &first,
126 const T &second,
127 typename EpsilonType<T>::Type epsilon)
128 {
129 return !eq<T, style>(first, second, epsilon);
130 }
131 template <class T, CmpStyle style>
132 bool gt(const T &first,
133 const T &second,
134 typename EpsilonType<T>::Type epsilon)
135 {
136 return first > second && ne<T, style>(first, second, epsilon);
137 }
138 template <class T, CmpStyle style>
139 bool lt(const T &first,
140 const T &second,
141 typename EpsilonType<T>::Type epsilon)
142 {
143 return first < second && ne<T, style>(first, second, epsilon);
144 }
145 template <class T, CmpStyle style>
146 bool ge(const T &first,
147 const T &second,
148 typename EpsilonType<T>::Type epsilon)
149 {
150 return first > second || eq<T, style>(first, second, epsilon);
151 }
152 template <class T, CmpStyle style>
153 bool le(const T &first,
154 const T &second,
155 typename EpsilonType<T>::Type epsilon)
156 {
157 return first < second || eq<T, style>(first, second, epsilon);
158 }
159
160 // default template arguments
161 template <class T>
162 bool eq(const T &first,
163 const T &second,
165 {
166 return eq<T, defaultCmpStyle>(first, second, epsilon);
167 }
168 template <class T>
169 bool ne(const T &first,
170 const T &second,
172 {
173 return ne<T, defaultCmpStyle>(first, second, epsilon);
174 }
175 template <class T>
176 bool gt(const T &first,
177 const T &second,
179 {
180 return gt<T, defaultCmpStyle>(first, second, epsilon);
181 }
182 template <class T>
183 bool lt(const T &first,
184 const T &second,
186 {
187 return lt<T, defaultCmpStyle>(first, second, epsilon);
188 }
189 template <class T>
190 bool ge(const T &first,
191 const T &second,
193 {
194 return ge<T, defaultCmpStyle>(first, second, epsilon);
195 }
196 template <class T>
197 bool le(const T &first,
198 const T &second,
200 {
201 return le<T, defaultCmpStyle>(first, second, epsilon);
202 }
203
204 // rounding operations
205 namespace Detail {
206 template<class I, class T, CmpStyle cstyle = defaultCmpStyle, RoundingStyle rstyle = defaultRoundingStyle>
207 struct round_t;
208 template<class I, class T, CmpStyle cstyle>
209 struct round_t<I, T, cstyle, downward> {
210 static I
211 round(const T &val,
213 // first get an approximation
214 I lower = I(val);
215 I upper;
216 if(eq<T, cstyle>(T(lower), val, epsilon)) return lower;
217 if(T(lower) > val) { upper = lower; lower--; }
218 else upper = lower+1;
219 if(le<T, cstyle>(val - T(lower), T(upper) - val, epsilon))
220 return lower;
221 else return upper;
222 }
223 };
224 template<class I, class T, CmpStyle cstyle>
225 struct round_t<I, T, cstyle, upward> {
226 static I
227 round(const T &val,
229 // first get an approximation
230 I lower = I(val);
231 I upper;
232 if(eq<T, cstyle>(T(lower), val, epsilon)) return lower;
233 if(T(lower) > val) { upper = lower; lower--; }
234 else upper = lower+1;
235 if(lt<T, cstyle>(val - T(lower), T(upper) - val, epsilon))
236 return lower;
237 else return upper;
238 }
239 };
240 template<class I, class T, CmpStyle cstyle>
241 struct round_t<I, T, cstyle, towardZero> {
242 static I
243 round(const T &val,
245 if(val > T(0))
246 return round_t<I, T, cstyle, downward>::round(val, epsilon);
247 else return round_t<I, T, cstyle, upward>::round(val, epsilon);
248 }
249 };
250 template<class I, class T, CmpStyle cstyle>
251 struct round_t<I, T, cstyle, towardInf> {
252 static I
253 round(const T &val,
255 if(val > T(0))
256 return round_t<I, T, cstyle, upward>::round(val, epsilon);
257 else return round_t<I, T, cstyle, downward>::round(val, epsilon);
258 }
259 };
260 template<class I, class T, CmpStyle cstyle, RoundingStyle rstyle>
261 struct round_t<std::vector<I>, std::vector<T>, cstyle, rstyle> {
262 static std::vector<I>
263 round(const T &val,
265 unsigned int size = val.size();
266 std::vector<I> res(size);
267 for(unsigned int i = 0; i < size; ++i)
268 res[i] = round_t<I, T, cstyle, rstyle>::round(val[i], epsilon);
269 return res;
270 }
271 };
272 template<class I, class T, int n, CmpStyle cstyle, RoundingStyle rstyle>
273 struct round_t<Dune::FieldVector<I, n>, Dune::FieldVector<T, n>, cstyle, rstyle> {
275 round(const T &val,
278 for(int i = 0; i < n; ++i)
279 res[i] = round_t<I, T, cstyle, rstyle>::round(val[i], epsilon);
280 return res;
281 }
282 };
283 } // namespace Detail
284 template<class I, class T, CmpStyle cstyle, RoundingStyle rstyle>
285 I round(const T &val, typename EpsilonType<T>::Type epsilon /*= DefaultEpsilon<T, cstyle>::value()*/)
286 {
288 }
289 template<class I, class T, CmpStyle cstyle>
290 I round(const T &val, typename EpsilonType<T>::Type epsilon = DefaultEpsilon<T, cstyle>::value())
291 {
292 return round<I, T, cstyle, defaultRoundingStyle>(val, epsilon);
293 }
294 template<class I, class T, RoundingStyle rstyle>
295 I round(const T &val, typename EpsilonType<T>::Type epsilon = DefaultEpsilon<T, defaultCmpStyle>::value())
296 {
297 return round<I, T, defaultCmpStyle, rstyle>(val, epsilon);
298 }
299 template<class I, class T>
301 {
302 return round<I, T, defaultCmpStyle>(val, epsilon);
303 }
304
305 // truncation
306 namespace Detail {
307 template<class I, class T, CmpStyle cstyle = defaultCmpStyle, RoundingStyle rstyle = defaultRoundingStyle>
308 struct trunc_t;
309 template<class I, class T, CmpStyle cstyle>
310 struct trunc_t<I, T, cstyle, downward> {
311 static I
312 trunc(const T &val,
314 // this should be optimized away unless needed
315 if(!std::numeric_limits<I>::is_signed)
316 // make sure this works for all useful cases even if I is an unsigned type
317 if(eq<T, cstyle>(val, T(0), epsilon)) return I(0);
318 // first get an approximation
319 I lower = I(val); // now |val-lower| < 1
320 // make sure we're really lower in case the cast truncated to an unexpected direction
321 if(T(lower) > val) lower--; // now val-lower < 1
322 // check whether lower + 1 is approximately val
323 if(eq<T, cstyle>(T(lower+1), val, epsilon))
324 return lower+1;
325 else return lower;
326 }
327 };
328 template<class I, class T, CmpStyle cstyle>
329 struct trunc_t<I, T, cstyle, upward> {
330 static I
331 trunc(const T &val,
333 I upper = trunc_t<I, T, cstyle, downward>::trunc(val, epsilon);
334 if(ne<T, cstyle>(T(upper), val, epsilon)) ++upper;
335 return upper;
336 }
337 };
338 template<class I, class T, CmpStyle cstyle>
339 struct trunc_t<I, T, cstyle, towardZero> {
340 static I
341 trunc(const T &val,
343 if(val > T(0)) return trunc_t<I, T, cstyle, downward>::trunc(val, epsilon);
344 else return trunc_t<I, T, cstyle, upward>::trunc(val, epsilon);
345 }
346 };
347 template<class I, class T, CmpStyle cstyle>
348 struct trunc_t<I, T, cstyle, towardInf> {
349 static I
350 trunc(const T &val,
352 if(val > T(0)) return trunc_t<I, T, cstyle, upward>::trunc(val, epsilon);
353 else return trunc_t<I, T, cstyle, downward>::trunc(val, epsilon);
354 }
355 };
356 template<class I, class T, CmpStyle cstyle, RoundingStyle rstyle>
357 struct trunc_t<std::vector<I>, std::vector<T>, cstyle, rstyle> {
358 static std::vector<I>
359 trunc(const std::vector<T> &val,
361 unsigned int size = val.size();
362 std::vector<I> res(size);
363 for(unsigned int i = 0; i < size; ++i)
364 res[i] = trunc_t<I, T, cstyle, rstyle>::trunc(val[i], epsilon);
365 return res;
366 }
367 };
368 template<class I, class T, int n, CmpStyle cstyle, RoundingStyle rstyle>
369 struct trunc_t<Dune::FieldVector<I, n>, Dune::FieldVector<T, n>, cstyle, rstyle> {
374 for(int i = 0; i < n; ++i)
375 res[i] = trunc_t<I, T, cstyle, rstyle>::trunc(val[i], epsilon);
376 return res;
377 }
378 };
379 } // namespace Detail
380 template<class I, class T, CmpStyle cstyle, RoundingStyle rstyle>
381 I trunc(const T &val, typename EpsilonType<T>::Type epsilon /*= DefaultEpsilon<T, cstyle>::value()*/)
382 {
384 }
385 template<class I, class T, CmpStyle cstyle>
386 I trunc(const T &val, typename EpsilonType<T>::Type epsilon = DefaultEpsilon<T, cstyle>::value())
387 {
388 return trunc<I, T, cstyle, defaultRoundingStyle>(val, epsilon);
389 }
390 template<class I, class T, RoundingStyle rstyle>
391 I trunc(const T &val, typename EpsilonType<T>::Type epsilon = DefaultEpsilon<T, defaultCmpStyle>::value())
392 {
393 return trunc<I, T, defaultCmpStyle, rstyle>(val, epsilon);
394 }
395 template<class I, class T>
397 {
398 return trunc<I, T, defaultCmpStyle>(val, epsilon);
399 }
400 } //namespace Dune
401
402 // oo interface
403 template<class T, FloatCmp::CmpStyle cstyle_, FloatCmp::RoundingStyle rstyle_>
406
407
408 template<class T, FloatCmp::CmpStyle cstyle_, FloatCmp::RoundingStyle rstyle_>
411 {
412 return epsilon_;
413 }
414
415 template<class T, FloatCmp::CmpStyle cstyle_, FloatCmp::RoundingStyle rstyle_>
416 void
418 {
419 epsilon_ = epsilon__;
420 }
421
422
423 template<class T, FloatCmp::CmpStyle cstyle_, FloatCmp::RoundingStyle rstyle_>
425 eq(const ValueType &first, const ValueType &second) const
426 {
427 return Dune::FloatCmp::eq<ValueType, cstyle>(first, second, epsilon_);
428 }
429
430 template<class T, FloatCmp::CmpStyle cstyle_, FloatCmp::RoundingStyle rstyle_>
432 ne(const ValueType &first, const ValueType &second) const
433 {
434 return Dune::FloatCmp::ne<ValueType, cstyle>(first, second, epsilon_);
435 }
436
437 template<class T, FloatCmp::CmpStyle cstyle_, FloatCmp::RoundingStyle rstyle_>
439 gt(const ValueType &first, const ValueType &second) const
440 {
441 return Dune::FloatCmp::gt<ValueType, cstyle>(first, second, epsilon_);
442 }
443
444 template<class T, FloatCmp::CmpStyle cstyle_, FloatCmp::RoundingStyle rstyle_>
446 lt(const ValueType &first, const ValueType &second) const
447 {
448 return Dune::FloatCmp::lt<ValueType, cstyle>(first, second, epsilon_);
449 }
450
451 template<class T, FloatCmp::CmpStyle cstyle_, FloatCmp::RoundingStyle rstyle_>
453 ge(const ValueType &first, const ValueType &second) const
454 {
455 return Dune::FloatCmp::ge<ValueType, cstyle>(first, second, epsilon_);
456 }
457
458 template<class T, FloatCmp::CmpStyle cstyle_, FloatCmp::RoundingStyle rstyle_>
460 le(const ValueType &first, const ValueType &second) const
461 {
462 return Dune::FloatCmp::le<ValueType, cstyle>(first, second, epsilon_);
463 }
464
465
466 template<class T, FloatCmp::CmpStyle cstyle_, FloatCmp::RoundingStyle rstyle_>
467 template<class I>
469 round(const ValueType &val) const
470 {
471 return Dune::FloatCmp::round<I, ValueType, cstyle, rstyle_>(val, epsilon_);
472 }
473
474 template<class T, FloatCmp::CmpStyle cstyle_, FloatCmp::RoundingStyle rstyle_>
475 template<class I>
477 trunc(const ValueType &val) const
478 {
479 return Dune::FloatCmp::trunc<I, ValueType, cstyle, rstyle_>(val, epsilon_);
480 }
481
482} //namespace Dune
Implements a vector constructed from a given type representing a field and a compile-time given size.
Various ways to compare floating-point numbers.
bool ne(const T &first, const T &second, typename EpsilonType< T >::Type epsilon)
test for inequality using epsilon
Definition float_cmp.cc:125
bool eq(const T &first, const T &second, typename EpsilonType< T >::Type epsilon)
test for equality using epsilon
Definition float_cmp.cc:118
I round(const T &val, typename EpsilonType< T >::Type epsilon)
round using epsilon
Definition float_cmp.cc:285
I trunc(const T &val, typename EpsilonType< T >::Type epsilon)
truncate using epsilon
Definition float_cmp.cc:381
bool lt(const T &first, const T &second, typename EpsilonType< T >::Type epsilon)
test if first lesser than second
Definition float_cmp.cc:139
bool gt(const T &first, const T &second, typename EpsilonType< T >::Type epsilon)
test if first greater than second
Definition float_cmp.cc:132
bool ge(const T &first, const T &second, typename EpsilonType< T >::Type epsilon)
test if first greater or equal second
Definition float_cmp.cc:146
bool le(const T &first, const T &second, typename EpsilonType< T >::Type epsilon)
test if first lesser or equal second
Definition float_cmp.cc:153
@ relativeStrong
|a-b|/|a| <= epsilon && |a-b|/|b| <= epsilon
Definition float_cmp.hh:106
@ relativeWeak
|a-b|/|a| <= epsilon || |a-b|/|b| <= epsilon
Definition float_cmp.hh:104
@ absolute
|a-b| <= epsilon
Definition float_cmp.hh:108
@ towardZero
always round toward 0
Definition float_cmp.hh:116
@ towardInf
always round away from 0
Definition float_cmp.hh:118
@ upward
round toward
Definition float_cmp.hh:122
@ downward
round toward
Definition float_cmp.hh:120
STL namespace.
Dune namespace.
Definition alignment.hh:11
vector space out of a tensor product of fields.
Definition fvector.hh:93
Mapping of value type to epsilon type.
Definition float_cmp.cc:21
T Type
The epsilon type corresponding to value type T.
Definition float_cmp.cc:23
EpsilonType< T > Type
The epsilon type corresponding to value type std::vector<T, A>
Definition float_cmp.cc:34
EpsilonType< T > Type
The epsilon type corresponding to value type Dune::FieldVector<T, n>
Definition float_cmp.cc:45
static EpsilonType< T >::Type value()
Definition float_cmp.cc:51
static EpsilonType< T >::Type value()
Definition float_cmp.cc:56
static EpsilonType< T >::Type value()
Definition float_cmp.cc:61
Definition float_cmp.cc:68
static bool eq(const T &first, const T &second, typename EpsilonType< T >::Type epsilon=DefaultEpsilon< T >::value())
Definition float_cmp.cc:71
static bool eq(const T &first, const T &second, typename EpsilonType< T >::Type epsilon=DefaultEpsilon< T >::value())
Definition float_cmp.cc:78
static bool eq(const T &first, const T &second, typename EpsilonType< T >::Type epsilon=DefaultEpsilon< T >::value())
Definition float_cmp.cc:85
static bool eq(const std::vector< T > &first, const std::vector< T > &second, typename EpsilonType< T >::Type epsilon=DefaultEpsilon< T >::value())
Definition float_cmp.cc:92
static bool eq(const Dune::FieldVector< T, n > &first, const Dune::FieldVector< T, n > &second, typename EpsilonType< T >::Type epsilon=DefaultEpsilon< T >::value())
Definition float_cmp.cc:105
Definition float_cmp.cc:207
static I round(const T &val, typename EpsilonType< T >::Type epsilon=(DefaultEpsilon< T, cstyle >::value()))
Definition float_cmp.cc:211
static I round(const T &val, typename EpsilonType< T >::Type epsilon=(DefaultEpsilon< T, cstyle >::value()))
Definition float_cmp.cc:227
static I round(const T &val, typename EpsilonType< T >::Type epsilon=(DefaultEpsilon< T, cstyle >::value()))
Definition float_cmp.cc:243
static I round(const T &val, typename EpsilonType< T >::Type epsilon=(DefaultEpsilon< T, cstyle >::value()))
Definition float_cmp.cc:253
static std::vector< I > round(const T &val, typename EpsilonType< T >::Type epsilon=(DefaultEpsilon< T, cstyle >::value()))
Definition float_cmp.cc:263
static Dune::FieldVector< I, n > round(const T &val, typename EpsilonType< T >::Type epsilon=(DefaultEpsilon< T, cstyle >::value()))
Definition float_cmp.cc:275
Definition float_cmp.cc:308
static I trunc(const T &val, typename EpsilonType< T >::Type epsilon=(DefaultEpsilon< T, cstyle >::value()))
Definition float_cmp.cc:312
static I trunc(const T &val, typename EpsilonType< T >::Type epsilon=(DefaultEpsilon< T, cstyle >::value()))
Definition float_cmp.cc:331
static I trunc(const T &val, typename EpsilonType< T >::Type epsilon=(DefaultEpsilon< T, cstyle >::value()))
Definition float_cmp.cc:341
static I trunc(const T &val, typename EpsilonType< T >::Type epsilon=(DefaultEpsilon< T, cstyle >::value()))
Definition float_cmp.cc:350
static std::vector< I > trunc(const std::vector< T > &val, typename EpsilonType< T >::Type epsilon=(DefaultEpsilon< T, cstyle >::value()))
Definition float_cmp.cc:359
static Dune::FieldVector< I, n > trunc(const Dune::FieldVector< T, n > &val, typename EpsilonType< T >::Type epsilon=(DefaultEpsilon< T, cstyle >::value()))
Definition float_cmp.cc:371
mapping from a value type and a compare style to a default epsilon
Definition float_cmp.hh:136
static EpsilonType< T >::Type value()
Returns the default epsilon for the given value type and compare style.
bool le(const ValueType &first, const ValueType &second) const
test if first lesser or equal second
Definition float_cmp.cc:460
FloatCmpOps(EpsilonType epsilon=DefaultEpsilon::value())
construct an operations object
Definition float_cmp.cc:405
bool eq(const ValueType &first, const ValueType &second) const
test for equality using epsilon
Definition float_cmp.cc:425
bool lt(const ValueType &first, const ValueType &second) const
test if first lesser than second
Definition float_cmp.cc:446
bool ge(const ValueType &first, const ValueType &second) const
test if first greater or equal second
Definition float_cmp.cc:453
FloatCmp::EpsilonType< T >::Type EpsilonType
Type of the epsilon.
Definition float_cmp.hh:302
bool ne(const ValueType &first, const ValueType &second) const
test for inequality using epsilon
Definition float_cmp.cc:432
T ValueType
Type of the values to compare.
Definition float_cmp.hh:297
EpsilonType epsilon() const
return the current epsilon
Definition float_cmp.cc:410
I round(const ValueType &val) const
round using epsilon
Definition float_cmp.cc:469
I trunc(const ValueType &val) const
truncate using epsilon
Definition float_cmp.cc:477
bool gt(const ValueType &first, const ValueType &second) const
test if first greater than second
Definition float_cmp.cc:439