]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | // Boost.Geometry (aka GGL, Generic Geometry Library) |
2 | ||
3 | // Copyright (c) 2007-2015 Barend Gehrels, Amsterdam, the Netherlands. | |
4 | // Copyright (c) 2008-2015 Bruno Lalande, Paris, France. | |
5 | // Copyright (c) 2009-2015 Mateusz Loskot, London, UK. | |
6 | ||
92f5a8d4 TL |
7 | // This file was modified by Oracle on 2014, 2015, 2018, 2019. |
8 | // Modifications copyright (c) 2014-2019, Oracle and/or its affiliates. | |
7c673cae FG |
9 | |
10 | // Contributed and/or modified by Menelaos Karavelas, on behalf of Oracle | |
11 | // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle | |
92f5a8d4 | 12 | // Contributed and/or modified by Adeel Ahmad, as part of Google Summer of Code 2018 program |
7c673cae FG |
13 | |
14 | // Parts of Boost.Geometry are redesigned from Geodan's Geographic Library | |
15 | // (geolib/GGL), copyright (c) 1995-2010 Geodan, Amsterdam, the Netherlands. | |
16 | ||
17 | // Use, modification and distribution is subject to the Boost Software License, | |
18 | // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at | |
19 | // http://www.boost.org/LICENSE_1_0.txt) | |
20 | ||
21 | #ifndef BOOST_GEOMETRY_UTIL_MATH_HPP | |
22 | #define BOOST_GEOMETRY_UTIL_MATH_HPP | |
23 | ||
24 | #include <cmath> | |
25 | #include <limits> | |
26 | ||
27 | #include <boost/core/ignore_unused.hpp> | |
28 | ||
29 | #include <boost/math/constants/constants.hpp> | |
30 | #include <boost/math/special_functions/fpclassify.hpp> | |
31 | //#include <boost/math/special_functions/round.hpp> | |
32 | #include <boost/numeric/conversion/cast.hpp> | |
33 | #include <boost/type_traits/is_fundamental.hpp> | |
34 | #include <boost/type_traits/is_integral.hpp> | |
35 | ||
36 | #include <boost/geometry/core/cs.hpp> | |
37 | ||
38 | #include <boost/geometry/util/select_most_precise.hpp> | |
39 | ||
40 | namespace boost { namespace geometry | |
41 | { | |
42 | ||
43 | namespace math | |
44 | { | |
45 | ||
46 | #ifndef DOXYGEN_NO_DETAIL | |
47 | namespace detail | |
48 | { | |
49 | ||
50 | template <typename T> | |
51 | inline T const& greatest(T const& v1, T const& v2) | |
52 | { | |
53 | return (std::max)(v1, v2); | |
54 | } | |
55 | ||
56 | template <typename T> | |
57 | inline T const& greatest(T const& v1, T const& v2, T const& v3) | |
58 | { | |
59 | return (std::max)(greatest(v1, v2), v3); | |
60 | } | |
61 | ||
62 | template <typename T> | |
63 | inline T const& greatest(T const& v1, T const& v2, T const& v3, T const& v4) | |
64 | { | |
65 | return (std::max)(greatest(v1, v2, v3), v4); | |
66 | } | |
67 | ||
68 | template <typename T> | |
69 | inline T const& greatest(T const& v1, T const& v2, T const& v3, T const& v4, T const& v5) | |
70 | { | |
71 | return (std::max)(greatest(v1, v2, v3, v4), v5); | |
72 | } | |
73 | ||
74 | ||
b32b8144 FG |
75 | template <typename T> |
76 | inline T bounded(T const& v, T const& lower, T const& upper) | |
77 | { | |
78 | return (std::min)((std::max)(v, lower), upper); | |
79 | } | |
80 | ||
81 | template <typename T> | |
82 | inline T bounded(T const& v, T const& lower) | |
83 | { | |
84 | return (std::max)(v, lower); | |
85 | } | |
86 | ||
87 | ||
7c673cae FG |
88 | template <typename T, |
89 | bool IsFloatingPoint = boost::is_floating_point<T>::value> | |
90 | struct abs | |
91 | { | |
92 | static inline T apply(T const& value) | |
93 | { | |
94 | T const zero = T(); | |
95 | return value < zero ? -value : value; | |
96 | } | |
97 | }; | |
98 | ||
99 | template <typename T> | |
100 | struct abs<T, true> | |
101 | { | |
102 | static inline T apply(T const& value) | |
103 | { | |
104 | using ::fabs; | |
105 | using std::fabs; // for long double | |
106 | ||
107 | return fabs(value); | |
108 | } | |
109 | }; | |
110 | ||
111 | ||
112 | struct equals_default_policy | |
113 | { | |
114 | template <typename T> | |
115 | static inline T apply(T const& a, T const& b) | |
116 | { | |
117 | // See http://www.parashift.com/c++-faq-lite/newbie.html#faq-29.17 | |
118 | return greatest(abs<T>::apply(a), abs<T>::apply(b), T(1)); | |
119 | } | |
120 | }; | |
121 | ||
122 | template <typename T, | |
123 | bool IsFloatingPoint = boost::is_floating_point<T>::value> | |
124 | struct equals_factor_policy | |
125 | { | |
126 | equals_factor_policy() | |
127 | : factor(1) {} | |
128 | explicit equals_factor_policy(T const& v) | |
129 | : factor(greatest(abs<T>::apply(v), T(1))) | |
130 | {} | |
131 | equals_factor_policy(T const& v0, T const& v1, T const& v2, T const& v3) | |
132 | : factor(greatest(abs<T>::apply(v0), abs<T>::apply(v1), | |
133 | abs<T>::apply(v2), abs<T>::apply(v3), | |
134 | T(1))) | |
135 | {} | |
136 | ||
137 | T const& apply(T const&, T const&) const | |
138 | { | |
139 | return factor; | |
140 | } | |
141 | ||
142 | T factor; | |
143 | }; | |
144 | ||
145 | template <typename T> | |
146 | struct equals_factor_policy<T, false> | |
147 | { | |
148 | equals_factor_policy() {} | |
149 | explicit equals_factor_policy(T const&) {} | |
150 | equals_factor_policy(T const& , T const& , T const& , T const& ) {} | |
151 | ||
152 | static inline T apply(T const&, T const&) | |
153 | { | |
154 | return T(1); | |
155 | } | |
156 | }; | |
157 | ||
158 | template <typename Type, | |
159 | bool IsFloatingPoint = boost::is_floating_point<Type>::value> | |
160 | struct equals | |
161 | { | |
162 | template <typename Policy> | |
163 | static inline bool apply(Type const& a, Type const& b, Policy const&) | |
164 | { | |
165 | return a == b; | |
166 | } | |
167 | }; | |
168 | ||
169 | template <typename Type> | |
170 | struct equals<Type, true> | |
171 | { | |
172 | template <typename Policy> | |
173 | static inline bool apply(Type const& a, Type const& b, Policy const& policy) | |
174 | { | |
175 | boost::ignore_unused(policy); | |
176 | ||
177 | if (a == b) | |
178 | { | |
179 | return true; | |
180 | } | |
181 | ||
182 | if (boost::math::isfinite(a) && boost::math::isfinite(b)) | |
183 | { | |
184 | // If a is INF and b is e.g. 0, the expression below returns true | |
185 | // but the values are obviously not equal, hence the condition | |
186 | return abs<Type>::apply(a - b) | |
187 | <= std::numeric_limits<Type>::epsilon() * policy.apply(a, b); | |
188 | } | |
189 | else | |
190 | { | |
191 | return a == b; | |
192 | } | |
193 | } | |
194 | }; | |
195 | ||
196 | template <typename T1, typename T2, typename Policy> | |
197 | inline bool equals_by_policy(T1 const& a, T2 const& b, Policy const& policy) | |
198 | { | |
199 | return detail::equals | |
200 | < | |
201 | typename select_most_precise<T1, T2>::type | |
202 | >::apply(a, b, policy); | |
203 | } | |
204 | ||
205 | template <typename Type, | |
206 | bool IsFloatingPoint = boost::is_floating_point<Type>::value> | |
207 | struct smaller | |
208 | { | |
209 | static inline bool apply(Type const& a, Type const& b) | |
210 | { | |
211 | return a < b; | |
212 | } | |
213 | }; | |
214 | ||
215 | template <typename Type> | |
216 | struct smaller<Type, true> | |
217 | { | |
218 | static inline bool apply(Type const& a, Type const& b) | |
219 | { | |
220 | if (!(a < b)) // a >= b | |
221 | { | |
222 | return false; | |
223 | } | |
224 | ||
225 | return ! equals<Type, true>::apply(b, a, equals_default_policy()); | |
226 | } | |
227 | }; | |
228 | ||
229 | template <typename Type, | |
230 | bool IsFloatingPoint = boost::is_floating_point<Type>::value> | |
231 | struct smaller_or_equals | |
232 | { | |
233 | static inline bool apply(Type const& a, Type const& b) | |
234 | { | |
235 | return a <= b; | |
236 | } | |
237 | }; | |
238 | ||
239 | template <typename Type> | |
240 | struct smaller_or_equals<Type, true> | |
241 | { | |
242 | static inline bool apply(Type const& a, Type const& b) | |
243 | { | |
244 | if (a <= b) | |
245 | { | |
246 | return true; | |
247 | } | |
248 | ||
249 | return equals<Type, true>::apply(a, b, equals_default_policy()); | |
250 | } | |
251 | }; | |
252 | ||
253 | ||
254 | template <typename Type, | |
255 | bool IsFloatingPoint = boost::is_floating_point<Type>::value> | |
256 | struct equals_with_epsilon | |
257 | : public equals<Type, IsFloatingPoint> | |
258 | {}; | |
259 | ||
260 | template | |
261 | < | |
262 | typename T, | |
263 | bool IsFundemantal = boost::is_fundamental<T>::value /* false */ | |
264 | > | |
265 | struct square_root | |
266 | { | |
267 | typedef T return_type; | |
268 | ||
269 | static inline T apply(T const& value) | |
270 | { | |
271 | // for non-fundamental number types assume that sqrt is | |
272 | // defined either: | |
273 | // 1) at T's scope, or | |
274 | // 2) at global scope, or | |
275 | // 3) in namespace std | |
276 | using ::sqrt; | |
277 | using std::sqrt; | |
278 | ||
279 | return sqrt(value); | |
280 | } | |
281 | }; | |
282 | ||
283 | template <typename FundamentalFP> | |
284 | struct square_root_for_fundamental_fp | |
285 | { | |
286 | typedef FundamentalFP return_type; | |
287 | ||
288 | static inline FundamentalFP apply(FundamentalFP const& value) | |
289 | { | |
290 | #ifdef BOOST_GEOMETRY_SQRT_CHECK_FINITENESS | |
291 | // This is a workaround for some 32-bit platforms. | |
292 | // For some of those platforms it has been reported that | |
293 | // std::sqrt(nan) and/or std::sqrt(-nan) returns a finite value. | |
294 | // For those platforms we need to define the macro | |
295 | // BOOST_GEOMETRY_SQRT_CHECK_FINITENESS so that the argument | |
296 | // to std::sqrt is checked appropriately before passed to std::sqrt | |
297 | if (boost::math::isfinite(value)) | |
298 | { | |
299 | return std::sqrt(value); | |
300 | } | |
301 | else if (boost::math::isinf(value) && value < 0) | |
302 | { | |
303 | return -std::numeric_limits<FundamentalFP>::quiet_NaN(); | |
304 | } | |
305 | return value; | |
306 | #else | |
307 | // for fundamental floating point numbers use std::sqrt | |
308 | return std::sqrt(value); | |
309 | #endif // BOOST_GEOMETRY_SQRT_CHECK_FINITENESS | |
310 | } | |
311 | }; | |
312 | ||
313 | template <> | |
314 | struct square_root<float, true> | |
315 | : square_root_for_fundamental_fp<float> | |
316 | { | |
317 | }; | |
318 | ||
319 | template <> | |
320 | struct square_root<double, true> | |
321 | : square_root_for_fundamental_fp<double> | |
322 | { | |
323 | }; | |
324 | ||
325 | template <> | |
326 | struct square_root<long double, true> | |
327 | : square_root_for_fundamental_fp<long double> | |
328 | { | |
329 | }; | |
330 | ||
331 | template <typename T> | |
332 | struct square_root<T, true> | |
333 | { | |
334 | typedef double return_type; | |
335 | ||
336 | static inline double apply(T const& value) | |
337 | { | |
338 | // for all other fundamental number types use also std::sqrt | |
339 | // | |
340 | // Note: in C++98 the only other possibility is double; | |
341 | // in C++11 there are also overloads for integral types; | |
342 | // this specialization works for those as well. | |
343 | return square_root_for_fundamental_fp | |
344 | < | |
345 | double | |
346 | >::apply(boost::numeric_cast<double>(value)); | |
347 | } | |
348 | }; | |
349 | ||
350 | ||
351 | ||
352 | template | |
353 | < | |
354 | typename T, | |
355 | bool IsFundemantal = boost::is_fundamental<T>::value /* false */ | |
356 | > | |
357 | struct modulo | |
358 | { | |
359 | typedef T return_type; | |
360 | ||
361 | static inline T apply(T const& value1, T const& value2) | |
362 | { | |
363 | // for non-fundamental number types assume that a free | |
364 | // function mod() is defined either: | |
365 | // 1) at T's scope, or | |
366 | // 2) at global scope | |
367 | return mod(value1, value2); | |
368 | } | |
369 | }; | |
370 | ||
371 | template | |
372 | < | |
373 | typename Fundamental, | |
374 | bool IsIntegral = boost::is_integral<Fundamental>::value | |
375 | > | |
376 | struct modulo_for_fundamental | |
377 | { | |
378 | typedef Fundamental return_type; | |
379 | ||
380 | static inline Fundamental apply(Fundamental const& value1, | |
381 | Fundamental const& value2) | |
382 | { | |
383 | return value1 % value2; | |
384 | } | |
385 | }; | |
386 | ||
387 | // specialization for floating-point numbers | |
388 | template <typename Fundamental> | |
389 | struct modulo_for_fundamental<Fundamental, false> | |
390 | { | |
391 | typedef Fundamental return_type; | |
392 | ||
393 | static inline Fundamental apply(Fundamental const& value1, | |
394 | Fundamental const& value2) | |
395 | { | |
396 | return std::fmod(value1, value2); | |
397 | } | |
398 | }; | |
399 | ||
400 | // specialization for fundamental number type | |
401 | template <typename Fundamental> | |
402 | struct modulo<Fundamental, true> | |
403 | : modulo_for_fundamental<Fundamental> | |
404 | {}; | |
405 | ||
406 | ||
407 | ||
408 | /*! | |
409 | \brief Short constructs to enable partial specialization for PI, 2*PI | |
410 | and PI/2, currently not possible in Math. | |
411 | */ | |
412 | template <typename T> | |
413 | struct define_pi | |
414 | { | |
415 | static inline T apply() | |
416 | { | |
417 | // Default calls Boost.Math | |
418 | return boost::math::constants::pi<T>(); | |
419 | } | |
420 | }; | |
421 | ||
422 | template <typename T> | |
423 | struct define_two_pi | |
424 | { | |
425 | static inline T apply() | |
426 | { | |
427 | // Default calls Boost.Math | |
428 | return boost::math::constants::two_pi<T>(); | |
429 | } | |
430 | }; | |
431 | ||
432 | template <typename T> | |
433 | struct define_half_pi | |
434 | { | |
435 | static inline T apply() | |
436 | { | |
437 | // Default calls Boost.Math | |
438 | return boost::math::constants::half_pi<T>(); | |
439 | } | |
440 | }; | |
441 | ||
442 | template <typename T> | |
443 | struct relaxed_epsilon | |
444 | { | |
445 | static inline T apply(const T& factor) | |
446 | { | |
447 | return factor * std::numeric_limits<T>::epsilon(); | |
448 | } | |
449 | }; | |
450 | ||
451 | // This must be consistent with math::equals. | |
452 | // By default math::equals() scales the error by epsilon using the greater of | |
453 | // compared values but here is only one value, though it should work the same way. | |
454 | // (a-a) <= max(a, a) * EPS -> 0 <= a*EPS | |
455 | // (a+da-a) <= max(a+da, a) * EPS -> da <= (a+da)*EPS | |
92f5a8d4 | 456 | template <typename T, bool IsIntegral = boost::is_integral<T>::value> |
7c673cae FG |
457 | struct scaled_epsilon |
458 | { | |
459 | static inline T apply(T const& val) | |
460 | { | |
461 | return (std::max)(abs<T>::apply(val), T(1)) | |
462 | * std::numeric_limits<T>::epsilon(); | |
463 | } | |
92f5a8d4 TL |
464 | |
465 | static inline T apply(T const& val, T const& eps) | |
466 | { | |
467 | return (std::max)(abs<T>::apply(val), T(1)) | |
468 | * eps; | |
469 | } | |
7c673cae FG |
470 | }; |
471 | ||
472 | template <typename T> | |
92f5a8d4 | 473 | struct scaled_epsilon<T, true> |
7c673cae FG |
474 | { |
475 | static inline T apply(T const&) | |
476 | { | |
477 | return T(0); | |
478 | } | |
92f5a8d4 TL |
479 | |
480 | static inline T apply(T const&, T const&) | |
481 | { | |
482 | return T(0); | |
483 | } | |
7c673cae FG |
484 | }; |
485 | ||
486 | // ItoF ItoI FtoF | |
487 | template <typename Result, typename Source, | |
488 | bool ResultIsInteger = std::numeric_limits<Result>::is_integer, | |
489 | bool SourceIsInteger = std::numeric_limits<Source>::is_integer> | |
490 | struct rounding_cast | |
491 | { | |
492 | static inline Result apply(Source const& v) | |
493 | { | |
494 | return boost::numeric_cast<Result>(v); | |
495 | } | |
496 | }; | |
497 | ||
498 | // TtoT | |
499 | template <typename Source, bool ResultIsInteger, bool SourceIsInteger> | |
500 | struct rounding_cast<Source, Source, ResultIsInteger, SourceIsInteger> | |
501 | { | |
502 | static inline Source apply(Source const& v) | |
503 | { | |
504 | return v; | |
505 | } | |
506 | }; | |
507 | ||
508 | // FtoI | |
509 | template <typename Result, typename Source> | |
510 | struct rounding_cast<Result, Source, true, false> | |
511 | { | |
512 | static inline Result apply(Source const& v) | |
513 | { | |
514 | return boost::numeric_cast<Result>(v < Source(0) ? | |
515 | v - Source(0.5) : | |
516 | v + Source(0.5)); | |
517 | } | |
518 | }; | |
519 | ||
520 | } // namespace detail | |
521 | #endif | |
522 | ||
523 | ||
524 | template <typename T> | |
525 | inline T pi() { return detail::define_pi<T>::apply(); } | |
526 | ||
527 | template <typename T> | |
528 | inline T two_pi() { return detail::define_two_pi<T>::apply(); } | |
529 | ||
530 | template <typename T> | |
531 | inline T half_pi() { return detail::define_half_pi<T>::apply(); } | |
532 | ||
533 | template <typename T> | |
534 | inline T relaxed_epsilon(T const& factor) | |
535 | { | |
536 | return detail::relaxed_epsilon<T>::apply(factor); | |
537 | } | |
538 | ||
539 | template <typename T> | |
540 | inline T scaled_epsilon(T const& value) | |
541 | { | |
542 | return detail::scaled_epsilon<T>::apply(value); | |
543 | } | |
544 | ||
92f5a8d4 TL |
545 | template <typename T> |
546 | inline T scaled_epsilon(T const& value, T const& eps) | |
547 | { | |
548 | return detail::scaled_epsilon<T>::apply(value, eps); | |
549 | } | |
7c673cae | 550 | |
b32b8144 | 551 | // Maybe replace this by boost equals or so |
7c673cae FG |
552 | |
553 | /*! | |
554 | \brief returns true if both arguments are equal. | |
555 | \ingroup utility | |
556 | \param a first argument | |
557 | \param b second argument | |
558 | \return true if a == b | |
559 | \note If both a and b are of an integral type, comparison is done by ==. | |
560 | If one of the types is floating point, comparison is done by abs and | |
561 | comparing with epsilon. If one of the types is non-fundamental, it might | |
562 | be a high-precision number and comparison is done using the == operator | |
563 | of that class. | |
564 | */ | |
565 | ||
566 | template <typename T1, typename T2> | |
567 | inline bool equals(T1 const& a, T2 const& b) | |
568 | { | |
569 | return detail::equals | |
570 | < | |
571 | typename select_most_precise<T1, T2>::type | |
572 | >::apply(a, b, detail::equals_default_policy()); | |
573 | } | |
574 | ||
575 | template <typename T1, typename T2> | |
576 | inline bool equals_with_epsilon(T1 const& a, T2 const& b) | |
577 | { | |
578 | return detail::equals_with_epsilon | |
579 | < | |
580 | typename select_most_precise<T1, T2>::type | |
581 | >::apply(a, b, detail::equals_default_policy()); | |
582 | } | |
583 | ||
584 | template <typename T1, typename T2> | |
585 | inline bool smaller(T1 const& a, T2 const& b) | |
586 | { | |
587 | return detail::smaller | |
588 | < | |
589 | typename select_most_precise<T1, T2>::type | |
590 | >::apply(a, b); | |
591 | } | |
592 | ||
593 | template <typename T1, typename T2> | |
594 | inline bool larger(T1 const& a, T2 const& b) | |
595 | { | |
596 | return detail::smaller | |
597 | < | |
598 | typename select_most_precise<T1, T2>::type | |
599 | >::apply(b, a); | |
600 | } | |
601 | ||
602 | template <typename T1, typename T2> | |
603 | inline bool smaller_or_equals(T1 const& a, T2 const& b) | |
604 | { | |
605 | return detail::smaller_or_equals | |
606 | < | |
607 | typename select_most_precise<T1, T2>::type | |
608 | >::apply(a, b); | |
609 | } | |
610 | ||
611 | template <typename T1, typename T2> | |
612 | inline bool larger_or_equals(T1 const& a, T2 const& b) | |
613 | { | |
614 | return detail::smaller_or_equals | |
615 | < | |
616 | typename select_most_precise<T1, T2>::type | |
617 | >::apply(b, a); | |
618 | } | |
619 | ||
620 | ||
621 | template <typename T> | |
622 | inline T d2r() | |
623 | { | |
624 | static T const conversion_coefficient = geometry::math::pi<T>() / T(180.0); | |
625 | return conversion_coefficient; | |
626 | } | |
627 | ||
628 | template <typename T> | |
629 | inline T r2d() | |
630 | { | |
631 | static T const conversion_coefficient = T(180.0) / geometry::math::pi<T>(); | |
632 | return conversion_coefficient; | |
633 | } | |
634 | ||
635 | ||
636 | #ifndef DOXYGEN_NO_DETAIL | |
637 | namespace detail { | |
638 | ||
639 | template <typename DegreeOrRadian> | |
640 | struct as_radian | |
641 | { | |
642 | template <typename T> | |
643 | static inline T apply(T const& value) | |
644 | { | |
645 | return value; | |
646 | } | |
647 | }; | |
648 | ||
649 | template <> | |
650 | struct as_radian<degree> | |
651 | { | |
652 | template <typename T> | |
653 | static inline T apply(T const& value) | |
654 | { | |
655 | return value * d2r<T>(); | |
656 | } | |
657 | }; | |
658 | ||
659 | template <typename DegreeOrRadian> | |
660 | struct from_radian | |
661 | { | |
662 | template <typename T> | |
663 | static inline T apply(T const& value) | |
664 | { | |
665 | return value; | |
666 | } | |
667 | }; | |
668 | ||
669 | template <> | |
670 | struct from_radian<degree> | |
671 | { | |
672 | template <typename T> | |
673 | static inline T apply(T const& value) | |
674 | { | |
675 | return value * r2d<T>(); | |
676 | } | |
677 | }; | |
678 | ||
679 | } // namespace detail | |
680 | #endif | |
681 | ||
682 | template <typename DegreeOrRadian, typename T> | |
683 | inline T as_radian(T const& value) | |
684 | { | |
685 | return detail::as_radian<DegreeOrRadian>::apply(value); | |
686 | } | |
687 | ||
688 | template <typename DegreeOrRadian, typename T> | |
689 | inline T from_radian(T const& value) | |
690 | { | |
691 | return detail::from_radian<DegreeOrRadian>::apply(value); | |
692 | } | |
693 | ||
694 | ||
695 | /*! | |
696 | \brief Calculates the haversine of an angle | |
697 | \ingroup utility | |
698 | \note See http://en.wikipedia.org/wiki/Haversine_formula | |
699 | haversin(alpha) = sin2(alpha/2) | |
700 | */ | |
701 | template <typename T> | |
702 | inline T hav(T const& theta) | |
703 | { | |
704 | T const half = T(0.5); | |
705 | T const sn = sin(half * theta); | |
706 | return sn * sn; | |
707 | } | |
708 | ||
709 | /*! | |
710 | \brief Short utility to return the square | |
711 | \ingroup utility | |
712 | \param value Value to calculate the square from | |
713 | \return The squared value | |
714 | */ | |
715 | template <typename T> | |
716 | inline T sqr(T const& value) | |
717 | { | |
718 | return value * value; | |
719 | } | |
720 | ||
721 | /*! | |
722 | \brief Short utility to return the square root | |
723 | \ingroup utility | |
724 | \param value Value to calculate the square root from | |
725 | \return The square root value | |
726 | */ | |
727 | template <typename T> | |
728 | inline typename detail::square_root<T>::return_type | |
729 | sqrt(T const& value) | |
730 | { | |
731 | return detail::square_root | |
732 | < | |
733 | T, boost::is_fundamental<T>::value | |
734 | >::apply(value); | |
735 | } | |
736 | ||
737 | /*! | |
738 | \brief Short utility to return the modulo of two values | |
739 | \ingroup utility | |
740 | \param value1 First value | |
741 | \param value2 Second value | |
742 | \return The result of the modulo operation on the (ordered) pair | |
743 | (value1, value2) | |
744 | */ | |
745 | template <typename T> | |
746 | inline typename detail::modulo<T>::return_type | |
747 | mod(T const& value1, T const& value2) | |
748 | { | |
749 | return detail::modulo | |
750 | < | |
751 | T, boost::is_fundamental<T>::value | |
752 | >::apply(value1, value2); | |
753 | } | |
754 | ||
755 | /*! | |
756 | \brief Short utility to workaround gcc/clang problem that abs is converting to integer | |
757 | and that older versions of MSVC does not support abs of long long... | |
758 | \ingroup utility | |
759 | */ | |
760 | template<typename T> | |
761 | inline T abs(T const& value) | |
762 | { | |
763 | return detail::abs<T>::apply(value); | |
764 | } | |
765 | ||
766 | /*! | |
767 | \brief Short utility to calculate the sign of a number: -1 (negative), 0 (zero), 1 (positive) | |
768 | \ingroup utility | |
769 | */ | |
770 | template <typename T> | |
771 | inline int sign(T const& value) | |
772 | { | |
773 | T const zero = T(); | |
774 | return value > zero ? 1 : value < zero ? -1 : 0; | |
775 | } | |
776 | ||
777 | /*! | |
778 | \brief Short utility to cast a value possibly rounding it to the nearest | |
779 | integral value. | |
780 | \ingroup utility | |
781 | \note If the source T is NOT an integral type and Result is an integral type | |
782 | the value is rounded towards the closest integral value. Otherwise it's | |
783 | casted without rounding. | |
784 | */ | |
785 | template <typename Result, typename T> | |
786 | inline Result rounding_cast(T const& v) | |
787 | { | |
788 | return detail::rounding_cast<Result, T>::apply(v); | |
789 | } | |
790 | ||
92f5a8d4 TL |
791 | /*! |
792 | \brief Evaluate the sine and cosine function with the argument in degrees | |
793 | \note The results obey exactly the elementary properties of the trigonometric | |
794 | functions, e.g., sin 9° = cos 81° = − sin 123456789°. | |
795 | If x = −0, then \e sinx = −0; this is the only case where | |
796 | −0 is returned. | |
797 | */ | |
798 | template<typename T> | |
799 | inline void sin_cos_degrees(T const& x, | |
800 | T & sinx, | |
801 | T & cosx) | |
802 | { | |
803 | // In order to minimize round-off errors, this function exactly reduces | |
804 | // the argument to the range [-45, 45] before converting it to radians. | |
805 | T remainder; int quotient; | |
806 | ||
807 | remainder = math::mod(x, T(360)); | |
808 | quotient = floor(remainder / 90 + T(0.5)); | |
809 | remainder -= 90 * quotient; | |
810 | ||
811 | // Convert to radians. | |
812 | remainder *= d2r<T>(); | |
813 | ||
814 | T s = sin(remainder), c = cos(remainder); | |
815 | ||
816 | switch (unsigned(quotient) & 3U) | |
817 | { | |
818 | case 0U: sinx = s; cosx = c; break; | |
819 | case 1U: sinx = c; cosx = -s; break; | |
820 | case 2U: sinx = -s; cosx = -c; break; | |
821 | default: sinx = -c; cosx = s; break; // case 3U | |
822 | } | |
823 | ||
824 | // Set sign of 0 results. -0 only produced for sin(-0). | |
825 | if (x != 0) | |
826 | { | |
827 | sinx += T(0); cosx += T(0); | |
828 | } | |
829 | } | |
830 | ||
831 | /*! | |
832 | \brief Round off a given angle | |
833 | */ | |
834 | template<typename T> | |
835 | inline T round_angle(T const& x) { | |
836 | static const T z = 1/T(16); | |
837 | ||
838 | if (x == 0) | |
839 | { | |
840 | return 0; | |
841 | } | |
842 | ||
843 | T y = math::abs(x); | |
844 | ||
845 | // z - (z - y) must not be simplified to y. | |
846 | y = y < z ? z - (z - y) : y; | |
847 | ||
848 | return x < 0 ? -y : y; | |
849 | } | |
850 | ||
851 | ||
852 | /*! | |
853 | \brief The error-free sum of two numbers. | |
854 | */ | |
855 | template<typename T> | |
856 | inline T sum_error(T const& u, T const& v, T& t) | |
857 | { | |
858 | volatile T s = u + v; | |
859 | volatile T up = s - v; | |
860 | volatile T vpp = s - up; | |
861 | ||
862 | up -= u; | |
863 | vpp -= v; | |
864 | t = -(up + vpp); | |
865 | ||
866 | return s; | |
867 | } | |
868 | ||
869 | /*! | |
870 | \brief Evaluate the polynomial in x using Horner's method. | |
871 | */ | |
872 | // TODO: adl1995 - Merge these functions with formulas/area_formulas.hpp | |
873 | // i.e. place them in one file. | |
874 | template <typename NT, typename IteratorType> | |
875 | inline NT horner_evaluate(NT const& x, | |
876 | IteratorType begin, | |
877 | IteratorType end) | |
878 | { | |
879 | NT result(0); | |
880 | IteratorType it = end; | |
881 | do | |
882 | { | |
883 | result = result * x + *--it; | |
884 | } | |
885 | while (it != begin); | |
886 | return result; | |
887 | } | |
888 | ||
889 | /*! | |
890 | \brief Evaluate the polynomial. | |
891 | */ | |
892 | template<typename IteratorType, typename CT> | |
893 | inline CT polyval(IteratorType first, | |
894 | IteratorType last, | |
895 | CT const& eps) | |
896 | { | |
897 | int N = std::distance(first, last) - 1; | |
898 | int index = 0; | |
899 | ||
900 | CT y = N < 0 ? 0 : *(first + (index++)); | |
901 | ||
902 | while (--N >= 0) | |
903 | { | |
904 | y = y * eps + *(first + (index++)); | |
905 | } | |
906 | ||
907 | return y; | |
908 | } | |
909 | ||
910 | /* | |
911 | \brief Short utility to calculate the power | |
912 | \ingroup utility | |
913 | */ | |
914 | template <typename T1, typename T2> | |
915 | inline T1 pow(T1 const& a, T2 const& b) | |
916 | { | |
917 | using std::pow; | |
918 | return pow(a, b); | |
919 | } | |
920 | ||
7c673cae FG |
921 | } // namespace math |
922 | ||
923 | ||
924 | }} // namespace boost::geometry | |
925 | ||
926 | #endif // BOOST_GEOMETRY_UTIL_MATH_HPP |