New python module to simulate arbitrary fixed and infinite precision binary floating point

[Rob Clewley]

Dear Pythonistas,

How many times have we seen posts recently along the lines of "why is
it that 0.1 appears as 0.10000000000000001 in python?" that lead to
posters being sent to the definition of the IEEE 754 standard and the
decimal.py module? I am teaching an introductory numerical analysis
class this fall, and I realized that the best way to teach this stuff
is to be able to play with the representations directly, in particular
to be able to see it in action on a simpler system than full 64-bit
precision, especially when str(f) or repr(f) won't show *all* of the
significant digits stored in a float. The decimal class deliberately
avoids binary representation issues, and I can't find what I want
online.

Consequently, I have written a module to simulate the machine
representation of binary floating point numbers and their arithmetic.
Values can be of arbitrary fixed precision or infinite precision,
along the same lines as python's in-built decimal class. The code is
here:  http://www2.gsu.edu/~matrhc/binary.html

The design is loosely based on that decimal module, although it
doesn't get in to threads, for instance. You can play with different
IEEE 754 representations with different precisions and rounding modes,
and compare with infinite precision Binary numbers. For instance, it
is easy to learn about machine epsilon, representation/rounding error
using a much simpler format such as a 4-bit exponent and 6-bit
mantissa. Such a format is easily defined in the new module and can be
manipulated easily:

>>> context = define_context(4, 6, ROUND_DOWN)

>>> zero = context(0)
Binary("0", (4, 6, ROUND_DOWN))

>>> print zero      # sign, characteristic, significand bits
0 0000 000000

>>>  zero.next()
Binary("0.001E-9", (4, 6, ROUND_DOWN))

>>> print zero.next()
0 0000 000001

>>> largest_denormalized = context('0 0000 111111')  # direct spec of the sign, characteristic, and significand bits
>>> largest_denormalized
Binary("0.111111E-6", (4, 6, ROUND_DOWN))

>>> largest_denormalized.as_decimal()
Decimal("0.015380859375")

>>> n01 = context(0.1)  # nearest representable is actually stored

>>> print n01, "    rounded to ", n01.as_decimal()
0 0011 100110      rounded to 0.099609375

>>> Binary('-10111.0000001').as_decimal()
Decimal("-23.0078125")

>>> Binary('-10111.0000001', context).as_decimal()     # not enough precision in this context
Decimal("-23.0000")

>>> diff = abs(Binary('-10111.0000001') - Binary('-10111.0000001', context))
>>> diff
Binary("0.1E-6", (4, 6, ROUND_DOWN))

The usual arithmetic operations are permitted on these objects, as
well as representations of their values in decimal or binary form.
Default contexts for half, single, double, and quadruple IEEE 754
precision floats are provided. Binary integer classes are also
provided, and some other utility functions for converting between
decimal and binary string representations. The module is compatible
with the numpy float classes and requires numpy to be installed.

The source code is released under the BSD license, but I am amenable
to other licensing ideas if there is interest in adapting the code for
some other purpose. Full details of the functionality and known issues
are in the module's docstring, and many examples of usage are in the
accompanying file binary_tests.py (which also acts to validate the
common representations against the built-in floating point types). I
look forward to hearing feedback, especially in case of bugs or
suggestions for improvements.

-Rob

-- 
Robert H. Clewley, Ph. D.
Assistant Professor
Department of Mathematics and Statistics
Georgia State University
720 COE, 30 Pryor St
Atlanta, GA 30303, USA

tel: 404-413-6420 fax: 404-413-6403
http://www2.gsu.edu/~matrhc
http://brainsbehavior.gsu.edu/