GeographicLib  1.51
Math.cpp
Go to the documentation of this file.
1 /**
2  * \file Math.cpp
3  * \brief Implementation for GeographicLib::Math class
4  *
5  * Copyright (c) Charles Karney (2015-2020) <charles@karney.com> and licensed
6  * under the MIT/X11 License. For more information, see
7  * https://geographiclib.sourceforge.io/
8  **********************************************************************/
9 
10 #include <GeographicLib/Math.hpp>
11 
12 #if defined(_MSC_VER)
13 // Squelch warnings about constant conditional expressions
14 # pragma warning (disable: 4127)
15 #endif
16 
17 namespace GeographicLib {
18 
19  using namespace std;
20 
21  void Math::dummy() {
22  static_assert(GEOGRAPHICLIB_PRECISION >= 1 && GEOGRAPHICLIB_PRECISION <= 5,
23  "Bad value of precision");
24  }
25 
26  int Math::digits() {
27 #if GEOGRAPHICLIB_PRECISION != 5
28  return numeric_limits<real>::digits;
29 #else
30  return numeric_limits<real>::digits();
31 #endif
32  }
33 
34  int Math::set_digits(int ndigits) {
35 #if GEOGRAPHICLIB_PRECISION != 5
36  (void)ndigits;
37 #else
38  mpfr::mpreal::set_default_prec(ndigits >= 2 ? ndigits : 2);
39 #endif
40  return digits();
41  }
42 
44 #if GEOGRAPHICLIB_PRECISION != 5
45  return numeric_limits<real>::digits10;
46 #else
47  return numeric_limits<real>::digits10();
48 #endif
49  }
50 
52  return
53  digits10() > numeric_limits<double>::digits10 ?
54  digits10() - numeric_limits<double>::digits10 : 0;
55  }
56 
57  template<typename T> T Math::hypot(T x, T y) {
58  using std::hypot; return hypot(x, y);
59  }
60 
61  template<typename T> T Math::expm1(T x) {
62  using std::expm1; return expm1(x);
63  }
64 
65  template<typename T> T Math::log1p(T x) {
66  using std::log1p; return log1p(x);
67  }
68 
69  template<typename T> T Math::asinh(T x) {
70  using std::asinh; return asinh(x);
71  }
72 
73  template<typename T> T Math::atanh(T x) {
74  using std::atanh; return atanh(x);
75  }
76 
77  template<typename T> T Math::copysign(T x, T y) {
78  using std::copysign; return copysign(x, y);
79  }
80 
81  template<typename T> T Math::cbrt(T x) {
82  using std::cbrt; return cbrt(x);
83  }
84 
85  template<typename T> T Math::remainder(T x, T y) {
86  using std::remainder; return remainder(x, y);
87  }
88 
89  template<typename T> T Math::remquo(T x, T y, int* n) {
90  using std::remquo; return remquo(x, y, n);
91  }
92 
93  template<typename T> T Math::round(T x) {
94  using std::round; return round(x);
95  }
96 
97  template<typename T> long Math::lround(T x) {
98  using std::lround; return lround(x);
99  }
100 
101  template<typename T> T Math::fma(T x, T y, T z) {
102  using std::fma; return fma(x, y, z);
103  }
104 
105  template<typename T> T Math::sum(T u, T v, T& t) {
106  GEOGRAPHICLIB_VOLATILE T s = u + v;
107  GEOGRAPHICLIB_VOLATILE T up = s - v;
108  GEOGRAPHICLIB_VOLATILE T vpp = s - up;
109  up -= u;
110  vpp -= v;
111  t = -(up + vpp);
112  // u + v = s + t
113  // = round(u + v) + t
114  return s;
115  }
116 
117  template<typename T> T Math::AngRound(T x) {
118  static const T z = 1/T(16);
119  if (x == 0) return 0;
120  GEOGRAPHICLIB_VOLATILE T y = abs(x);
121  // The compiler mustn't "simplify" z - (z - y) to y
122  y = y < z ? z - (z - y) : y;
123  return x < 0 ? -y : y;
124  }
125 
126  template<typename T> void Math::sincosd(T x, T& sinx, T& cosx) {
127  // In order to minimize round-off errors, this function exactly reduces
128  // the argument to the range [-45, 45] before converting it to radians.
129  using std::remquo;
130  T r; int q;
131  // N.B. the implementation of remquo in glibc pre 2.22 were buggy. See
132  // https://sourceware.org/bugzilla/show_bug.cgi?id=17569
133  // This was fixed in version 2.22 on 2015-08-05
134  r = remquo(x, T(90), &q); // now abs(r) <= 45
135  r *= degree<T>();
136  // g++ -O turns these two function calls into a call to sincos
137  T s = sin(r), c = cos(r);
138 #if defined(_MSC_VER) && _MSC_VER < 1900
139  // Before version 14 (2015), Visual Studio had problems dealing
140  // with -0.0. Specifically
141  // VC 10,11,12 and 32-bit compile: fmod(-0.0, 360.0) -> +0.0
142  // VC 12 and 64-bit compile: sin(-0.0) -> +0.0
143  // AngNormalize has a similar fix.
144  // python 2.7 on Windows 32-bit machines has the same problem.
145  if (x == 0) s = x;
146 #endif
147  switch (unsigned(q) & 3U) {
148  case 0U: sinx = s; cosx = c; break;
149  case 1U: sinx = c; cosx = -s; break;
150  case 2U: sinx = -s; cosx = -c; break;
151  default: sinx = -c; cosx = s; break; // case 3U
152  }
153  // Set sign of 0 results. -0 only produced for sin(-0)
154  if (x != 0) { sinx += T(0); cosx += T(0); }
155  }
156 
157  template<typename T> T Math::sind(T x) {
158  // See sincosd
159  using std::remquo;
160  T r; int q;
161  r = remquo(x, T(90), &q); // now abs(r) <= 45
162  r *= degree<T>();
163  unsigned p = unsigned(q);
164  r = p & 1U ? cos(r) : sin(r);
165  if (p & 2U) r = -r;
166  if (x != 0) r += T(0);
167  return r;
168  }
169 
170  template<typename T> T Math::cosd(T x) {
171  // See sincosd
172  using std::remquo;
173  T r; int q;
174  r = remquo(x, T(90), &q); // now abs(r) <= 45
175  r *= degree<T>();
176  unsigned p = unsigned(q + 1);
177  r = p & 1U ? cos(r) : sin(r);
178  if (p & 2U) r = -r;
179  return T(0) + r;
180  }
181 
182  template<typename T> T Math::tand(T x) {
183  static const T overflow = 1 / sq(numeric_limits<T>::epsilon());
184  T s, c;
185  sincosd(x, s, c);
186  return c != 0 ? s / c : (s < 0 ? -overflow : overflow);
187  }
188 
189  template<typename T> T Math::atan2d(T y, T x) {
190  // In order to minimize round-off errors, this function rearranges the
191  // arguments so that result of atan2 is in the range [-pi/4, pi/4] before
192  // converting it to degrees and mapping the result to the correct
193  // quadrant.
194  int q = 0;
195  if (abs(y) > abs(x)) { swap(x, y); q = 2; }
196  if (x < 0) { x = -x; ++q; }
197  // here x >= 0 and x >= abs(y), so angle is in [-pi/4, pi/4]
198  T ang = atan2(y, x) / degree<T>();
199  switch (q) {
200  // Note that atan2d(-0.0, 1.0) will return -0. However, we expect that
201  // atan2d will not be called with y = -0. If need be, include
202  //
203  // case 0: ang = 0 + ang; break;
204  //
205  // and handle mpfr as in AngRound.
206  case 1: ang = (y >= 0 ? 180 : -180) - ang; break;
207  case 2: ang = 90 - ang; break;
208  case 3: ang = -90 + ang; break;
209  }
210  return ang;
211  }
212 
213  template<typename T> T Math::atand(T x)
214  { return atan2d(x, T(1)); }
215 
216  template<typename T> T Math::eatanhe(T x, T es) {
217  using std::atanh;
218  return es > T(0) ? es * atanh(es * x) : -es * atan(es * x);
219  }
220 
221  template<typename T> T Math::taupf(T tau, T es) {
222  // Need this test, otherwise tau = +/-inf gives taup = nan.
223  using std::isfinite; using std::hypot;
224  if (isfinite(tau)) {
225  T tau1 = hypot(T(1), tau),
226  sig = sinh( eatanhe(tau / tau1, es ) );
227  return hypot(T(1), sig) * tau - sig * tau1;
228  } else
229  return tau;
230  }
231 
232  template<typename T> T Math::tauf(T taup, T es) {
233  using std::hypot;
234  static const int numit = 5;
235  // min iterations = 1, max iterations = 2; mean = 1.95
236  static const T tol = sqrt(numeric_limits<T>::epsilon()) / 10;
237  static const T taumax = 2 / sqrt(numeric_limits<T>::epsilon());
238  T e2m = T(1) - sq(es),
239  // To lowest order in e^2, taup = (1 - e^2) * tau = _e2m * tau; so use
240  // tau = taup/e2m as a starting guess. Only 1 iteration is needed for
241  // |lat| < 3.35 deg, otherwise 2 iterations are needed. If, instead, tau
242  // = taup is used the mean number of iterations increases to 1.999 (2
243  // iterations are needed except near tau = 0).
244  //
245  // For large tau, taup = exp(-es*atanh(es)) * tau. Use this as for the
246  // initial guess for |taup| > 70 (approx |phi| > 89deg). Then for
247  // sufficiently large tau (such that sqrt(1+tau^2) = |tau|), we can exit
248  // with the intial guess and avoid overflow problems. This also reduces
249  // the mean number of iterations slightly from 1.963 to 1.954.
250  tau = abs(taup) > 70 ? taup * exp(eatanhe(T(1), es)) : taup/e2m,
251  stol = tol * max(T(1), abs(taup));
252  if (!(abs(tau) < taumax)) return tau; // handles +/-inf and nan
253  for (int i = 0; i < numit || GEOGRAPHICLIB_PANIC; ++i) {
254  T taupa = taupf(tau, es),
255  dtau = (taup - taupa) * (1 + e2m * sq(tau)) /
256  ( e2m * hypot(T(1), tau) * hypot(T(1), taupa) );
257  tau += dtau;
258  if (!(abs(dtau) >= stol))
259  break;
260  }
261  return tau;
262  }
263 
264  template<typename T> bool Math::isfinite(T x) {
265  using std::isfinite; return isfinite(x);
266  }
267 
268  template<typename T> T Math::NaN() {
269 #if defined(_MSC_VER)
270  return numeric_limits<T>::has_quiet_NaN ?
271  numeric_limits<T>::quiet_NaN() :
272  (numeric_limits<T>::max)();
273 #else
274  return numeric_limits<T>::has_quiet_NaN ?
275  numeric_limits<T>::quiet_NaN() :
276  numeric_limits<T>::max();
277 #endif
278  }
279 
280  template<typename T> bool Math::isnan(T x) {
281  using std::isnan; return isnan(x);
282  }
283 
284  template<typename T> T Math::infinity() {
285 #if defined(_MSC_VER)
286  return numeric_limits<T>::has_infinity ?
287  numeric_limits<T>::infinity() :
288  (numeric_limits<T>::max)();
289 #else
290  return numeric_limits<T>::has_infinity ?
291  numeric_limits<T>::infinity() :
292  numeric_limits<T>::max();
293 #endif
294  }
295 
296  /// \cond SKIP
297  // Instantiate
298 #define GEOGRAPHICLIB_MATH_INSTANTIATE(T) \
299  template T GEOGRAPHICLIB_EXPORT Math::hypot <T>(T, T); \
300  template T GEOGRAPHICLIB_EXPORT Math::expm1 <T>(T); \
301  template T GEOGRAPHICLIB_EXPORT Math::log1p <T>(T); \
302  template T GEOGRAPHICLIB_EXPORT Math::asinh <T>(T); \
303  template T GEOGRAPHICLIB_EXPORT Math::atanh <T>(T); \
304  template T GEOGRAPHICLIB_EXPORT Math::cbrt <T>(T); \
305  template T GEOGRAPHICLIB_EXPORT Math::remainder<T>(T, T); \
306  template T GEOGRAPHICLIB_EXPORT Math::remquo <T>(T, T, int*); \
307  template T GEOGRAPHICLIB_EXPORT Math::round <T>(T); \
308  template long GEOGRAPHICLIB_EXPORT Math::lround <T>(T); \
309  template T GEOGRAPHICLIB_EXPORT Math::copysign <T>(T, T); \
310  template T GEOGRAPHICLIB_EXPORT Math::fma <T>(T, T, T); \
311  template T GEOGRAPHICLIB_EXPORT Math::sum <T>(T, T, T&); \
312  template T GEOGRAPHICLIB_EXPORT Math::AngRound <T>(T); \
313  template void GEOGRAPHICLIB_EXPORT Math::sincosd <T>(T, T&, T&); \
314  template T GEOGRAPHICLIB_EXPORT Math::sind <T>(T); \
315  template T GEOGRAPHICLIB_EXPORT Math::cosd <T>(T); \
316  template T GEOGRAPHICLIB_EXPORT Math::tand <T>(T); \
317  template T GEOGRAPHICLIB_EXPORT Math::atan2d <T>(T, T); \
318  template T GEOGRAPHICLIB_EXPORT Math::atand <T>(T); \
319  template T GEOGRAPHICLIB_EXPORT Math::eatanhe <T>(T, T); \
320  template T GEOGRAPHICLIB_EXPORT Math::taupf <T>(T, T); \
321  template T GEOGRAPHICLIB_EXPORT Math::tauf <T>(T, T); \
322  template bool GEOGRAPHICLIB_EXPORT Math::isfinite <T>(T); \
323  template T GEOGRAPHICLIB_EXPORT Math::NaN <T>(); \
324  template bool GEOGRAPHICLIB_EXPORT Math::isnan <T>(T); \
325  template T GEOGRAPHICLIB_EXPORT Math::infinity <T>();
326 
327  // Instantiate with the standard floating type
328  GEOGRAPHICLIB_MATH_INSTANTIATE(float)
329  GEOGRAPHICLIB_MATH_INSTANTIATE(double)
330 #if GEOGRAPHICLIB_HAVE_LONG_DOUBLE
331  // Instantiate if long double is distinct from double
332  GEOGRAPHICLIB_MATH_INSTANTIATE(long double)
333 #endif
334 #if GEOGRAPHICLIB_PRECISION > 3
335  // Instantiate with the high precision type
336  GEOGRAPHICLIB_MATH_INSTANTIATE(Math::real)
337 #endif
338 
339 #undef GEOGRAPHICLIB_MATH_INSTANTIATE
340 
341  // Also we need int versions for Utility::nummatch
342  template int GEOGRAPHICLIB_EXPORT Math::NaN <int>();
343  template int GEOGRAPHICLIB_EXPORT Math::infinity<int>();
344  /// \endcond
345 
346 } // namespace GeographicLib
static T atand(T x)
Definition: Math.cpp:213
static T tand(T x)
Definition: Math.cpp:182
#define GEOGRAPHICLIB_EXPORT
Definition: Constants.hpp:66
static bool isfinite(T x)
Definition: Math.cpp:264
static T infinity()
Definition: Math.cpp:284
static T log1p(T x)
Definition: Math.cpp:65
static long lround(T x)
Definition: Math.cpp:97
#define GEOGRAPHICLIB_PRECISION
Definition: Math.hpp:35
#define GEOGRAPHICLIB_VOLATILE
Definition: Math.hpp:58
static int extra_digits()
Definition: Math.cpp:51
static T atan2d(T y, T x)
Definition: Math.cpp:189
Header for GeographicLib::Math class.
static int set_digits(int ndigits)
Definition: Math.cpp:34
static T cbrt(T x)
Definition: Math.cpp:81
static T round(T x)
Definition: Math.cpp:93
Namespace for GeographicLib.
Definition: Accumulator.cpp:12
static T NaN()
Definition: Math.cpp:268
static T sum(T u, T v, T &t)
Definition: Math.cpp:105
static T atanh(T x)
Definition: Math.cpp:73
static bool isnan(T x)
Definition: Math.cpp:280
static T asinh(T x)
Definition: Math.cpp:69
void swap(GeographicLib::NearestNeighbor< dist_t, pos_t, distfun_t > &a, GeographicLib::NearestNeighbor< dist_t, pos_t, distfun_t > &b)
static T copysign(T x, T y)
Definition: Math.cpp:77
static T sind(T x)
Definition: Math.cpp:157
static T cosd(T x)
Definition: Math.cpp:170
static T remquo(T x, T y, int *n)
Definition: Math.cpp:89
static T hypot(T x, T y)
Definition: Math.cpp:57
static int digits10()
Definition: Math.cpp:43
static T tauf(T taup, T es)
Definition: Math.cpp:232
static void sincosd(T x, T &sinx, T &cosx)
Definition: Math.cpp:126
static T eatanhe(T x, T es)
Definition: Math.cpp:216
static int digits()
Definition: Math.cpp:26
static T fma(T x, T y, T z)
Definition: Math.cpp:101
static T expm1(T x)
Definition: Math.cpp:61
static T taupf(T tau, T es)
Definition: Math.cpp:221
static T remainder(T x, T y)
Definition: Math.cpp:85
#define GEOGRAPHICLIB_PANIC
Definition: Math.hpp:61
static T AngRound(T x)
Definition: Math.cpp:117