2832-Principles:.sagews

Lab 1: Representing real numbers

Principles:

• aid the students in doing hand computations, 

• check their answers,

• visualize their answers,

• help them generalize to world-size problems

(1.0).nextabove().str(truncate=False)

'1.0000000000000002'

x=1.0

x.nextabove().str(truncate=False)

'1.0000000000000002'

x.str(base=2)

'1.0000000000000000000000000000000000000000000000000000'

x.nextabove().str(base=2)

'1.0000000000000000000000000000000000000000000000000001'

x.nextbelow().str(truncate=False)

'0.99999999999999989'

x.nextbelow().str(base=2)

'0.11111111111111111111111111111111111111111111111111111'

x.nextabove()-x

2.22044604925031e-16

This means that our number $x$ above actually represented a whole range of numbers.  The range of numbers between 1.0 and the next number above it is related to the machine epsilon, which is used in a number of contexts to understand how precise a machine can be and estimate roundoff errors.

A finite system

Let's imagine we have a 2-bit mantissa and a 2-bit exponent.  Then all of our machine numbers look like $1.b_1b_2\times 2^{k-2}$, where $b_1$ and $b_2$ can each be $0$ or $1$, and $k$ can be $0,1,2,3$.

Question: List all possible numbers in our set.

Answer: (double-click to edit)

nums=[RR('1.%d%d'%(b1,b2),base=2)*2^(k-2) for b1,b2,k in CartesianProduct([0,1],[0,1],[0,1,2,3])]
sorted(nums)

[0.250000000000000, 0.312500000000000, 0.375000000000000, 0.437500000000000, 0.500000000000000, 0.625000000000000, 0.750000000000000, 0.875000000000000, 1.00000000000000, 1.25000000000000, 1.50000000000000, 1.75000000000000, 2.00000000000000, 2.50000000000000, 3.00000000000000, 3.50000000000000]

or as fractions

sorted([i.exact_rational() for i in nums])

[1/4, 5/16, 3/8, 7/16, 1/2, 5/8, 3/4, 7/8, 1, 5/4, 3/2, 7/4, 2, 5/2, 3, 7/2]

The distribution of these numbers is interesting.  Let's plot them on the real line.

points(zip(nums,[0]*len(nums)),pointsize=25).show(xmin=0,xmax=4,figsize=[10,1])

Note that the numbers are unevenly spaced.  In fact, count how many numbers are in $[2^{-2},2^{-1})$, and how many are in $[2^{-1},2^0)$, and between $[2^0,2^1)$.  What do you notice?  Is there a pattern, and why?

Answer:



Note that there is no machine number representing 1.1.  In order to deal with machine arithmetic, the machine needs to decide how to represent 1.1.  We notate the machine's version of 1.1 as $fl(1.1)$ ($fl$ stands for "float", as in "floating point approximation").

What are the reasonable possibilities for $fl(1.1)$?  What would you choose, and why?

Answer:

The gap around zero is called the "hole at zero".  Where does this come from? (Subnormal numbers fill this gap in IEEE 754-2008.)

Answer:

What happens if we do $5/16-1/4=1/16$ on the machine? In other words, what is $fl(fl(5/16)-fl(1/4))$?

Answer:

What is $fl(fl(1)+fl(5/4))$?

Answer:

What is $1.1+1.1+1.1$, as the machine does it?  (first write down all applicable $fl$, then calculate the result.)

Answer:

 

Rounding Modes

Look up in Wikipedia what the possible rounding modes are for IEEE 754-2008.  For each of these rounding modes, what is $fl(3.25)$? (To answer this for some of the rounding modes, you'll need to know what the mantissas for 3 and 3.5 are.)

Answer:

Relative Error

What is the relative error of the approximation $fl(1.1)$ and 1.1? (See the definition from last class or at the bottom of p. 19 of the text.)

Answer:

What is the maximum relative error for a number $x$ such that $fl(x)=1$?  This relative error is called the machine epsilon, and is often denoted with the Greek epsilon $\epsilon$ or $\epsilon_{\textrm{mach}}$. 

Answer:

What is the maximum relative error for a number $y$ such that $fl(y)=2$? 

Answer:

What about the maximum relative error for a number $z$ such that $fl(z)=1/2$?

Answer:

Do you notice a pattern? What about the maximum relative error for a number $w$ such that $fl(z)=3$?

Answer:

What is the machine epsilon for a 53-bit number? (Work it out by hand, from what you know about the exponent range, and verify the next computation.)

((RR(1).nextabove()-RR(1))/2*RR(1))

1.1102230246251565e-16

So intuitively, what does it mean if you do several computations using 53-bit precision floating point numbers (i.e., double precision) and the results differ by about $10^{-16}$, or if you get an answer that is about $10^{-16}$?

2.0^(-53) # Here's a hint for the by-hand part above...

1.1102230246251565e-16

You can read about two disasters in history coming from ignoring numerical issues here: http://www.ima.umn.edu/~arnold/455.f96/disasters.html