machar.py 11 KB

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