""" Machine arithmetics - determine the parameters of the floating-point arithmetic system """ # Author: Pearu Peterson, September 2003 __all__ = ['MachAr'] from numpy.core.fromnumeric import any from numpy.core.numeric import seterr # Need to speed this up...especially for longfloat class MachAr(object): """Diagnosing machine parameters. The following attributes are available: ibeta - radix in which numbers are represented it - number of base-ibeta digits in the floating point mantissa M machep - exponent of the smallest (most negative) power of ibeta that, added to 1.0, gives something different from 1.0 eps - floating-point number beta**machep (floating point precision) negep - exponent of the smallest power of ibeta that, substracted from 1.0, gives something different from 1.0 epsneg - floating-point number beta**negep iexp - number of bits in the exponent (including its sign and bias) minexp - smallest (most negative) power of ibeta consistent with there being no leading zeros in the mantissa xmin - floating point number beta**minexp (the smallest (in magnitude) usable floating value) maxexp - smallest (positive) power of ibeta that causes overflow xmax - (1-epsneg)* beta**maxexp (the largest (in magnitude) usable floating value) irnd - in range(6), information on what kind of rounding is done in addition, and on how underflow is handled ngrd - number of 'guard digits' used when truncating the product of two mantissas to fit the representation epsilon - same as eps tiny - same as xmin huge - same as xmax precision - int(-log10(eps)) resolution - 10**(-precision) Reference: Numerical Recipies. """ def __init__(self, float_conv=float,int_conv=int, float_to_float=float, float_to_str = lambda v:'%24.16e' % v, title = 'Python floating point number'): """ float_conv - convert integer to float (array) int_conv - convert float (array) to integer float_to_float - convert float array to float float_to_str - convert array float to str title - description of used floating point numbers """ # We ignore all errors here because we are purposely triggering # underflow to detect the properties of the runninng arch. saverrstate = seterr(under='ignore') try: self._do_init(float_conv, int_conv, float_to_float, float_to_str, title) finally: seterr(**saverrstate) def _do_init(self, float_conv, int_conv, float_to_float, float_to_str, title): max_iterN = 10000 msg = "Did not converge after %d tries with %s" one = float_conv(1) two = one + one zero = one - one # Do we really need to do this? Aren't they 2 and 2.0? # Determine ibeta and beta a = one for _ in xrange(max_iterN): a = a + a temp = a + one temp1 = temp - a if any(temp1 - one != zero): break else: raise RuntimeError, msg % (_, one.dtype) b = one for _ in xrange(max_iterN): b = b + b temp = a + b itemp = int_conv(temp-a) if any(itemp != 0): break else: raise RuntimeError, msg % (_, one.dtype) ibeta = itemp beta = float_conv(ibeta) # Determine it and irnd it = -1 b = one for _ in xrange(max_iterN): it = it + 1 b = b * beta temp = b + one temp1 = temp - b if any(temp1 - one != zero): break else: raise RuntimeError, msg % (_, one.dtype) betah = beta / two a = one for _ in xrange(max_iterN): a = a + a temp = a + one temp1 = temp - a if any(temp1 - one != zero): break else: raise RuntimeError, msg % (_, one.dtype) temp = a + betah irnd = 0 if any(temp-a != zero): irnd = 1 tempa = a + beta temp = tempa + betah if irnd==0 and any(temp-tempa != zero): irnd = 2 # Determine negep and epsneg negep = it + 3 betain = one / beta a = one for i in range(negep): a = a * betain b = a for _ in xrange(max_iterN): temp = one - a if any(temp-one != zero): break a = a * beta negep = negep - 1 # Prevent infinite loop on PPC with gcc 4.0: if negep < 0: raise RuntimeError, "could not determine machine tolerance " \ "for 'negep', locals() -> %s" % (locals()) else: raise RuntimeError, msg % (_, one.dtype) negep = -negep epsneg = a # Determine machep and eps machep = - it - 3 a = b for _ in xrange(max_iterN): temp = one + a if any(temp-one != zero): break a = a * beta machep = machep + 1 else: raise RuntimeError, msg % (_, one.dtype) eps = a # Determine ngrd ngrd = 0 temp = one + eps if irnd==0 and any(temp*one - one != zero): ngrd = 1 # Determine iexp i = 0 k = 1 z = betain t = one + eps nxres = 0 for _ in xrange(max_iterN): y = z z = y*y a = z*one # Check here for underflow temp = z*t if any(a+a == zero) or any(abs(z)>=y): break temp1 = temp * betain if any(temp1*beta == z): break i = i + 1 k = k + k else: raise RuntimeError, msg % (_, one.dtype) if ibeta != 10: iexp = i + 1 mx = k + k else: iexp = 2 iz = ibeta while k >= iz: iz = iz * ibeta iexp = iexp + 1 mx = iz + iz - 1 # Determine minexp and xmin for _ in xrange(max_iterN): xmin = y y = y * betain a = y * one temp = y * t if any(a+a != zero) and any(abs(y) < xmin): k = k + 1 temp1 = temp * betain if any(temp1*beta == y) and any(temp != y): nxres = 3 xmin = y break else: break else: raise RuntimeError, msg % (_, one.dtype) minexp = -k # Determine maxexp, xmax if mx <= k + k - 3 and ibeta != 10: mx = mx + mx iexp = iexp + 1 maxexp = mx + minexp irnd = irnd + nxres if irnd >= 2: maxexp = maxexp - 2 i = maxexp + minexp if ibeta == 2 and not i: maxexp = maxexp - 1 if i > 20: maxexp = maxexp - 1 if any(a != y): maxexp = maxexp - 2 xmax = one - epsneg if any(xmax*one != xmax): xmax = one - beta*epsneg xmax = xmax / (xmin*beta*beta*beta) i = maxexp + minexp + 3 for j in range(i): if ibeta==2: xmax = xmax + xmax else: xmax = xmax * beta self.ibeta = ibeta self.it = it self.negep = negep self.epsneg = float_to_float(epsneg) self._str_epsneg = float_to_str(epsneg) self.machep = machep self.eps = float_to_float(eps) self._str_eps = float_to_str(eps) self.ngrd = ngrd self.iexp = iexp self.minexp = minexp self.xmin = float_to_float(xmin) self._str_xmin = float_to_str(xmin) self.maxexp = maxexp self.xmax = float_to_float(xmax) self._str_xmax = float_to_str(xmax) self.irnd = irnd self.title = title # Commonly used parameters self.epsilon = self.eps self.tiny = self.xmin self.huge = self.xmax import math self.precision = int(-math.log10(float_to_float(self.eps))) ten = two + two + two + two + two resolution = ten ** (-self.precision) self.resolution = float_to_float(resolution) self._str_resolution = float_to_str(resolution) def __str__(self): return '''\ Machine parameters for %(title)s --------------------------------------------------------------------- ibeta=%(ibeta)s it=%(it)s iexp=%(iexp)s ngrd=%(ngrd)s irnd=%(irnd)s machep=%(machep)s eps=%(_str_eps)s (beta**machep == epsilon) negep =%(negep)s epsneg=%(_str_epsneg)s (beta**epsneg) minexp=%(minexp)s xmin=%(_str_xmin)s (beta**minexp == tiny) maxexp=%(maxexp)s xmax=%(_str_xmax)s ((1-epsneg)*beta**maxexp == huge) --------------------------------------------------------------------- ''' % self.__dict__ if __name__ == '__main__': print MachAr()