dune-common 2.9.0
Loading...
Searching...
No Matches
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// SPDX-FileCopyrightInfo: Copyright (C) DUNE Project contributors, see file LICENSE.md in module root
4// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
5#include "float_cmp.hh"
6
7#include <vector>
8#include <limits>
9#include <algorithm>
10#include <cstdlib>
12
13namespace Dune {
14
15
16 namespace FloatCmp {
17 // traits
19
23 template<class T> struct EpsilonType {
25 typedef T Type;
26 };
28
33 template<class T, typename A>
34 struct EpsilonType<std::vector<T, A> > {
36 typedef typename EpsilonType<T>::Type Type;
37 };
39
44 template<class T, int n>
45 struct EpsilonType<FieldVector<T, n> > {
47 typedef typename EpsilonType<T>::Type Type;
48 };
49
50 // default epsilon
51 template<class T>
53 static typename EpsilonType<T>::Type value()
54 { return std::numeric_limits<typename EpsilonType<T>::Type>::epsilon()*8.; }
55 };
56 template<class T>
58 static typename EpsilonType<T>::Type value()
59 { return std::numeric_limits<typename EpsilonType<T>::Type>::epsilon()*8.; }
60 };
61 template<class T>
63 static typename EpsilonType<T>::Type value()
64 { return std::max(std::numeric_limits<typename EpsilonType<T>::Type>::epsilon(), 1e-6); }
65 };
66
67 namespace Impl {
68 // basic comparison
69 template<class T, CmpStyle style = defaultCmpStyle>
70 struct eq_t;
71 template<class T>
72 struct eq_t<T, relativeWeak> {
73 static bool eq(const T &first,
74 const T &second,
76 {
77 using std::abs;
78 return abs(first - second) <= epsilon*std::max(abs(first), abs(second));
79 }
80 };
81 template<class T>
82 struct eq_t<T, relativeStrong> {
83 static bool eq(const T &first,
84 const T &second,
86 {
87 using std::abs;
88 return abs(first - second) <= epsilon*std::min(abs(first), abs(second));
89 }
90 };
91 template<class T>
92 struct eq_t<T, absolute> {
93 static bool eq(const T &first,
94 const T &second,
96 {
97 using std::abs;
98 return abs(first-second) <= epsilon;
99 }
100 };
101 template<class T, CmpStyle cstyle>
102 struct eq_t_std_vec {
103 typedef std::vector<T> V;
104 static bool eq(const V &first,
105 const V &second,
106 typename EpsilonType<V>::Type epsilon = DefaultEpsilon<V>::value()) {
107 auto size = first.size();
108 if(size != second.size()) return false;
109 for(unsigned int i = 0; i < size; ++i)
110 if(!eq_t<T, cstyle>::eq(first[i], second[i], epsilon))
111 return false;
112 return true;
113 }
114 };
115 template< class T>
116 struct eq_t<std::vector<T>, relativeWeak> : eq_t_std_vec<T, relativeWeak> {};
117 template< class T>
118 struct eq_t<std::vector<T>, relativeStrong> : eq_t_std_vec<T, relativeStrong> {};
119 template< class T>
120 struct eq_t<std::vector<T>, absolute> : eq_t_std_vec<T, absolute> {};
121
122 template<class T, int n, CmpStyle cstyle>
123 struct eq_t_fvec {
124 typedef Dune::FieldVector<T, n> V;
125 static bool eq(const V &first,
126 const V &second,
127 typename EpsilonType<V>::Type epsilon = DefaultEpsilon<V>::value()) {
128 for(int i = 0; i < n; ++i)
129 if(!eq_t<T, cstyle>::eq(first[i], second[i], epsilon))
130 return false;
131 return true;
132 }
133 };
134 template< class T, int n >
135 struct eq_t< Dune::FieldVector<T, n>, relativeWeak> : eq_t_fvec<T, n, relativeWeak> {};
136 template< class T, int n >
137 struct eq_t< Dune::FieldVector<T, n>, relativeStrong> : eq_t_fvec<T, n, relativeStrong> {};
138 template< class T, int n >
139 struct eq_t< Dune::FieldVector<T, n>, absolute> : eq_t_fvec<T, n, absolute> {};
140 } // namespace Impl
141
142 // operations in functional style
143 template <class T, CmpStyle style>
144 bool eq(const T &first,
145 const T &second,
146 typename EpsilonType<T>::Type epsilon)
147 {
148 return Impl::eq_t<T, style>::eq(first, second, epsilon);
149 }
150 template <class T, CmpStyle style>
151 bool ne(const T &first,
152 const T &second,
153 typename EpsilonType<T>::Type epsilon)
154 {
155 return !eq<T, style>(first, second, epsilon);
156 }
157 template <class T, CmpStyle style>
158 bool gt(const T &first,
159 const T &second,
160 typename EpsilonType<T>::Type epsilon)
161 {
162 return first > second && ne<T, style>(first, second, epsilon);
163 }
164 template <class T, CmpStyle style>
165 bool lt(const T &first,
166 const T &second,
167 typename EpsilonType<T>::Type epsilon)
168 {
169 return first < second && ne<T, style>(first, second, epsilon);
170 }
171 template <class T, CmpStyle style>
172 bool ge(const T &first,
173 const T &second,
174 typename EpsilonType<T>::Type epsilon)
175 {
176 return first > second || eq<T, style>(first, second, epsilon);
177 }
178 template <class T, CmpStyle style>
179 bool le(const T &first,
180 const T &second,
181 typename EpsilonType<T>::Type epsilon)
182 {
183 return first < second || eq<T, style>(first, second, epsilon);
184 }
185
186 // default template arguments
187 template <class T>
188 bool eq(const T &first,
189 const T &second,
191 {
192 return eq<T, defaultCmpStyle>(first, second, epsilon);
193 }
194 template <class T>
195 bool ne(const T &first,
196 const T &second,
198 {
199 return ne<T, defaultCmpStyle>(first, second, epsilon);
200 }
201 template <class T>
202 bool gt(const T &first,
203 const T &second,
205 {
206 return gt<T, defaultCmpStyle>(first, second, epsilon);
207 }
208 template <class T>
209 bool lt(const T &first,
210 const T &second,
212 {
213 return lt<T, defaultCmpStyle>(first, second, epsilon);
214 }
215 template <class T>
216 bool ge(const T &first,
217 const T &second,
219 {
220 return ge<T, defaultCmpStyle>(first, second, epsilon);
221 }
222 template <class T>
223 bool le(const T &first,
224 const T &second,
226 {
227 return le<T, defaultCmpStyle>(first, second, epsilon);
228 }
229
230 // rounding operations
231 namespace Impl {
232 template<class I, class T, CmpStyle cstyle = defaultCmpStyle, RoundingStyle rstyle = defaultRoundingStyle>
233 struct round_t;
234 template<class I, class T, CmpStyle cstyle>
235 struct round_t<I, T, cstyle, downward> {
236 static I
237 round(const T &val,
239 // first get an approximation
240 I lower = I(val);
241 I upper;
242 if(eq<T, cstyle>(T(lower), val, epsilon)) return lower;
243 if(T(lower) > val) { upper = lower; lower--; }
244 else upper = lower+1;
245 if(le<T, cstyle>(val - T(lower), T(upper) - val, epsilon))
246 return lower;
247 else return upper;
248 }
249 };
250 template<class I, class T, CmpStyle cstyle>
251 struct round_t<I, T, cstyle, upward> {
252 static I
253 round(const T &val,
255 // first get an approximation
256 I lower = I(val);
257 I upper;
258 if(eq<T, cstyle>(T(lower), val, epsilon)) return lower;
259 if(T(lower) > val) { upper = lower; lower--; }
260 else upper = lower+1;
261 if(lt<T, cstyle>(val - T(lower), T(upper) - val, epsilon))
262 return lower;
263 else return upper;
264 }
265 };
266 template<class I, class T, CmpStyle cstyle>
267 struct round_t<I, T, cstyle, towardZero> {
268 static I
269 round(const T &val,
271 if(val > T(0))
272 return round_t<I, T, cstyle, downward>::round(val, epsilon);
273 else return round_t<I, T, cstyle, upward>::round(val, epsilon);
274 }
275 };
276 template<class I, class T, CmpStyle cstyle>
277 struct round_t<I, T, cstyle, towardInf> {
278 static I
279 round(const T &val,
281 if(val > T(0))
282 return round_t<I, T, cstyle, upward>::round(val, epsilon);
283 else return round_t<I, T, cstyle, downward>::round(val, epsilon);
284 }
285 };
286 template<class I, class T, CmpStyle cstyle, RoundingStyle rstyle>
287 struct round_t<std::vector<I>, std::vector<T>, cstyle, rstyle> {
288 static std::vector<I>
289 round(const T &val,
291 unsigned int size = val.size();
292 std::vector<I> res(size);
293 for(unsigned int i = 0; i < size; ++i)
294 res[i] = round_t<I, T, cstyle, rstyle>::round(val[i], epsilon);
295 return res;
296 }
297 };
298 template<class I, class T, int n, CmpStyle cstyle, RoundingStyle rstyle>
299 struct round_t<Dune::FieldVector<I, n>, Dune::FieldVector<T, n>, cstyle, rstyle> {
301 round(const T &val,
304 for(int i = 0; i < n; ++i)
305 res[i] = round_t<I, T, cstyle, rstyle>::round(val[i], epsilon);
306 return res;
307 }
308 };
309 } // end namespace Impl
310 template<class I, class T, CmpStyle cstyle, RoundingStyle rstyle>
311 I round(const T &val, typename EpsilonType<T>::Type epsilon /*= DefaultEpsilon<T, cstyle>::value()*/)
312 {
313 return Impl::round_t<I, T, cstyle, rstyle>::round(val, epsilon);
314 }
315 template<class I, class T, CmpStyle cstyle>
316 I round(const T &val, typename EpsilonType<T>::Type epsilon = DefaultEpsilon<T, cstyle>::value())
317 {
318 return round<I, T, cstyle, defaultRoundingStyle>(val, epsilon);
319 }
320 template<class I, class T, RoundingStyle rstyle>
321 I round(const T &val, typename EpsilonType<T>::Type epsilon = DefaultEpsilon<T, defaultCmpStyle>::value())
322 {
323 return round<I, T, defaultCmpStyle, rstyle>(val, epsilon);
324 }
325 template<class I, class T>
327 {
328 return round<I, T, defaultCmpStyle>(val, epsilon);
329 }
330
331 // truncation
332 namespace Impl {
333 template<class I, class T, CmpStyle cstyle = defaultCmpStyle, RoundingStyle rstyle = defaultRoundingStyle>
334 struct trunc_t;
335 template<class I, class T, CmpStyle cstyle>
336 struct trunc_t<I, T, cstyle, downward> {
337 static I
338 trunc(const T &val,
340 // this should be optimized away unless needed
341 if(!std::numeric_limits<I>::is_signed)
342 // make sure this works for all useful cases even if I is an unsigned type
343 if(eq<T, cstyle>(val, T(0), epsilon)) return I(0);
344 // first get an approximation
345 I lower = I(val); // now |val-lower| < 1
346 // make sure we're really lower in case the cast truncated to an unexpected direction
347 if(T(lower) > val) lower--; // now val-lower < 1
348 // check whether lower + 1 is approximately val
349 if(eq<T, cstyle>(T(lower+1), val, epsilon))
350 return lower+1;
351 else return lower;
352 }
353 };
354 template<class I, class T, CmpStyle cstyle>
355 struct trunc_t<I, T, cstyle, upward> {
356 static I
357 trunc(const T &val,
359 I upper = trunc_t<I, T, cstyle, downward>::trunc(val, epsilon);
360 if(ne<T, cstyle>(T(upper), val, epsilon)) ++upper;
361 return upper;
362 }
363 };
364 template<class I, class T, CmpStyle cstyle>
365 struct trunc_t<I, T, cstyle, towardZero> {
366 static I
367 trunc(const T &val,
369 if(val > T(0)) return trunc_t<I, T, cstyle, downward>::trunc(val, epsilon);
370 else return trunc_t<I, T, cstyle, upward>::trunc(val, epsilon);
371 }
372 };
373 template<class I, class T, CmpStyle cstyle>
374 struct trunc_t<I, T, cstyle, towardInf> {
375 static I
376 trunc(const T &val,
378 if(val > T(0)) return trunc_t<I, T, cstyle, upward>::trunc(val, epsilon);
379 else return trunc_t<I, T, cstyle, downward>::trunc(val, epsilon);
380 }
381 };
382 template<class I, class T, CmpStyle cstyle, RoundingStyle rstyle>
383 struct trunc_t<std::vector<I>, std::vector<T>, cstyle, rstyle> {
384 static std::vector<I>
385 trunc(const std::vector<T> &val,
387 unsigned int size = val.size();
388 std::vector<I> res(size);
389 for(unsigned int i = 0; i < size; ++i)
390 res[i] = trunc_t<I, T, cstyle, rstyle>::trunc(val[i], epsilon);
391 return res;
392 }
393 };
394 template<class I, class T, int n, CmpStyle cstyle, RoundingStyle rstyle>
395 struct trunc_t<Dune::FieldVector<I, n>, Dune::FieldVector<T, n>, cstyle, rstyle> {
400 for(int i = 0; i < n; ++i)
401 res[i] = trunc_t<I, T, cstyle, rstyle>::trunc(val[i], epsilon);
402 return res;
403 }
404 };
405 } // namespace Impl
406 template<class I, class T, CmpStyle cstyle, RoundingStyle rstyle>
407 I trunc(const T &val, typename EpsilonType<T>::Type epsilon /*= DefaultEpsilon<T, cstyle>::value()*/)
408 {
409 return Impl::trunc_t<I, T, cstyle, rstyle>::trunc(val, epsilon);
410 }
411 template<class I, class T, CmpStyle cstyle>
412 I trunc(const T &val, typename EpsilonType<T>::Type epsilon = DefaultEpsilon<T, cstyle>::value())
413 {
414 return trunc<I, T, cstyle, defaultRoundingStyle>(val, epsilon);
415 }
416 template<class I, class T, RoundingStyle rstyle>
417 I trunc(const T &val, typename EpsilonType<T>::Type epsilon = DefaultEpsilon<T, defaultCmpStyle>::value())
418 {
419 return trunc<I, T, defaultCmpStyle, rstyle>(val, epsilon);
420 }
421 template<class I, class T>
423 {
424 return trunc<I, T, defaultCmpStyle>(val, epsilon);
425 }
426 } //namespace Dune
427
428 // oo interface
429 template<class T, FloatCmp::CmpStyle cstyle_, FloatCmp::RoundingStyle rstyle_>
432
433
434 template<class T, FloatCmp::CmpStyle cstyle_, FloatCmp::RoundingStyle rstyle_>
437 {
438 return epsilon_;
439 }
440
441 template<class T, FloatCmp::CmpStyle cstyle_, FloatCmp::RoundingStyle rstyle_>
442 void
444 {
445 epsilon_ = epsilon__;
446 }
447
448
449 template<class T, FloatCmp::CmpStyle cstyle_, FloatCmp::RoundingStyle rstyle_>
451 eq(const ValueType &first, const ValueType &second) const
452 {
453 return Dune::FloatCmp::eq<ValueType, cstyle>(first, second, epsilon_);
454 }
455
456 template<class T, FloatCmp::CmpStyle cstyle_, FloatCmp::RoundingStyle rstyle_>
458 ne(const ValueType &first, const ValueType &second) const
459 {
460 return Dune::FloatCmp::ne<ValueType, cstyle>(first, second, epsilon_);
461 }
462
463 template<class T, FloatCmp::CmpStyle cstyle_, FloatCmp::RoundingStyle rstyle_>
465 gt(const ValueType &first, const ValueType &second) const
466 {
467 return Dune::FloatCmp::gt<ValueType, cstyle>(first, second, epsilon_);
468 }
469
470 template<class T, FloatCmp::CmpStyle cstyle_, FloatCmp::RoundingStyle rstyle_>
472 lt(const ValueType &first, const ValueType &second) const
473 {
474 return Dune::FloatCmp::lt<ValueType, cstyle>(first, second, epsilon_);
475 }
476
477 template<class T, FloatCmp::CmpStyle cstyle_, FloatCmp::RoundingStyle rstyle_>
479 ge(const ValueType &first, const ValueType &second) const
480 {
481 return Dune::FloatCmp::ge<ValueType, cstyle>(first, second, epsilon_);
482 }
483
484 template<class T, FloatCmp::CmpStyle cstyle_, FloatCmp::RoundingStyle rstyle_>
486 le(const ValueType &first, const ValueType &second) const
487 {
488 return Dune::FloatCmp::le<ValueType, cstyle>(first, second, epsilon_);
489 }
490
491
492 template<class T, FloatCmp::CmpStyle cstyle_, FloatCmp::RoundingStyle rstyle_>
493 template<class I>
495 round(const ValueType &val) const
496 {
497 return Dune::FloatCmp::round<I, ValueType, cstyle, rstyle_>(val, epsilon_);
498 }
499
500 template<class T, FloatCmp::CmpStyle cstyle_, FloatCmp::RoundingStyle rstyle_>
501 template<class I>
503 trunc(const ValueType &val) const
504 {
505 return Dune::FloatCmp::trunc<I, ValueType, cstyle, rstyle_>(val, epsilon_);
506 }
507
508} //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:151
bool eq(const T &first, const T &second, typename EpsilonType< T >::Type epsilon)
test for equality using epsilon
Definition float_cmp.cc:144
I round(const T &val, typename EpsilonType< T >::Type epsilon)
round using epsilon
Definition float_cmp.cc:311
I trunc(const T &val, typename EpsilonType< T >::Type epsilon)
truncate using epsilon
Definition float_cmp.cc:407
bool lt(const T &first, const T &second, typename EpsilonType< T >::Type epsilon)
test if first lesser than second
Definition float_cmp.cc:165
bool gt(const T &first, const T &second, typename EpsilonType< T >::Type epsilon)
test if first greater than second
Definition float_cmp.cc:158
bool ge(const T &first, const T &second, typename EpsilonType< T >::Type epsilon)
test if first greater or equal second
Definition float_cmp.cc:172
bool le(const T &first, const T &second, typename EpsilonType< T >::Type epsilon)
test if first lesser or equal second
Definition float_cmp.cc:179
@ relativeStrong
|a-b|/|a| <= epsilon && |a-b|/|b| <= epsilon
Definition float_cmp.hh:108
@ relativeWeak
|a-b|/|a| <= epsilon || |a-b|/|b| <= epsilon
Definition float_cmp.hh:106
@ absolute
|a-b| <= epsilon
Definition float_cmp.hh:110
@ towardZero
always round toward 0
Definition float_cmp.hh:118
@ towardInf
always round away from 0
Definition float_cmp.hh:120
@ upward
round toward
Definition float_cmp.hh:124
@ downward
round toward
Definition float_cmp.hh:122
STL namespace.
Dune namespace.
Definition alignedallocator.hh:13
vector space out of a tensor product of fields.
Definition fvector.hh:95
Mapping of value type to epsilon type.
Definition float_cmp.cc:23
T Type
The epsilon type corresponding to value type T.
Definition float_cmp.cc:25
EpsilonType< T >::Type Type
The epsilon type corresponding to value type std::vector<T, A>
Definition float_cmp.cc:36
EpsilonType< T >::Type Type
The epsilon type corresponding to value type Dune::FieldVector<T, n>
Definition float_cmp.cc:47
static EpsilonType< T >::Type value()
Definition float_cmp.cc:53
static EpsilonType< T >::Type value()
Definition float_cmp.cc:58
static EpsilonType< T >::Type value()
Definition float_cmp.cc:63
mapping from a value type and a compare style to a default epsilon
Definition float_cmp.hh:138
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:486
FloatCmpOps(EpsilonType epsilon=DefaultEpsilon::value())
construct an operations object
Definition float_cmp.cc:431
bool eq(const ValueType &first, const ValueType &second) const
test for equality using epsilon
Definition float_cmp.cc:451
bool lt(const ValueType &first, const ValueType &second) const
test if first lesser than second
Definition float_cmp.cc:472
bool ge(const ValueType &first, const ValueType &second) const
test if first greater or equal second
Definition float_cmp.cc:479
FloatCmp::EpsilonType< T >::Type EpsilonType
Type of the epsilon.
Definition float_cmp.hh:304
bool ne(const ValueType &first, const ValueType &second) const
test for inequality using epsilon
Definition float_cmp.cc:458
T ValueType
Type of the values to compare.
Definition float_cmp.hh:299
EpsilonType epsilon() const
return the current epsilon
Definition float_cmp.cc:436
I round(const ValueType &val) const
round using epsilon
Definition float_cmp.cc:495
I trunc(const ValueType &val) const
truncate using epsilon
Definition float_cmp.cc:503
bool gt(const ValueType &first, const ValueType &second) const
test if first greater than second
Definition float_cmp.cc:465