Return to Python Lectures

When arithmetic operations are performed on operands of differing types, including Numeric arrays, Python promotes all types to the highest and stores the result as the resulting highest type. An operation on a Float32 and a Float (which is double precision) results in a Float. For operations involving arrays alone, the behavior is usually what one would anticipate, but some surprising behavior results from the fact that Python floating point variables are always double precision. There is no single precision Python type.

After adding 1. to the array x, the result z is now a default float, which is double precision. Similarly, when we increment x by 1, x has been changed to double precision. Note that assignment of Numeric arrays just changes a reference, rather than copying data, the statement "x=x+1." actually creates a new double precision array, and then discards the reference to the old array (and also the storage, if nothing else refers to it), redefining x as the reference to the new double precision array.

The promotion to double is a kind of social engineering in Python, encouraging people to think twice (or more) before deciding to do single-precision floating arithmetic. However, this often creates undesired behavior, especially is one is interacting with imported data or compiled routines that produce or expect single precision. Consider the following simple array cross section assignment, which often occurs in the course of accumulating averages into some slice of an array:

>>> T = Numeric.ones((2,5),Numeric.Float32)
>>> T[0] = T[0]/12.
Traceback (most recent call last):
File "<pyshell#80>", line 1, in ?
T[0] = T[0]/12.
TypeError: Array can not be safely cast to required type


This innocuous assignment produces a cryptic error message. What is happening is that the right hand side produces a double precision array, which one is trying to stuff into a single precision array cross-section. Some of us think that the desired behavior for Numeric in a case like this would be to downcast the result to single precision before copying back into the array cross section, but Numeric doesn't work that way. Fortunately, there are some workarounds that are not too cumbersome.

One approach is to simply do everything in double precision, and then downcast to single precision before invoking a compiled routine (e.g. a time stepper) that expects single precision. In Numeric, one can convert types using the astype(...) method of an array. Consider the following examples:

>>> z = Numeric.zeros(5,Numeric.Float)
>>> z
array([ 0., 0., 0., 0., 0.])
>>> z = z.astype(Numeric.Float32)
>>> z
array([ 0., 0., 0., 0., 0.],'f')
>>> z = z.astype(Numeric.Complex)
>>> z
array([ 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j])

Since a type cast creates new storage and does a copy, this may cause some inefficiencies if it has to be done too often. An alternative is to set the "savespace" flag of an array. This flag tells Numeric not to promote when doing operations. The savespace flag setting is inherited by the result of the operation, as shown in the following example:

>>> T = Numeric.ones((2,5),Numeric.Float32)
>>> T.savespace(1)
>>> T[0] = T[0]/12.
>>> T[0] = T[0]*12.
>>> T[0]
array([ 1., 1., 1., 1., 1.],'f')


Now, T has remained single precision, and the array cross section assignment works. Note, however, that in operations involving two arrays, Numeric has no choice but to promote to the higher type, even if the savespace flag is set:

>>> x = Numeric.ones(5,Numeric.Float)
>>> y = Numeric.ones(5,Numeric.Float32)
>>> x.savespace(1)
>>> y.savespace(1)
>>> z = x+y
>>> z
array([ 2., 2., 2., 2., 2.])

In a case like this, if one really wanted z to be single precision, one would have to cast x to single precision using astype before the operation. This can be done in-line using the following, of one wants to retain the original double precision array:

>>> z = x.astype(Numeric.Float32) + y
>>> z
array([ 2., 2., 2., 2., 2.],'f')

Somewhat cumbersome, to be sure, but at least one is forced to be aware of what one is doing whenever there is a loss of precision. The general rule of thumb is that Python will force you to explicitly give permission whenever there is a risk of loss of precision in mixed operations.

Return to Python Lectures