"how numbers are stored and used in computers"

Units in last place

The measurement ulp stands for Unit in the Last Place, and is a common measurement of mathematical accuracy. Kahan's original definition of ulp in the 1960s was

is the gap between the two floating-point numbers nearest , even if is one of them.

That is, if you take the binary representation of a floating point number and increment it, the difference is 1 ulp.

"Ulp" was coined to characterize the accuracies of implementations of run-time math library functions on the IBM 7090 and 7094, which were not as accurate nor fast as they should have been in 1960. By 1964 most of their deficiencies had been repaired by contributions to SHARE (the organized users of IBM mainframes), and especially by the work of Hirondo Kuki at the University of Chicago.

Kahan, 2004
From time to time, well-intentioned attempts to speed up a math library introduced inaccuries so large that Fred Gustavson at IBM Yorktown Heights coined the word "GULP" (Giga-Ulp) to describe them. Today, there is no excuse for a math library less accurate than the Freely Distributed Math Library fdlibm, promulgated from Sun Microsystems mostly by graduates from the University of California at Berkeley, who created the libraries distributed with 4.3 BSD Berkeley UNIX in the 1980s.

By then the adoption of IEEE Standard 754 had made infinities and NaNs so ubiquitous that the definition of an ulp had to be changed:

is the gap between the two finite floating-point numbers nearest , even if is one of them, and is .

E.g., for IEEE 754 Double Precision (8 bytes wide) as in MATLAB, ulp(1) = eps/2, ulp(1.5) = eps, ulp(inf) = 2^971, ulp(0) = 2^-1074 .

Why finite?

Without this word we would have ulp(Infinity) = Infinity , which is not too bad but not so good in situations like the following:

For any finite floating-point number , we expect the computed values and to satisfy the inequality . However, we cannot always expect the strict inequality to hold. For instance, when , the IEEE 754 standard rounds midway cases to the nearest even number, resulting in . These inferences are also applicable to decimal floating-point numbers, albeit for different reasons.

By defining as finite, we maintain these inferences even when is infinite; otherwise, they would be invalidated. When is finite, tighter inferences are possible: either , , or both.

If is infinite, both strict inequalities are false, yet no invalid operation occurs. If were infinite, subtracting it from would result in an invalid operation, producing , and further invalid comparisons between and would occur. Such gratuitous invalid operations could disrupt a program significantly, except in environments like MATLAB where this is handled gracefully.

Here is a vectorized MATLAB program to compute ulp(x) on machines whose arithmetics conform to IEEE Standard 754, as almost all do now.

code.txt
1function  u = ulp(x)
2%  ulp(x)  is the gap between the two finite floating-point numbers
3%  nearest  x ,  even if  x  is one of them.  ulp(NaN)  is  NaN .  E.g.,
4%  ulp(1) = eps/2,  ulp(1.5) = eps,  ulp(inf) = 2^971,  ulp(0) = 2^-1074
5
6u = abs(real(x)) ;                  %... to propagate  NaN
7k = ( u < 5.0e-308 ) ;              %... copes with  x  too near  0
8if any(k(:)),  v = k*5.0e-324 ;  u(k) = v(k) ;  end
9j = ( 9.0e307 < u ) ;               %... copes with  x  near  Infinity
10if any(j(:)),  v = j*(2^(1024-53)) ;  u(j) = v(j) ;  end
11j = ~(j|k) ;                        %... copes with the rest
12if any(j(:))
13    v = u(j) ;  w = (0.7*eps)*v ;   %... "0.7"  works for  IEEE 754
14    u(j) = min( (v+w)-v, (w-v)+v ) ;  end
15if any(imag(x(:))),  u = u + i*ulp(imag(x)) ;  end % recursive call !