1 // Copyright (c) 1999,2000 INRIA Sophia-Antipolis (France).
2 // All rights reserved.
3 //
4 // This file is part of CGAL (www.cgal.org)
5 //
6 // $URL: https://github.com/CGAL/cgal/blob/v5.3/Filtered_kernel/include/CGAL/internal/Static_filters/Static_filter_error.h $
7 // $Id: Static_filter_error.h 0779373 2020-03-26T13:31:46+01:00 Sébastien Loriot
8 // SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial
9 //
10 //
11 // Author(s) : Sylvain Pion
12
13 #ifndef CGAL_INTERNAL_STATIC_FILTER_ERROR_H
14 #define CGAL_INTERNAL_STATIC_FILTER_ERROR_H
15
16 // This file contains the description of the class Static_filter_error.
17 // The goal of this class is to be run by some overloaded predicates,
18 // to compute error bound done in these functions.
19 //
20 // The original idea is from Olivier Devillers.
21
22 // It is still EXPERIMENTAL.
23
24 // TODO:
25 // - I need to add some missing operators and functions, min/max...
26 // - Remove the degree stuff, it's only meant for debug (?).
27 // - Add __attribute__((const)) for optimizing ?
28
29 #include <CGAL/config.h>
30 #include <CGAL/FPU.h>
31
32 #include <boost/math/special_functions/next.hpp>
33
34 namespace CGAL { namespace internal {
35
36 struct Static_filter_error
37 {
38 typedef Static_filter_error Sfe;
39
Static_filter_errorStatic_filter_error40 Static_filter_error () {}
41
42 Static_filter_error (const int &i, const double &e = 0, const int &d = 1)
_bStatic_filter_error43 : _b(i), _e(e), _d(d) {}
44
45 Static_filter_error (const double &b, const double &e = 0, const int &d = 1)
_bStatic_filter_error46 : _b(b), _e(e), _d(d) {}
47
ulpStatic_filter_error48 static double ulp ()
49 {
50 FPU_CW_t backup = FPU_get_and_set_cw(CGAL_FE_UPWARD);
51 double e = ulp(1);
52 FPU_set_cw(backup);
53 return e;
54 }
55
ulpStatic_filter_error56 static double ulp (double d)
57 {
58 // You are supposed to call this function with rounding towards
59 // +infinity, and on a positive number.
60 d = CGAL_IA_FORCE_TO_DOUBLE(d); // stop constant propagation.
61 CGAL_assertion(d>=0);
62 double u;
63 double nd = boost::math::float_next(d);
64 u = nd - d;
65 // Then add extra bonus, because of Intel's extended precision feature.
66 // (ulp can be 2^-53 + 2^-64)
67 u += u / (1<<11);
68 CGAL_assertion(u!=0);
69 return u;
70 }
71
72 Sfe operator+ (const Sfe &f) const
73 {
74 CGAL_warning_msg( _d == f._d ,
75 "you are adding variables of different homogeneous degree");
76 // We have to add an ulp, since the homogeneization could induce such
77 // an error.
78 FPU_CW_t backup = FPU_get_and_set_cw(CGAL_FE_UPWARD);
79 double b = _b + f._b;
80 double u = ulp(b) / 2;
81 b += u;
82 double e = u + _e + f._e;
83 FPU_set_cw(backup);
84 return Sfe(b, e, _d);
85 }
86
87 Sfe operator* (const Sfe &f) const
88 {
89 // We have to add an ulp, since the homogeneization could induce such
90 // an error.
91 FPU_CW_t backup = FPU_get_and_set_cw(CGAL_FE_UPWARD);
92 double b = _b * f._b;
93 double u = ulp(b) / 2;
94 b += u;
95 double e = u + _e * f._e + _e * f._b + _b * f._e;
96 FPU_set_cw(backup);
97 return Sfe(b, e, _d+f._d);
98 }
99
100 Sfe operator- (const Sfe &f) const { return *this + f; }
101 Sfe operator- () const { return *this; }
102 // Sfe operator/ (const Sfe &) const { CGAL_error(); }
103 // Division not supported.
104
105 Sfe& operator+=(const Sfe &f) { return *this = *this + f; }
106 Sfe& operator-=(const Sfe &f) { return *this = *this - f; }
107 Sfe& operator*=(const Sfe &f) { return *this = *this * f; }
108 // Sfe& operator/=(const Sfe &f) { return *this = *this / f; }
109
errorStatic_filter_error110 double error() const { return _e; }
boundStatic_filter_error111 double bound() const { return _b; }
degreeStatic_filter_error112 int degree() const { return _d; }
113
114 bool operator< (const Sfe &f) const
115 {
116 Sfe e = *this + f;
117 std::cerr << "Static error is : " << e.error() << std::endl;
118 CGAL_error();
119 return false;
120 }
121 bool operator> (const Sfe &f) const { return *this < f; }
122 bool operator<=(const Sfe &f) const { return *this < f; }
123 bool operator>=(const Sfe &f) const { return *this < f; }
124 bool operator==(const Sfe &f) const { return *this < f; }
125 bool operator!=(const Sfe &f) const { return *this < f; }
126
127 private:
128 // _b is a bound on the absolute value of the _double_ value of the
129 // variable.
130 // _e is a bound on the absolute error (difference between _b and the
131 // _real_ value of the variable.
132 // _d is the degree of the variable, it allows some additionnal checks.
133 double _b, _e;
134 int _d;
135 };
136
137 inline
138 Static_filter_error
sqrt(const Static_filter_error & f)139 sqrt(const Static_filter_error &f)
140 {
141 CGAL_warning_msg( (f.degree() & 1) == 0,
142 "Do you really want a non integral degree ???");
143 // We have to add an ulp, since the homogeneization could induce such
144 // an error.
145 FPU_CW_t backup = FPU_get_and_set_cw(CGAL_FE_UPWARD);
146 double b = std::sqrt(f.bound());
147 double u = Static_filter_error::ulp(b) / 2;
148 b += u;
149 double e = std::sqrt(f.error()) + u;
150 FPU_set_cw(backup);
151 return Static_filter_error(b, e, f.degree()/2);
152 }
153
154 } } // namespace CGAL::internal
155
156 #endif // CGAL_INTERNAL_STATIC_FILTER_ERROR_H
157