1"""
2Machine arithmetic - determine the parameters of the
3floating-point arithmetic system
4
5Author: Pearu Peterson, September 2003
6
7"""
8__all__ = ['MachAr']
9
10from .fromnumeric import any
11from ._ufunc_config import errstate
12from .._utils import set_module
13
14# Need to speed this up...especially for longdouble
15
16# Deprecated 2021-10-20, NumPy 1.22
17class MachAr:
18 """
19 Diagnosing machine parameters.
20
21 Attributes
22 ----------
23 ibeta : int
24 Radix in which numbers are represented.
25 it : int
26 Number of base-`ibeta` digits in the floating point mantissa M.
27 machep : int
28 Exponent of the smallest (most negative) power of `ibeta` that,
29 added to 1.0, gives something different from 1.0
30 eps : float
31 Floating-point number ``beta**machep`` (floating point precision)
32 negep : int
33 Exponent of the smallest power of `ibeta` that, subtracted
34 from 1.0, gives something different from 1.0.
35 epsneg : float
36 Floating-point number ``beta**negep``.
37 iexp : int
38 Number of bits in the exponent (including its sign and bias).
39 minexp : int
40 Smallest (most negative) power of `ibeta` consistent with there
41 being no leading zeros in the mantissa.
42 xmin : float
43 Floating-point number ``beta**minexp`` (the smallest [in
44 magnitude] positive floating point number with full precision).
45 maxexp : int
46 Smallest (positive) power of `ibeta` that causes overflow.
47 xmax : float
48 ``(1-epsneg) * beta**maxexp`` (the largest [in magnitude]
49 usable floating value).
50 irnd : int
51 In ``range(6)``, information on what kind of rounding is done
52 in addition, and on how underflow is handled.
53 ngrd : int
54 Number of 'guard digits' used when truncating the product
55 of two mantissas to fit the representation.
56 epsilon : float
57 Same as `eps`.
58 tiny : float
59 An alias for `smallest_normal`, kept for backwards compatibility.
60 huge : float
61 Same as `xmax`.
62 precision : float
63 ``- int(-log10(eps))``
64 resolution : float
65 ``- 10**(-precision)``
66 smallest_normal : float
67 The smallest positive floating point number with 1 as leading bit in
68 the mantissa following IEEE-754. Same as `xmin`.
69 smallest_subnormal : float
70 The smallest positive floating point number with 0 as leading bit in
71 the mantissa following IEEE-754.
72
73 Parameters
74 ----------
75 float_conv : function, optional
76 Function that converts an integer or integer array to a float
77 or float array. Default is `float`.
78 int_conv : function, optional
79 Function that converts a float or float array to an integer or
80 integer array. Default is `int`.
81 float_to_float : function, optional
82 Function that converts a float array to float. Default is `float`.
83 Note that this does not seem to do anything useful in the current
84 implementation.
85 float_to_str : function, optional
86 Function that converts a single float to a string. Default is
87 ``lambda v:'%24.16e' %v``.
88 title : str, optional
89 Title that is printed in the string representation of `MachAr`.
90
91 See Also
92 --------
93 finfo : Machine limits for floating point types.
94 iinfo : Machine limits for integer types.
95
96 References
97 ----------
98 .. [1] Press, Teukolsky, Vetterling and Flannery,
99 "Numerical Recipes in C++," 2nd ed,
100 Cambridge University Press, 2002, p. 31.
101
102 """
103
104 def __init__(self, float_conv=float,int_conv=int,
105 float_to_float=float,
106 float_to_str=lambda v:'%24.16e' % v,
107 title='Python floating point number'):
108 """
109
110 float_conv - convert integer to float (array)
111 int_conv - convert float (array) to integer
112 float_to_float - convert float array to float
113 float_to_str - convert array float to str
114 title - description of used floating point numbers
115
116 """
117 # We ignore all errors here because we are purposely triggering
118 # underflow to detect the properties of the running arch.
119 with errstate(under='ignore'):
120 self._do_init(float_conv, int_conv, float_to_float, float_to_str, title)
121
122 def _do_init(self, float_conv, int_conv, float_to_float, float_to_str, title):
123 max_iterN = 10000
124 msg = "Did not converge after %d tries with %s"
125 one = float_conv(1)
126 two = one + one
127 zero = one - one
128
129 # Do we really need to do this? Aren't they 2 and 2.0?
130 # Determine ibeta and beta
131 a = one
132 for _ in range(max_iterN):
133 a = a + a
134 temp = a + one
135 temp1 = temp - a
136 if any(temp1 - one != zero):
137 break
138 else:
139 raise RuntimeError(msg % (_, one.dtype))
140 b = one
141 for _ in range(max_iterN):
142 b = b + b
143 temp = a + b
144 itemp = int_conv(temp-a)
145 if any(itemp != 0):
146 break
147 else:
148 raise RuntimeError(msg % (_, one.dtype))
149 ibeta = itemp
150 beta = float_conv(ibeta)
151
152 # Determine it and irnd
153 it = -1
154 b = one
155 for _ in range(max_iterN):
156 it = it + 1
157 b = b * beta
158 temp = b + one
159 temp1 = temp - b
160 if any(temp1 - one != zero):
161 break
162 else:
163 raise RuntimeError(msg % (_, one.dtype))
164
165 betah = beta / two
166 a = one
167 for _ in range(max_iterN):
168 a = a + a
169 temp = a + one
170 temp1 = temp - a
171 if any(temp1 - one != zero):
172 break
173 else:
174 raise RuntimeError(msg % (_, one.dtype))
175 temp = a + betah
176 irnd = 0
177 if any(temp-a != zero):
178 irnd = 1
179 tempa = a + beta
180 temp = tempa + betah
181 if irnd == 0 and any(temp-tempa != zero):
182 irnd = 2
183
184 # Determine negep and epsneg
185 negep = it + 3
186 betain = one / beta
187 a = one
188 for i in range(negep):
189 a = a * betain
190 b = a
191 for _ in range(max_iterN):
192 temp = one - a
193 if any(temp-one != zero):
194 break
195 a = a * beta
196 negep = negep - 1
197 # Prevent infinite loop on PPC with gcc 4.0:
198 if negep < 0:
199 raise RuntimeError("could not determine machine tolerance "
200 "for 'negep', locals() -> %s" % (locals()))
201 else:
202 raise RuntimeError(msg % (_, one.dtype))
203 negep = -negep
204 epsneg = a
205
206 # Determine machep and eps
207 machep = - it - 3
208 a = b
209
210 for _ in range(max_iterN):
211 temp = one + a
212 if any(temp-one != zero):
213 break
214 a = a * beta
215 machep = machep + 1
216 else:
217 raise RuntimeError(msg % (_, one.dtype))
218 eps = a
219
220 # Determine ngrd
221 ngrd = 0
222 temp = one + eps
223 if irnd == 0 and any(temp*one - one != zero):
224 ngrd = 1
225
226 # Determine iexp
227 i = 0
228 k = 1
229 z = betain
230 t = one + eps
231 nxres = 0
232 for _ in range(max_iterN):
233 y = z
234 z = y*y
235 a = z*one # Check here for underflow
236 temp = z*t
237 if any(a+a == zero) or any(abs(z) >= y):
238 break
239 temp1 = temp * betain
240 if any(temp1*beta == z):
241 break
242 i = i + 1
243 k = k + k
244 else:
245 raise RuntimeError(msg % (_, one.dtype))
246 if ibeta != 10:
247 iexp = i + 1
248 mx = k + k
249 else:
250 iexp = 2
251 iz = ibeta
252 while k >= iz:
253 iz = iz * ibeta
254 iexp = iexp + 1
255 mx = iz + iz - 1
256
257 # Determine minexp and xmin
258 for _ in range(max_iterN):
259 xmin = y
260 y = y * betain
261 a = y * one
262 temp = y * t
263 if any((a + a) != zero) and any(abs(y) < xmin):
264 k = k + 1
265 temp1 = temp * betain
266 if any(temp1*beta == y) and any(temp != y):
267 nxres = 3
268 xmin = y
269 break
270 else:
271 break
272 else:
273 raise RuntimeError(msg % (_, one.dtype))
274 minexp = -k
275
276 # Determine maxexp, xmax
277 if mx <= k + k - 3 and ibeta != 10:
278 mx = mx + mx
279 iexp = iexp + 1
280 maxexp = mx + minexp
281 irnd = irnd + nxres
282 if irnd >= 2:
283 maxexp = maxexp - 2
284 i = maxexp + minexp
285 if ibeta == 2 and not i:
286 maxexp = maxexp - 1
287 if i > 20:
288 maxexp = maxexp - 1
289 if any(a != y):
290 maxexp = maxexp - 2
291 xmax = one - epsneg
292 if any(xmax*one != xmax):
293 xmax = one - beta*epsneg
294 xmax = xmax / (xmin*beta*beta*beta)
295 i = maxexp + minexp + 3
296 for j in range(i):
297 if ibeta == 2:
298 xmax = xmax + xmax
299 else:
300 xmax = xmax * beta
301
302 smallest_subnormal = abs(xmin / beta ** (it))
303
304 self.ibeta = ibeta
305 self.it = it
306 self.negep = negep
307 self.epsneg = float_to_float(epsneg)
308 self._str_epsneg = float_to_str(epsneg)
309 self.machep = machep
310 self.eps = float_to_float(eps)
311 self._str_eps = float_to_str(eps)
312 self.ngrd = ngrd
313 self.iexp = iexp
314 self.minexp = minexp
315 self.xmin = float_to_float(xmin)
316 self._str_xmin = float_to_str(xmin)
317 self.maxexp = maxexp
318 self.xmax = float_to_float(xmax)
319 self._str_xmax = float_to_str(xmax)
320 self.irnd = irnd
321
322 self.title = title
323 # Commonly used parameters
324 self.epsilon = self.eps
325 self.tiny = self.xmin
326 self.huge = self.xmax
327 self.smallest_normal = self.xmin
328 self._str_smallest_normal = float_to_str(self.xmin)
329 self.smallest_subnormal = float_to_float(smallest_subnormal)
330 self._str_smallest_subnormal = float_to_str(smallest_subnormal)
331
332 import math
333 self.precision = int(-math.log10(float_to_float(self.eps)))
334 ten = two + two + two + two + two
335 resolution = ten ** (-self.precision)
336 self.resolution = float_to_float(resolution)
337 self._str_resolution = float_to_str(resolution)
338
339 def __str__(self):
340 fmt = (
341 'Machine parameters for %(title)s\n'
342 '---------------------------------------------------------------------\n'
343 'ibeta=%(ibeta)s it=%(it)s iexp=%(iexp)s ngrd=%(ngrd)s irnd=%(irnd)s\n'
344 'machep=%(machep)s eps=%(_str_eps)s (beta**machep == epsilon)\n'
345 'negep =%(negep)s epsneg=%(_str_epsneg)s (beta**epsneg)\n'
346 'minexp=%(minexp)s xmin=%(_str_xmin)s (beta**minexp == tiny)\n'
347 'maxexp=%(maxexp)s xmax=%(_str_xmax)s ((1-epsneg)*beta**maxexp == huge)\n'
348 'smallest_normal=%(smallest_normal)s '
349 'smallest_subnormal=%(smallest_subnormal)s\n'
350 '---------------------------------------------------------------------\n'
351 )
352 return fmt % self.__dict__
353
354
355if __name__ == '__main__':
356 print(MachAr())