Welcome, guest | Sign In | My Account | Store | Cart

Arbitrary precision binary floating point arithmetic (Python recipe) by Raymond Hettinger
ActiveState Code (http://code.activestate.com/recipes/412428/)

Achieves exactness by manipulating values stored as a long integer mantissa with an integer exponent (base two).

Python, 112 lines
Download
  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)

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.

Tags: algorithms