Fixed precision emulation in python

Question

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

One of the difficulties we have is to emulate fixed precision in Python. Maple allows the user to specify decimal precision (for example, 10 or 20 digits), and from now on, each calculation is done with such precision that you can see the effect of rounding errors. In Python, we tried to do several ways:

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

This, however, is not what we are looking for. These parameters calculate the exact result and round the exact result to the specified number of digits. We are looking for a solution that performs each intermediate calculation with a given accuracy. Something that can show, for example, rounding errors that can occur when dividing two very small numbers.

The best solution seems to be a custom data type in Numpy. Using float16, float32, and float64, we could at least give an indication of what might 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.

Is there something more flexible for our purpose? Or is this what we are looking for hidden somewhere in the mpmath documentation? Of course, workarounds wrap each element of the calculation in a rounding function, but this hides the code for students.

Answer 1

You can use decimal . There are several uses, for example localcontext or getcontext .

Example with getcontext from the documentation:

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

An example of using localcontext :

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

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

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