geomath.py
Go to the documentation of this file.
1 """geomath.py: transcription of GeographicLib::Math class."""
2 # geomath.py
3 #
4 # This is a rather literal translation of the GeographicLib::Math class to
5 # python. See the documentation for the C++ class for more information at
6 #
7 # https://geographiclib.sourceforge.io/html/annotated.html
8 #
9 # Copyright (c) Charles Karney (2011-2017) <charles@karney.com> and
10 # licensed under the MIT/X11 License. For more information, see
11 # https://geographiclib.sourceforge.io/
12 
13 
14 import sys
15 import math
16 
17 class Math(object):
18  """
19  Additional math routines for GeographicLib.
20 
21  This defines constants:
22  epsilon, difference between 1 and the next bigger number
23  digits, the number of digits in the fraction of a real number
24  minval, minimum normalized positive number
25  maxval, maximum finite number
26  nan, not a number
27  inf, infinity
28  """
29 
30  digits = 53
31  epsilon = math.pow(2.0, 1-digits)
32  minval = math.pow(2.0, -1022)
33  maxval = math.pow(2.0, 1023) * (2 - epsilon)
34  inf = float("inf") if sys.version_info > (2, 6) else 2 * maxval
35  nan = float("nan") if sys.version_info > (2, 6) else inf - inf
36 
37  def sq(x):
38  """Square a number"""
39 
40  return x * x
41  sq = staticmethod(sq)
42 
43  def cbrt(x):
44  """Real cube root of a number"""
45 
46  y = math.pow(abs(x), 1/3.0)
47  return y if x >= 0 else -y
48  cbrt = staticmethod(cbrt)
49 
50  def log1p(x):
51  """log(1 + x) accurate for small x (missing from python 2.5.2)"""
52 
53  if sys.version_info > (2, 6):
54  return math.log1p(x)
55 
56  y = 1 + x
57  z = y - 1
58  # Here's the explanation for this magic: y = 1 + z, exactly, and z
59  # approx x, thus log(y)/z (which is nearly constant near z = 0) returns
60  # a good approximation to the true log(1 + x)/x. The multiplication x *
61  # (log(y)/z) introduces little additional error.
62  return x if z == 0 else x * math.log(y) / z
63  log1p = staticmethod(log1p)
64 
65  def atanh(x):
66  """atanh(x) (missing from python 2.5.2)"""
67 
68  if sys.version_info > (2, 6):
69  return math.atanh(x)
70 
71  y = abs(x) # Enforce odd parity
72  y = Math.log1p(2 * y/(1 - y))/2
73  return -y if x < 0 else y
74  atanh = staticmethod(atanh)
75 
76  def copysign(x, y):
77  """return x with the sign of y (missing from python 2.5.2)"""
78 
79  if sys.version_info > (2, 6):
80  return math.copysign(x, y)
81 
82  return math.fabs(x) * (-1 if y < 0 or (y == 0 and 1/y < 0) else 1)
83  copysign = staticmethod(copysign)
84 
85  def norm(x, y):
86  """Private: Normalize a two-vector."""
87  r = math.hypot(x, y)
88  return x/r, y/r
89  norm = staticmethod(norm)
90 
91  def sum(u, v):
92  """Error free transformation of a sum."""
93  # Error free transformation of a sum. Note that t can be the same as one
94  # of the first two arguments.
95  s = u + v
96  up = s - v
97  vpp = s - up
98  up -= u
99  vpp -= v
100  t = -(up + vpp)
101  # u + v = s + t
102  # = round(u + v) + t
103  return s, t
104  sum = staticmethod(sum)
105 
106  def polyval(N, p, s, x):
107  """Evaluate a polynomial."""
108  y = float(0 if N < 0 else p[s]) # make sure the returned value is a float
109  while N > 0:
110  N -= 1; s += 1
111  y = y * x + p[s]
112  return y
113  polyval = staticmethod(polyval)
114 
115  def AngRound(x):
116  """Private: Round an angle so that small values underflow to zero."""
117  # The makes the smallest gap in x = 1/16 - nextafter(1/16, 0) = 1/2^57
118  # for reals = 0.7 pm on the earth if x is an angle in degrees. (This
119  # is about 1000 times more resolution than we get with angles around 90
120  # degrees.) We use this to avoid having to deal with near singular
121  # cases when x is non-zero but tiny (e.g., 1.0e-200).
122  z = 1/16.0
123  y = abs(x)
124  # The compiler mustn't "simplify" z - (z - y) to y
125  if y < z: y = z - (z - y)
126  return 0.0 if x == 0 else (-y if x < 0 else y)
127  AngRound = staticmethod(AngRound)
128 
129  def AngNormalize(x):
130  """reduce angle to (-180,180]"""
131 
132  y = math.fmod(x, 360)
133  # On Windows 32-bit with python 2.7, math.fmod(-0.0, 360) = +0.0
134  # This fixes this bug. See also Math::AngNormalize in the C++ library.
135  # sincosd has a similar fix.
136  y = x if x == 0 else y
137  return (y + 360 if y <= -180 else
138  (y if y <= 180 else y - 360))
139  AngNormalize = staticmethod(AngNormalize)
140 
141  def LatFix(x):
142  """replace angles outside [-90,90] by NaN"""
143 
144  return Math.nan if abs(x) > 90 else x
145  LatFix = staticmethod(LatFix)
146 
147  def AngDiff(x, y):
148  """compute y - x and reduce to [-180,180] accurately"""
149 
150  d, t = Math.sum(Math.AngNormalize(-x), Math.AngNormalize(y))
151  d = Math.AngNormalize(d)
152  return Math.sum(-180 if d == 180 and t > 0 else d, t)
153  AngDiff = staticmethod(AngDiff)
154 
155  def sincosd(x):
156  """Compute sine and cosine of x in degrees."""
157 
158  r = math.fmod(x, 360)
159  q = Math.nan if Math.isnan(r) else int(math.floor(r / 90 + 0.5))
160  r -= 90 * q; r = math.radians(r)
161  s = math.sin(r); c = math.cos(r)
162  q = q % 4
163  if q == 1:
164  s, c = c, -s
165  elif q == 2:
166  s, c = -s, -c
167  elif q == 3:
168  s, c = -c, s
169  # Remove the minus sign on -0.0 except for sin(-0.0).
170  # On Windows 32-bit with python 2.7, math.fmod(-0.0, 360) = +0.0
171  # (x, c) here fixes this bug. See also Math::sincosd in the C++ library.
172  # AngNormalize has a similar fix.
173  s, c = (x, c) if x == 0 else (0.0+s, 0.0+c)
174  return s, c
175  sincosd = staticmethod(sincosd)
176 
177  def atan2d(y, x):
178  """compute atan2(y, x) with the result in degrees"""
179 
180  if abs(y) > abs(x):
181  q = 2; x, y = y, x
182  else:
183  q = 0
184  if x < 0:
185  q += 1; x = -x
186  ang = math.degrees(math.atan2(y, x))
187  if q == 1:
188  ang = (180 if y >= 0 else -180) - ang
189  elif q == 2:
190  ang = 90 - ang
191  elif q == 3:
192  ang = -90 + ang
193  return ang
194  atan2d = staticmethod(atan2d)
195 
196  def isfinite(x):
197  """Test for finiteness"""
198 
199  return abs(x) <= Math.maxval
200  isfinite = staticmethod(isfinite)
201 
202  def isnan(x):
203  """Test if nan"""
204 
205  return math.isnan(x) if sys.version_info > (2, 6) else x != x
206  isnan = staticmethod(isnan)
#define abs(x)
Definition: datatypes.h:17


gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:34:18