Section 14. Floating Point

14.1:   When I set a float variable to, say, 3.1, why is printf printing
	it as 3.0999999?

A:      Most computers use base 2 for floating-point numbers as well as
	for integers.  In base 2, one divided by ten is an infinitely-
	repeating fraction (0.0001100110011...), so fractions such as
	3.1 (which look like they can be exactly represented in decimal)
	cannot be represented exactly in binary.  Depending on how
	carefully your compiler's binary/decimal conversion routines
	(such as those used by printf) have been written, you may see
	discrepancies when numbers (especially low-precision floats) not
	exactly representable in base 2 are assigned or read in and then
	printed (i.e. converted from base 10 to base 2 and back again).
	See also question 14.6.

14.2:   I'm trying to take some square roots, but I'm getting crazy
	numbers.

A:      Make sure that you have #included <math.h>, and correctly
	declared other functions returning double.  (Another library
	function to be careful with is atof(), which is declared in
	<stdlib.h>.)  See also question 14.3 below.

	References: CT&P Sec. 4.5 pp. 65-6.

14.3:   I'm trying to do some simple trig, and I am #including <math.h>,
	but I keep getting "undefined: sin" compilation errors.

A:      Make sure you're actually linking with the math library.  For
	instance, under Unix, you usually need to use the -lm option, at
	the *end* of the command line, when compiling/linking.  See also
	questions 13.25, 13.26, and 14.2.

14.4:   My floating-point calculations are acting strangely and giving
	me different answers on different machines.

A:      First, see question 14.2 above.

	If the problem isn't that simple, recall that digital computers
	usually use floating-point formats which provide a close but by
	no means exact simulation of real number arithmetic.  Underflow,
	cumulative precision loss, and other anomalies are often
	troublesome.

	Don't assume that floating-point results will be exact, and
	especially don't assume that floating-point values can be
	compared for equality.  (Don't throw haphazard "fuzz factors"
	in, either; see question 14.5.)

	These problems are no worse for C than they are for any other
	computer language.  Certain aspects of floating-point are
	usually defined as "however the processor does them" (see also
	question 11.34), otherwise a compiler for a machine without the
	"right" model would have to do prohibitively expensive
	emulations.

	This article cannot begin to list the pitfalls associated with,
	and workarounds appropriate for, floating-point work.  A good
	numerical programming text should cover the basics; see also the
	references below.

	References: Kernighan and Plauger, _The Elements of Programming
	Style_ Sec. 6 pp. 115-8; Knuth, Volume 2 chapter 4; David
	Goldberg, "What Every Computer Scientist Should Know about
	Floating-Point Arithmetic".

14.5:   What's a good way to check for "close enough" floating-point
	equality?

A:      Since the absolute accuracy of floating point values varies, by
	definition, with their magnitude, the best way of comparing two
	floating point values is to use an accuracy threshold which is
	relative to the magnitude of the numbers being compared.  Rather
	than

		double a, b;
		...
		if(a == b)      /* WRONG */

	use something like

		#include <math.h>

		if(fabs(a - b) <= epsilon * fabs(a))

	for some suitably-chosen degree of closeness epsilon (as long as
	a is nonzero!).

	References: Knuth Sec. 4.2.2 pp. 217-8.

14.6:   How do I round numbers?

A:      The simplest and most straightforward way is with code like

		(int)(x + 0.5)

	This technique won't work properly for negative numbers,
	though (for which you could use something like
	(int)(x < 0 ? x - 0.5 : x + 0.5)).

14.7:   Why doesn't C have an exponentiation operator?

A:      Because few processors have an exponentiation instruction.
	C has a pow() function, declared in <math.h>, although explicit
	multiplication is usually better for small positive integral
	exponents.

	References: ISO Sec. 7.5.5.1; H&S Sec. 17.6 p. 393.

14.8:   The predefined constant M_PI seems to be missing from my
	machine's copy of <math.h>.

A:      That constant (which is apparently supposed to be the value of
	pi, accurate to the machine's precision), is not standard.  If
	you need pi, you'll have to define it yourself, or compute it
	with 4*atan(1.0).

	References: PCS Sec. 13 p. 237.

14.9:   How do I test for IEEE NaN and other special values?

A:      Many systems with high-quality IEEE floating-point
	implementations provide facilities (e.g. predefined constants,
	and functions like isnan(), either as nonstandard extensions in
	<math.h> or perhaps in <ieee.h> or <nan.h>) to deal with these
	values cleanly, and work is being done to formally standardize
	such facilities.  A crude but usually effective test for NaN is
	exemplified by

		#define isnan(x) ((x) != (x))

	although non-IEEE-aware compilers may optimize the test away.

	C9X will provide isnan(), fpclassify(), and several other
	classification routines.

	Another possibility is to to format the value in question using
	sprintf(): on many systems it generates strings like "NaN" and
	"Inf" which you could compare for in a pinch.

	See also question 19.39.

	References: C9X Sec. 7.7.3.

14.11:  What's a good way to implement complex numbers in C?

A:      It is straightforward to define a simple structure and some
	arithmetic functions to manipulate them.  C9X will support
	complex as a standard type.  See also questions 2.7, 2.10, and
	14.12.

	References: C9X Sec. 6.1.2.5, Sec. 7.8.

14.12:  I'm looking for some code to do:
		Fast Fourier Transforms (FFT's)
		matrix arithmetic (multiplication, inversion, etc.)
		complex arithmetic

A:      Ajay Shah has prepared a nice index of free numerical
	software which has been archived pretty widely; one URL
	is ftp://ftp.math.psu.edu/pub/FAQ/numcomp-free-c .
	See also questions 18.13, 18.15c, and 18.16.

14.13:  I'm having trouble with a Turbo C program which crashes and says
	something like "floating point formats not linked."

A:      Some compilers for small machines, including Borland's
	(and Ritchie's original PDP-11 compiler), leave out certain
	floating point support if it looks like it will not be needed.
	In particular, the non-floating-point versions of printf()
	and scanf() save space by not including code to handle %e, %f,
	and %g.  It happens that Borland's heuristics for determining
	whether the program uses floating point are insufficient,
	and the programmer must sometimes insert a dummy call to a
	floating-point library function (such as sqrt(); any will
	do) to force loading of floating-point support.  (See the
	comp.os.msdos.programmer FAQ list for more information.)