9. Using Floating-Point Numbers

To say that the implementation of numerical algorithms that exhibit the same behavior across a wide variety of platforms is difficult, is an understatement. This section provides very little help but we hope it is worth reading. Any additional suggestions and information are very much appreciated as we would like to expand this section.

9.1. Machine Constants

One problem when writing numerical algorithms is obtaining machine constants. Typical values one needs are:

• The radix of the floating-point representation.

• The number of digits in the floating-point significand expressed in terms of the radix of the representation.

• The number of bits reserved for the representation of the exponent.

• The smallest positive floating-point number epsilon such that 1.0 + epsilon != 1.0.

• The smallest non-vanishing normalized floating-point power of the radix.

• The largest finite (12) floating-point number.

On Suns, they can be obtained in <values.h>. The ANSI C Standard recommends that such constants be defined in the header file <float.h>.

Suns and standards apart, these values are not always readily available, e.g., in Tektronix workstations running UTek. One solution is to use a modified version of a program that can be obtained from the network which is called machar. Machar is described in Algorithm 665, MACHAR: A Subroutine to Dynamically Determine Machine Parameters [2] and can obtained by anonymous FTP from the netlib. (13)

It is straightforward to modify the C version of machar to generate a C preprocessor file that can be included directly by C programs.

There is also a publicly available program called config.c that attempts to determine many properties of the C compiler and machine that it is run on. It can generate the ANSI C header files <float.h> and <limits.h> among other useful features. This program was submitted to comp.sources.misc. (14) The latest version, 4.2, is available by FTP from mcsun.eu.net in directory misc and is called config42.c (the next version, 4.3, will be called enquire.c). Version 4.2 is also distributed with gcc, where it is called hard-params.c.

9.2. Floating-Point Arguments

In the days of K&R [9] one was "encouraged" to use float and double interchangeably (15) since all expressions with such data types where always evaluated using the double representation -- a real nightmare for those implementing efficient numerical algorithms in C. This rule applied, in particular, to floating-point arguments and for most compilers around, it does not matter whether one defines the argument as float or double.

According to the ANSI C Standard, such programs will continue to exhibit the same behavior as long as one does not prototype the function. Therefore, when prototyping functions, make sure that the prototype is included when the function definition is compiled so the compiler can check if the arguments match.

9.3. Floating-Point Arithmetic

Be careful when using the == and != operators to compare floating-point types. Expressions such as

	if (float_expr1 == float_expr2)

	10^50 + 812 - 10^50 + 10^55 + 511 - 10^55 = 812 + 511 = 1323

Most computers will produce zero regardless of whether one uses float or double. Although the absolute error is large, the relative error is quite small and probably acceptable for many applications.

It is rather better to use expressions such as

	|float_expr1 - float_expr2| <= K

	||float_expr1/float_expr2| - 1.0| <= K

1. The floating type, e.g., float or double,

2. the machine architecture (the machine constants defined in the previous section), and

3. the precision of the input values and the rounding errors introduced by the numerical method used.

Other possibilities exist and the choice depends on the application.

The development of reliable and robust numerical algorithms is a very difficult undertaking. Methods for certifying that the results are correct within reasonable bounds must usually be implemented. A reference such as NUMERICAL RECIPES in C [12] is always useful.

• Keep in mind that the double representation does not necessarily increase the precision. Actually, in some implementations the precision decreases, but the range increases.

• Do not use double unnecessarily, since in many cases there is a large performance penalty. Furthermore, there is no point in using higher precision, if the additional bits that would be computed are garbage anyway. The precision one needs depends mostly on the precision of the input data and the numerical method used.

9.4. Exceptions

Floating-point exceptions (overflow, underflow, division by zero, etc) are not signaled automatically in some systems. In that case, they must be explicitly enabled.

Always enable floating-point exceptions, since they may be an indication that the method is unstable. Otherwise, one must be sure that such events do not affect the output.

12. Some representations have reserved values for +inf and -inf.

13. Email (Internet) address is netlib@ornl.gov. For more information, send a message containing the line send index to that address.

14. The archive site of comp.sources.misc is uunet.uu.net.

15. In fact one wonders why they even bothered to define two representations for floating-point numbers considering the rules applied to them.

contents

Run-time Library

VMS