Emulating fixed precision in python

Question

For a university course in numerical analysis we are transitioning from Maple to a combination of Numpy and Sympy for various illustrations of the course material. This is because the students already learn Python the semester before.

One of the difficulties we have is in emulating fixed precision in Python. Maple allows the user to specify a decimal precision (say 10 or 20 digits) and from then on every calculation is made with that precision so you can see the effect of rounding errors. In Python we tried some ways to achieve this:

• Sympy has a rounding function to a specified number of digits.
• Mpmath supports custom precision.

This is however not what we're looking for. These options calculate the exact result and round the exact result to the specified number of digits. We are looking for a solution that does every intermediate calculation in the specified precision. Something that can show, for example, the rounding errors that can happen when dividing two very small numbers.

The best solution so far seems to be the custom data types in Numpy. Using float16, float32 and float64 we were able to al least give an indication of what could go wrong. The problem here is that we always need to use arrays of one element and that we are limited to these three data types.

Does anything more flexible exist for our purpose? Or is the very thing we're looking for hidden somewhere in the mpmath documentation? Of course there are workarounds by wrapping every element of a calculation in a rounding function but this obscures the code to the students.

Answer 1

You can use decimal. There are several ways of usage, for example, localcontext or getcontext.

Example with getcontext from documentation:

>>> from decimal import *
>>> getcontext().prec = 6
>>> Decimal(1) / Decimal(7)
Decimal('0.142857')

Example of localcontext usage:

>>> from decimal import Decimal, localcontext
>>> with localcontext() as ctx:
...     ctx.prec = 4
...     print Decimal(1) / Decimal(3)
... 
0.3333

To reduce typing you can abbreviate the constructor (example from documentation):

>>> D = decimal.Decimal
>>> D('1.23') + D('3.45')
Decimal('4.68')