Controlling FPU rounding modes with Python
Hey folks! In the previous post I talked a little about the theory behind the floating-point arithmetic and rounding methods. This time I’ll try a more practical approach showing how to control the rounding modes that FPUs (Floating-Point Units) employ to perform their operations. For that I’ll use Python. But, first let me introduce the <fenv.h> header of the ANSI C99 standard.
The ANSI C99 standard establishes the <fenv.h> header to control the floating-point environment of the machine. A lot of functions and constants allowing to configure the properties of the floating-point system (including rounding modes) are defined in this header. In our case, the important definitions are:
Functions
int fegetround()
– Returns the current rounding modeint fesetround(int mode)
– Sets the rounding mode returning 0 if all went ok and some other value otherwise
Constants
FE_TOWARDZERO
– Flag for round to zeroFE_DOWNWARD
– Flag for round toward minus infinityFE_UPWARD
– Flag for round toward plus infinityFE_TONEAREST
– Flag for round to nearest
Unfortunately, not all compilers fully implement the ANSI C99 standard. Therefore, there is no guarantee on the portability of code that uses the <fenv.h> header. Considering this caveat, the GCC compiler supports it, which makes it possible to build a Python extension wrapping the features of <fenv.h>. However, this is not a smart way to reach our goal since it would require forcing users to compile the source code of the extension. An alternative is to use the ctypes Python module to instantiate the libm (part of the standard C library used by GCC, typically glibc, that implements <fenv.h>):
from ctypes import cdll from ctypes.util import find_library libm = cdll.LoadLibrary(find_library('m'))
The constants that identify the rounding modes are usually defined as macros in <fenv.h> and, therefore, are not accessible via ctypes. We need to redefine them in Python. The problem is that they vary according the processor so that it is necessary to establish some logic on the result of the function platform.processor()
. Right now, however, I just have a x86 processor to test, so I will use the constants for it only:
FE_TOWARDZERO = 0xc00 FE_DOWNWARD = 0x400 FE_UPWARD = 0x800 FE_TONEAREST = 0
Now, just call the appropriate functions:
>>> back_rounding_mode = libm.fegetround() >>> libm.fesetround(FE_DOWNWARD) 0 >>> 1.0/10.0 0.099999999999999992 >>> libm.fesetround(FE_UPWARD) 0 >>> 1.0/10.0 0.10000000000000001 >>> libm.fesetround(back_rounding_mode)
But, what about the MS Visual Studio? Well… the standard C library of this compiler, msvcrt (MS Visual Studio C Run-Time Library), does not support <fenv.h>. Nonetheless, MS Visual Studio implements a set of non-standard (quite typical…) constants and functions specialized in manipulating the properties of the machine FPU. These constants and functions are defined in the <float.h> header distributed with that compiler. The following definitions are of particular interest to us:
Functions
unsigned int _controlfp(unsigned int new, unsigned int mask)
– Sets/gets the control vector of the floating-point system (more information here)
Constants
_MCW_RC = 0x300
– Control vector mask for information about rounding modes_RC_CHOP = 0x300
– Control vector value for round to zero mode_RC_UP = 0x200
– Control vector value for round toward plus infinity mode_RC_DOWN = 0x100
– Control vector value for round toward minus infinity mode_RC_NEAR = 0
– Control vector value for round to nearest mode
Analogous to the previous case:
>>> _MCW_RC = 0x300 >>> _RC_UP = 0x200 >>> _RC_DOWN = 0x100 >>> from ctypes import cdll >>> msvcrt = cdll.msvcrt >>> back_rounding_mode = msvcrt._controlfp(0, 0) >>> msvcrt._controlfp(_RC_DOWN, _MCW_RC) 590111 >>> 1.0/10.0 0.099999999999999992 >>> msvcrt._controlfp(_RC_UP, _MCW_RC) 590367 >>> 1.0/10.0 0.10000000000000001 >>> msvcrt._controlfp(back_rounding_mode, _MCW_RC) 589855
For sure, it is not convenient to type all this every time you want to change the rounding mode. Nor it is appropriate to guess the standard C library of the system. But you can always create functions encapsulating this behavior. Something like this:
def _start_libm(): global TO_ZERO, TOWARD_MINUS_INF, TOWARD_PLUS_INF, TO_NEAREST global set_rounding, get_rounding from ctypes import cdll from ctypes.util import find_library libm = cdll.LoadLibrary(find_library('m')) set_rounding, get_rounding = libm.fesetround, libm.fegetround # x86 TO_ZERO = 0xc00 TOWARD_MINUS_INF = 0x400 TOWARD_PLUS_INF = 0x800 TO_NEAREST = 0 def _start_msvcrt(): global TO_ZERO, TOWARD_MINUS_INF, TOWARD_PLUS_INF, TO_NEAREST global set_rounding, get_rounding from ctypes import cdll msvcrt = cdll.msvcrt set_rounding = lambda mode: msvcrt._controlfp(mode, 0x300) get_rounding = lambda: msvcrt._controlfp(0, 0) TO_ZERO = 0x300 TOWARD_MINUS_INF = 0x100 TOWARD_PLUS_INF = 0x200 TO_NEAREST = 0 for _start_rounding in _start_libm, _start_msvcrt: try: _start_rounding() break except: pass else: print "ERROR: You couldn't start the FPU module"
In this case, the constants and functions to be called, as well as the standard C library instantiated, are abstracted. Just use the constants TO_ZERO
, TOWARD_MINUS_INF
, TOWARD_PLUS_INF
, TO_NEAREST
and the functions set_rounding
, get_rounding
.
So, that’s it… I hope this has been useful to you somehow (I doubt it… :P) or, at least, that it was interesting… Until next time…
Posted on March 30, 2009, in Floating-Point Numbers, Python, Scientific Computing and tagged floating-point, fpu, python, rounding. Bookmark the permalink. 2 Comments.
Excelente post!
mas como eu ainda to engatinhando em python, um post mais basico seria de maior ajuda para mim.
abs
oi rafa, ta massa teu blog! parabéns pela iniciativa. continua postando, principalmente sobre python, e outras coisas legais tb.