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
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
|
// Boost.Geometry (aka GGL, Generic Geometry Library)
// Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands.
// Copyright (c) 2008-2012 Bruno Lalande, Paris, France.
// Copyright (c) 2009-2012 Mateusz Loskot, London, UK.
// Parts of Boost.Geometry are redesigned from Geodan's Geographic Library
// (geolib/GGL), copyright (c) 1995-2010 Geodan, Amsterdam, the Netherlands.
// Use, modification and distribution is subject to the Boost Software License,
// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)
#ifndef BOOST_GEOMETRY_UTIL_MATH_HPP
#define BOOST_GEOMETRY_UTIL_MATH_HPP
#include <cmath>
#include <limits>
#include <boost/math/constants/constants.hpp>
#include <boost/geometry/util/select_most_precise.hpp>
namespace boost { namespace geometry
{
namespace math
{
#ifndef DOXYGEN_NO_DETAIL
namespace detail
{
template <typename Type, bool IsFloatingPoint>
struct equals
{
static inline bool apply(Type const& a, Type const& b)
{
return a == b;
}
};
template <typename Type>
struct equals<Type, true>
{
static inline Type get_max(Type const& a, Type const& b, Type const& c)
{
return (std::max)((std::max)(a, b), c);
}
static inline bool apply(Type const& a, Type const& b)
{
if (a == b)
{
return true;
}
// See http://www.parashift.com/c++-faq-lite/newbie.html#faq-29.17,
// FUTURE: replace by some boost tool or boost::test::close_at_tolerance
return std::abs(a - b) <= std::numeric_limits<Type>::epsilon() * get_max(std::abs(a), std::abs(b), 1.0);
}
};
template <typename Type, bool IsFloatingPoint>
struct smaller
{
static inline bool apply(Type const& a, Type const& b)
{
return a < b;
}
};
template <typename Type>
struct smaller<Type, true>
{
static inline bool apply(Type const& a, Type const& b)
{
if (equals<Type, true>::apply(a, b))
{
return false;
}
return a < b;
}
};
template <typename Type, bool IsFloatingPoint>
struct equals_with_epsilon : public equals<Type, IsFloatingPoint> {};
/*!
\brief Short construct to enable partial specialization for PI, currently not possible in Math.
*/
template <typename T>
struct define_pi
{
static inline T apply()
{
// Default calls Boost.Math
return boost::math::constants::pi<T>();
}
};
} // namespace detail
#endif
template <typename T>
inline T pi() { return detail::define_pi<T>::apply(); }
// Maybe replace this by boost equals or boost ublas numeric equals or so
/*!
\brief returns true if both arguments are equal.
\ingroup utility
\param a first argument
\param b second argument
\return true if a == b
\note If both a and b are of an integral type, comparison is done by ==.
If one of the types is floating point, comparison is done by abs and
comparing with epsilon. If one of the types is non-fundamental, it might
be a high-precision number and comparison is done using the == operator
of that class.
*/
template <typename T1, typename T2>
inline bool equals(T1 const& a, T2 const& b)
{
typedef typename select_most_precise<T1, T2>::type select_type;
return detail::equals
<
select_type,
boost::is_floating_point<select_type>::type::value
>::apply(a, b);
}
template <typename T1, typename T2>
inline bool equals_with_epsilon(T1 const& a, T2 const& b)
{
typedef typename select_most_precise<T1, T2>::type select_type;
return detail::equals_with_epsilon
<
select_type,
boost::is_floating_point<select_type>::type::value
>::apply(a, b);
}
template <typename T1, typename T2>
inline bool smaller(T1 const& a, T2 const& b)
{
typedef typename select_most_precise<T1, T2>::type select_type;
return detail::smaller
<
select_type,
boost::is_floating_point<select_type>::type::value
>::apply(a, b);
}
template <typename T1, typename T2>
inline bool larger(T1 const& a, T2 const& b)
{
typedef typename select_most_precise<T1, T2>::type select_type;
return detail::smaller
<
select_type,
boost::is_floating_point<select_type>::type::value
>::apply(b, a);
}
double const d2r = geometry::math::pi<double>() / 180.0;
double const r2d = 1.0 / d2r;
/*!
\brief Calculates the haversine of an angle
\ingroup utility
\note See http://en.wikipedia.org/wiki/Haversine_formula
haversin(alpha) = sin2(alpha/2)
*/
template <typename T>
inline T hav(T const& theta)
{
T const half = T(0.5);
T const sn = sin(half * theta);
return sn * sn;
}
/*!
\brief Short utility to return the square
\ingroup utility
\param value Value to calculate the square from
\return The squared value
*/
template <typename T>
inline T sqr(T const& value)
{
return value * value;
}
/*!
\brief Short utility to workaround gcc/clang problem that abs is converting to integer
\ingroup utility
*/
template<typename T>
inline T abs(const T& t)
{
using std::abs;
return abs(t);
}
} // namespace math
}} // namespace boost::geometry
#endif // BOOST_GEOMETRY_UTIL_MATH_HPP
|
