_machar.py 11 KB

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