00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #include <GeographicLib/Constants.hpp>
00013
00014 #if !defined(GEOGRAPHICLIB_MATH_HPP)
00015 #define GEOGRAPHICLIB_MATH_HPP 1
00016
00017
00018
00019
00020 #if !defined(GEOGRAPHICLIB_CXX11_MATH)
00021
00022
00023
00024
00025
00026 # if defined(__GNUC__) && __GNUC__ == 4 && __GNUC_MINOR__ >= 7 \
00027 && __cplusplus >= 201103 && \
00028 !(defined(__ANDROID__) || defined(ANDROID) || defined(__CYGWIN__))
00029 # define GEOGRAPHICLIB_CXX11_MATH 1
00030
00031 # elif defined(_MSC_VER) && _MSC_VER >= 1800
00032 # define GEOGRAPHICLIB_CXX11_MATH 1
00033 # else
00034 # define GEOGRAPHICLIB_CXX11_MATH 0
00035 # endif
00036 #endif
00037
00038 #if !defined(GEOGRAPHICLIB_WORDS_BIGENDIAN)
00039 # define GEOGRAPHICLIB_WORDS_BIGENDIAN 0
00040 #endif
00041
00042 #if !defined(GEOGRAPHICLIB_HAVE_LONG_DOUBLE)
00043 # define GEOGRAPHICLIB_HAVE_LONG_DOUBLE 0
00044 #endif
00045
00046 #if !defined(GEOGRAPHICLIB_PRECISION)
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056 # define GEOGRAPHICLIB_PRECISION 2
00057 #endif
00058
00059 #include <cmath>
00060 #include <algorithm>
00061 #include <limits>
00062
00063 #if GEOGRAPHICLIB_PRECISION == 4
00064 #include <boost/multiprecision/float128.hpp>
00065 #include <boost/math/special_functions/hypot.hpp>
00066 #include <boost/math/special_functions/expm1.hpp>
00067 #include <boost/math/special_functions/log1p.hpp>
00068 #include <boost/math/special_functions/atanh.hpp>
00069 #include <boost/math/special_functions/asinh.hpp>
00070 #include <boost/math/special_functions/cbrt.hpp>
00071 #elif GEOGRAPHICLIB_PRECISION == 5
00072 #include <mpreal.h>
00073 #endif
00074
00075 #if GEOGRAPHICLIB_PRECISION > 3
00076
00077 #define GEOGRAPHICLIB_VOLATILE
00078
00079
00080 #define GEOGRAPHICLIB_PANIC \
00081 (throw GeographicLib::GeographicErr("Convergence failure"), false)
00082 #else
00083 #define GEOGRAPHICLIB_VOLATILE volatile
00084
00085
00086 #define GEOGRAPHICLIB_PANIC false
00087 #endif
00088
00089 namespace GeographicLib {
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101 class GEOGRAPHICLIB_EXPORT Math {
00102 private:
00103 void dummy() {
00104 GEOGRAPHICLIB_STATIC_ASSERT(GEOGRAPHICLIB_PRECISION >= 1 &&
00105 GEOGRAPHICLIB_PRECISION <= 5,
00106 "Bad value of precision");
00107 }
00108 Math();
00109 public:
00110
00111 #if GEOGRAPHICLIB_HAVE_LONG_DOUBLE
00112
00113
00114
00115
00116 typedef long double extended;
00117 #else
00118 typedef double extended;
00119 #endif
00120
00121 #if GEOGRAPHICLIB_PRECISION == 2
00122
00123
00124
00125
00126
00127
00128 typedef double real;
00129 #elif GEOGRAPHICLIB_PRECISION == 1
00130 typedef float real;
00131 #elif GEOGRAPHICLIB_PRECISION == 3
00132 typedef extended real;
00133 #elif GEOGRAPHICLIB_PRECISION == 4
00134 typedef boost::multiprecision::float128 real;
00135 #elif GEOGRAPHICLIB_PRECISION == 5
00136 typedef mpfr::mpreal real;
00137 #else
00138 typedef double real;
00139 #endif
00140
00141
00142
00143
00144 static inline int digits() {
00145 #if GEOGRAPHICLIB_PRECISION != 5
00146 return std::numeric_limits<real>::digits;
00147 #else
00148 return std::numeric_limits<real>::digits();
00149 #endif
00150 }
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160 static inline int set_digits(int ndigits) {
00161 #if GEOGRAPHICLIB_PRECISION != 5
00162 (void)ndigits;
00163 #else
00164 mpfr::mpreal::set_default_prec(ndigits >= 2 ? ndigits : 2);
00165 #endif
00166 return digits();
00167 }
00168
00169
00170
00171
00172 static inline int digits10() {
00173 #if GEOGRAPHICLIB_PRECISION != 5
00174 return std::numeric_limits<real>::digits10;
00175 #else
00176 return std::numeric_limits<real>::digits10();
00177 #endif
00178 }
00179
00180
00181
00182
00183
00184 static inline int extra_digits() {
00185 return
00186 digits10() > std::numeric_limits<double>::digits10 ?
00187 digits10() - std::numeric_limits<double>::digits10 : 0;
00188 }
00189
00190 #if GEOGRAPHICLIB_PRECISION <= 3
00191
00192
00193
00194
00195
00196
00197 static const int extradigits =
00198 std::numeric_limits<real>::digits10 >
00199 std::numeric_limits<double>::digits10 ?
00200 std::numeric_limits<real>::digits10 -
00201 std::numeric_limits<double>::digits10 : 0;
00202 #endif
00203
00204
00205
00206
00207 static const bool bigendian = GEOGRAPHICLIB_WORDS_BIGENDIAN;
00208
00209
00210
00211
00212
00213 template<typename T> static inline T pi() {
00214 using std::atan2;
00215 static const T pi = atan2(T(0), T(-1));
00216 return pi;
00217 }
00218
00219
00220
00221 static inline real pi() { return pi<real>(); }
00222
00223
00224
00225
00226
00227 template<typename T> static inline T degree() {
00228 static const T degree = pi<T>() / 180;
00229 return degree;
00230 }
00231
00232
00233
00234 static inline real degree() { return degree<real>(); }
00235
00236
00237
00238
00239
00240
00241
00242
00243 template<typename T> static inline T sq(T x)
00244 { return x * x; }
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254 template<typename T> static inline T hypot(T x, T y) {
00255 #if GEOGRAPHICLIB_CXX11_MATH
00256 using std::hypot; return hypot(x, y);
00257 #else
00258 using std::abs; using std::sqrt;
00259 x = abs(x); y = abs(y);
00260 if (x < y) std::swap(x, y);
00261 y /= (x ? x : 1);
00262 return x * sqrt(1 + y * y);
00263
00264
00265
00266 #endif
00267 }
00268
00269
00270
00271
00272
00273
00274
00275
00276 template<typename T> static inline T expm1(T x) {
00277 #if GEOGRAPHICLIB_CXX11_MATH
00278 using std::expm1; return expm1(x);
00279 #else
00280 using std::exp; using std::abs; using std::log;
00281 volatile T
00282 y = exp(x),
00283 z = y - 1;
00284
00285
00286
00287
00288 return abs(x) > 1 ? z : (z == 0 ? x : x * z / log(y));
00289 #endif
00290 }
00291
00292
00293
00294
00295
00296
00297
00298
00299 template<typename T> static inline T log1p(T x) {
00300 #if GEOGRAPHICLIB_CXX11_MATH
00301 using std::log1p; return log1p(x);
00302 #else
00303 using std::log;
00304 volatile T
00305 y = 1 + x,
00306 z = y - 1;
00307
00308
00309
00310
00311 return z == 0 ? x : x * log(y) / z;
00312 #endif
00313 }
00314
00315
00316
00317
00318
00319
00320
00321
00322 template<typename T> static inline T asinh(T x) {
00323 #if GEOGRAPHICLIB_CXX11_MATH
00324 using std::asinh; return asinh(x);
00325 #else
00326 using std::abs; T y = abs(x);
00327 y = log1p(y * (1 + y/(hypot(T(1), y) + 1)));
00328 return x < 0 ? -y : y;
00329 #endif
00330 }
00331
00332
00333
00334
00335
00336
00337
00338
00339 template<typename T> static inline T atanh(T x) {
00340 #if GEOGRAPHICLIB_CXX11_MATH
00341 using std::atanh; return atanh(x);
00342 #else
00343 using std::abs; T y = abs(x);
00344 y = log1p(2 * y/(1 - y))/2;
00345 return x < 0 ? -y : y;
00346 #endif
00347 }
00348
00349
00350
00351
00352
00353
00354
00355
00356 template<typename T> static inline T cbrt(T x) {
00357 #if GEOGRAPHICLIB_CXX11_MATH
00358 using std::cbrt; return cbrt(x);
00359 #else
00360 using std::abs; using std::pow;
00361 T y = pow(abs(x), 1/T(3));
00362 return x < 0 ? -y : y;
00363 #endif
00364 }
00365
00366
00367
00368
00369
00370
00371
00372
00373
00374
00375
00376
00377
00378 template<typename T> static inline T sum(T u, T v, T& t) {
00379 GEOGRAPHICLIB_VOLATILE T s = u + v;
00380 GEOGRAPHICLIB_VOLATILE T up = s - v;
00381 GEOGRAPHICLIB_VOLATILE T vpp = s - up;
00382 up -= u;
00383 vpp -= v;
00384 t = -(up + vpp);
00385
00386
00387 return s;
00388 }
00389
00390
00391
00392
00393
00394
00395
00396
00397
00398
00399 template<typename T> static inline T AngNormalize(T x)
00400 { return x >= 180 ? x - 360 : (x < -180 ? x + 360 : x); }
00401
00402
00403
00404
00405
00406
00407
00408
00409
00410
00411 template<typename T> static inline T AngNormalize2(T x)
00412 { using std::fmod; return AngNormalize<T>(fmod(x, T(360))); }
00413
00414
00415
00416
00417
00418
00419
00420
00421
00422
00423
00424
00425
00426
00427
00428
00429 template<typename T> static inline T AngDiff(T x, T y) {
00430 T t, d = sum(-x, y, t);
00431 if ((d - T(180)) + t > T(0))
00432 d -= T(360);
00433 else if ((d + T(180)) + t <= T(0))
00434 d += T(360);
00435 return d + t;
00436 }
00437
00438
00439
00440
00441
00442
00443
00444
00445 template<typename T> static inline bool isfinite(T x) {
00446 #if GEOGRAPHICLIB_CXX11_MATH
00447 using std::isfinite; return isfinite(x);
00448 #else
00449 using std::abs;
00450 return abs(x) <= (std::numeric_limits<T>::max)();
00451 #endif
00452 }
00453
00454
00455
00456
00457
00458
00459
00460 template<typename T> static inline T NaN() {
00461 return std::numeric_limits<T>::has_quiet_NaN ?
00462 std::numeric_limits<T>::quiet_NaN() :
00463 (std::numeric_limits<T>::max)();
00464 }
00465
00466
00467
00468 static inline real NaN() { return NaN<real>(); }
00469
00470
00471
00472
00473
00474
00475
00476
00477 template<typename T> static inline bool isnan(T x) {
00478 #if GEOGRAPHICLIB_CXX11_MATH
00479 using std::isnan; return isnan(x);
00480 #else
00481 return x != x;
00482 #endif
00483 }
00484
00485
00486
00487
00488
00489
00490
00491 template<typename T> static inline T infinity() {
00492 return std::numeric_limits<T>::has_infinity ?
00493 std::numeric_limits<T>::infinity() :
00494 (std::numeric_limits<T>::max)();
00495 }
00496
00497
00498
00499 static inline real infinity() { return infinity<real>(); }
00500
00501
00502
00503
00504
00505
00506
00507
00508 template<typename T> static inline T swab(T x) {
00509 union {
00510 T r;
00511 unsigned char c[sizeof(T)];
00512 } b;
00513 b.r = x;
00514 for (int i = sizeof(T)/2; i--; )
00515 std::swap(b.c[i], b.c[sizeof(T) - 1 - i]);
00516 return b.r;
00517 }
00518
00519 #if GEOGRAPHICLIB_PRECISION == 4
00520 typedef boost::math::policies::policy
00521 < boost::math::policies::domain_error
00522 <boost::math::policies::errno_on_error>,
00523 boost::math::policies::pole_error
00524 <boost::math::policies::errno_on_error>,
00525 boost::math::policies::overflow_error
00526 <boost::math::policies::errno_on_error>,
00527 boost::math::policies::evaluation_error
00528 <boost::math::policies::errno_on_error> >
00529 boost_special_functions_policy;
00530
00531 static inline real hypot(real x, real y)
00532 { return boost::math::hypot(x, y, boost_special_functions_policy()); }
00533
00534 static inline real expm1(real x)
00535 { return boost::math::expm1(x, boost_special_functions_policy()); }
00536
00537 static inline real log1p(real x)
00538 { return boost::math::log1p(x, boost_special_functions_policy()); }
00539
00540 static inline real asinh(real x)
00541 { return boost::math::asinh(x, boost_special_functions_policy()); }
00542
00543 static inline real atanh(real x)
00544 { return boost::math::atanh(x, boost_special_functions_policy()); }
00545
00546 static inline real cbrt(real x)
00547 { return boost::math::cbrt(x, boost_special_functions_policy()); }
00548
00549 static inline bool isnan(real x) { return boost::math::isnan(x); }
00550
00551 static inline bool isfinite(real x) { return boost::math::isfinite(x); }
00552 #endif
00553 };
00554
00555 }
00556
00557 #endif // GEOGRAPHICLIB_MATH_HPP