1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
|
/****************************************************************************
* Core Library Version 1.7, August 2004
* Copyright (c) 1995-2004 Exact Computation Project
* All rights reserved.
*
* This file is part of CORE (http://cs.nyu.edu/exact/core/); you may
* redistribute it under the terms of the Q Public License version 1.0.
* See the file LICENSE.QPL distributed with CORE.
*
* Licensees holding a valid commercial license may use this file in
* accordance with the commercial license agreement provided with the
* software.
*
* This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
* WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
*
*
* File: Filter.h
* Synopsis:
* This is a simple filtered floating point number,
* represented by the main class, FilterFp.
* based on the Burnikel-Funke-Schirra (BFS) filter scheme.
* We do not use IEEE exception mechanism here.
* It is used by the Expr class.
*
* Written by
* Zilin Du <zilin@cs.nyu.edu>
* Chee Yap <yap@cs.nyu.edu>
*
* WWW URL: http://cs.nyu.edu/exact/
* Email: exact@cs.nyu.edu
*
* $URL: svn+ssh://scm.gforge.inria.fr/svn/cgal/branches/CGAL-3.2-branch/Core/include/CORE/Filter.h $
* $Id: Filter.h 29485 2006-03-14 11:52:49Z efif $
***************************************************************************/
#ifndef _CORE_FILTER_H_
#define _CORE_FILTER_H_
#include <CORE/Real.h>
#include <math.h>
#if defined (_MSC_VER) || defined (__MINGW32__) // add support for MinGW
#define finite(x) _finite(x)
#define ilogb(x) (int)_logb(x)
#endif
#if defined(sun) || defined(__sun)
#include <ieeefp.h>
#endif
CORE_BEGIN_NAMESPACE
const int POWTWO_26 = (1 << 26); ///< constant 2^26
/// \class filteredFp Filter.h
/// \brief filteredFp represents filtered floating point
/// numbers, based on BFS filter
class filteredFp {
double fpVal; // approximate double value for some "real value"
double maxAbs; // if (|fpVal| > maxAbs * ind * 2^{-53}) then
int ind; // sign of value is sign(fpVal). Else, don't know.
// REFERENCE: Burnikel, Funke, Schirra (BFS) filter
// Chee: in isOK(), you used the test "|fpVal| >= maxAbs * ind * 2^{-53}"
// which seems to be correct (i.e., not |fpVal| > maxAbs * ind * 2^{-53})
public:
/// \name Constructors
//@{
/// constructor
filteredFp (double val = 0.0)
: fpVal(val), maxAbs(core_abs(val)), ind(0) {}
/// constructor
filteredFp (double val, double m, int i)
: fpVal(val), maxAbs(m), ind(i) {}
/// construct a filteredFp from Real v.
/** if v causes an overflow, fpVal = +/- Infty
if v causes an underflow, fpVal = ...? */
filteredFp (const Real & value) : fpVal(0.0), maxAbs(0.0), ind(0) {
if (value != CORE_REAL_ZERO) {
ind = 1;
fpVal = value.doubleValue();
if (value.MSB() <= -1075)
maxAbs = 1;
else
maxAbs = core_abs(fpVal); // NaN are propagated correctly by core_abs.
}
}
//@}
/// \name Help Functions
//@{
/// return filtered value (for debug)
double getValue() const {
return fpVal;
}
/// check whether the sign (!) of the filtered value is OK
bool isOK() const {
return (fpFilterFlag && // To disable filter
finite(fpVal) && // Test for infinite and NaNs
(core_abs(fpVal) >= maxAbs*ind*CORE_EPS));
}
/// return the sign of fitered value.
/** (Note: must call isOK() to check whether the sign is ok
before call this function.) */
int sign() const {
#ifdef CORE_DEBUG
assert(isOK());
#endif
if (fpVal == 0.0)
return 0;
else
return fpVal > 0.0 ? 1: -1;
}
/// lower bound on MSB
/** defined to be cel(lg(real value));
ilogb(x) is floor(log_2(|x|)).
Also, ilogb(0) = -INT_MAX. ilogb(NaN) = ilogb(+/-Inf) = INT_MAX */
extLong lMSB() const {
return extLong(ilogb(core_abs(fpVal)-maxAbs*ind*CORE_EPS));
}
/// upper bound on MSB
extLong uMSB() const {
return extLong(ilogb(core_abs(fpVal)+maxAbs*ind*CORE_EPS));
}
//@}
/// \name Operators
//@{
/// unary minus
filteredFp operator -() const {
return filteredFp(-fpVal, maxAbs, ind);
}
/// addition
filteredFp operator+ (const filteredFp& x) const {
return filteredFp(fpVal+x.fpVal, maxAbs+x.maxAbs, 1+core_max(ind, x.ind));
}
/// subtraction
filteredFp operator- (const filteredFp& x) const {
return filteredFp(fpVal-x.fpVal, maxAbs+x.maxAbs, 1+core_max(ind, x.ind));
}
/// multiplication
filteredFp operator* (const filteredFp& x) const {
return filteredFp(fpVal*x.fpVal, maxAbs*x.maxAbs+DBL_MIN, 1+ind+x.ind);
}
/// division
filteredFp operator/ (const filteredFp& x) const {
if (x.fpVal == 0.0)
core_error("possible zero divisor!", __FILE__, __LINE__, false);
double xxx = core_abs(x.fpVal) / x.maxAbs - (x.ind+1)*CORE_EPS + DBL_MIN;
if (xxx > 0) {
double val = fpVal / x.fpVal;
double maxVal = ( core_abs(val) + maxAbs / x.maxAbs) / xxx + DBL_MIN;
return filteredFp(val, maxVal, 1 + core_max(ind, x.ind + 1));
} else
return filteredFp(getDoubleInfty(), 0.0, 0);
}
/// square root
filteredFp sqrt () const {
if (fpVal < 0.0)
core_error("possible negative sqrt!", __FILE__, __LINE__, false);
if (fpVal > 0.0) {
double val = std::sqrt(fpVal);
return filteredFp(val, ( maxAbs / fpVal ) * val, 1 + ind);
} else
return filteredFp(0.0, std::sqrt(maxAbs) * POWTWO_26, 1 + ind);
}
/// dump function
void dump (std::ostream&os) const {
os << "Filter=[fpVal=" << fpVal << ",maxAbs=" << maxAbs << ",ind=" << ind << "]";
}
/// helper function (to avoid warning under some compilers)
static double getDoubleInfty() {
static double d = DBL_MAX;
return 2*d;
}
//@}
}; //filteredFp class
inline std::ostream & operator<< (std::ostream & os, const filteredFp& fp) {
fp.dump(os);
return os;
}
CORE_END_NAMESPACE
#endif // _CORE_FILTER_H_
|
