machar.py 11 KB

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