Achieves exactness by manipulating values stored as a long integer mantissa with an integer exponent (base two).
| Python |
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 | from math import frexp
class LF(object):
"""Long Float -- a class for high-precision floating point arithmetic.
All arithmetic operations are exact except for division which has
regular floating point precision.
Construct with:
* integers: LF(10)
* floats: LF(math.pi)
* tuples: LF((323, -5)) # equal to 323 * 2.0 ** -5
"""
__slots__ = ['mant', 'exp'] # long integer mantissa and int exponent
def __new__(cls, value):
if isinstance(value, LF):
return value
self = object.__new__(cls)
if isinstance(value, tuple):
mant, exp = value
assert int(mant)==mant and int(exp)==exp
else:
mant, exp = frexp(value)
mant, exp = int(mant * 2.0 ** 53), exp-53
while mant and not mant & 1:
mant >>= 1
exp += 1
self.mant, self.exp = mant, exp
return self
def __float__(self):
return float(str(self.mant)) * 2.0 ** self.exp
def __int__(self):
m, e = self.mant, self.exp
if e >= 0:
value = m << e
else:
value = m >> -e
return int(value)
def __long__(self):
return long(int(self))
def __repr__(self):
return '%s((%d, %d))' % (self.__class__.__name__, self.mant, self.exp)
def __str__(self):
return str(float(self))
def _sync(self, other):
other = LF(other)
smant, sexp, omant, oexp = self.mant, self.exp, other.mant, other.exp
if sexp < oexp:
omant <<= oexp - sexp
oexp = sexp
else:
smant <<= sexp - oexp
return smant, omant, oexp
def __add__(self, other):
smant, omant, exp = self._sync(other)
return LF((smant + omant, exp))
def __radd__(self, other):
return LF(self) + other
def __sub__(self, other):
smant, omant, exp = self._sync(other)
return LF((smant - omant, exp))
def __rsub__(self, other):
return LF(self) - other
def __mul__(self, other):
other = LF(other)
return LF((self.mant*other.mant, self.exp+other.exp))
def __rmul__(self, other):
return LF(self) * other
def __cmp__(self, other):
smant, omant, exp = self._sync(other)
return cmp(smant, omant)
def __hash__(self):
return hash(float(self))
def __pow__(self, b):
assert b == int(b) and b >= 0
return LF((self.mant ** b, self.exp * b))
def __abs__(self):
return LF((abs(self.mant), self.exp))
def __pos__(self):
return self
def __neg__(self):
return LF((-self.mant, self.exp))
def __div__(self, other):
return LF(float(self) / float(other))
def __rdiv__(self, other):
return LF(float(self) / float(other))
def __nonzero__(self):
return bool(self.mant)
|
Discussion
All arithmetic is exact except for division (which is carried out using regular floating point arithmetic).
Construction from integers and floats is exact (i.e. LF(10) or LF(math.pi)).
Coercion back to a regular float is rounded to within about 1/2 unit in the last decimal place.
Sign in to comment