"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 every finite floating-point number y we expect computed values x := ( y - ulp(y) rounded ) and z := ( y + ulp(y) rounded ) to satisfy . We cannot expect in every case; consider the case y := 1 and remember that by default IEEE 754 rounds midway cases to nearest even so z = 1 . These inferences hold for decimal floating-point too though for very different reasons. By defining ulp(Infinity) to be finite we sustain these inferences at infinite y too; they would be invalidated there otherwise. Tighter inferences are available just when y is finite; then

or both

Both strict inequalities are false if y is infinite; even so, no INVALID OPERATON occurs. Were ulp(Infinity) = Infinity , an INVALID OPERATION would occur at x := y - ulp(y) producing NaN, and again at each comparison of x and y . Gratuitous INVALID OPERATIONs like these seem to punish a program cruelly and unusually: they would wrest control from it in many environmemts, though not MATLAB's.

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 !