\label{se*se}

In scientific computing, we never expect to get the exact answer. 
Inexactness is practically the definition of scientific computing.  
Getting the exact answer, generally with integers or rational numbers, is 
{\em symbolic computing}, an interesting but distinct subject.

To put this technically, suppose we are trying to compute the number
$A$.  The computer will produce an approximation, which we call
$\widehat{A}$.  This notation is borrowed from statistics, where an 
estimate for a quantity $A$ is often called $\widehat{A}$.
Whether $A$ and $\widehat{A}$ agree to $16$ decimal
places or differ by $30\%$, we never expect the identity
$A = \widehat{A}$ to be true in the mathematical sense, if only because
we do not have an exact representation for $A$.  For
example, if we need to find $x$ so that $x^2 = 175$, we might
get $13$ or $13.22876$, depending on our computational method, but we 
cannot represent the exact $\sqrt{175}$ if we only allow finitely many 
digits.   

Since the answer is never exactly right, it can be hard to know whether
errors in computed results are due to a poor algorithm, an impossible problem
(technically, {\em ill conditioned}), or a software bug.  Understanding the
sources of error will help us know their behavior well enough to recognize
them and tell them apart.  We will use this constantly in devising algorithms
and building and testing software.

Of the four sources of error mentioned, we discuss only roundoff error in 
detail here.  The others are discussed at length in later chapters.
To explain roundoff error, we pause to discuss computer floating point 
arithmetic in some detail.  Thanks to the IEEE Floating Point 
Standard [], it we can predict the outcome of inexact computer arithmetic
in most cases on most computers.

\section{Relative and absolute error} \label{RelAbsErr*se}

The {\em absolute error} in approximating $A$ by $\widehat{A}$  is 
$e = \widehat{A} - A$. The {\em relative error}, which is $\epsilon = e/A$, 
is usaually more meaningful.  Some algebra puts these definitions into the
form
\begin{equation}
\widehat{A} = A + e 
\;\;\;\mbox{(absolute error)}\;\;,\;\;\;\;\;\;\;\;
\widehat{A} = A \cdot (1 + \epsilon) 
\;\;\;\mbox{(relative error).}
\label{error*se}    \end{equation}
For example, the
absolute error in approximating $A=\sqrt{175}$ by $\widehat{A}=13$ is
$e=13.22876\cdots - 13 \approx .229$.  The corresponding relative error
is $e/A \approx .229/13.23 \approx .017 < 2\%$.  Saying that the error
is less than $2\%$ is probably more informative than saying that the error
is less than $.25=1/4$.  

Relative error is a {\em dimensionless} measure of error.  In practical 
situations, the desired $A$ probably has units, seconds,
meters, etc.  Knowing that the error is, say, $.229$ meters does not tell 
you whether that is large or small.  If the correct length is half a meter,
then $.229$ meters is a large error.  If the correct length is 
$13.22876\cdots$
meters, the approximate answer is off by less than $2\%$.  If we measure
lengths in centimeters, then the error becomes $22.9\mbox{cm}$.  Is $22.9$
a large error?  It is less than $2\%$ of the exact length, 
$1,322.876\cdots\mbox{cm}$.

\section{Computer arithmetic}     \label{CompArith*se}
One of the many sources of error in scientific computing is inexact
computer arithmetic, which is called {\em roundoff error}.  Roundoff
error is inevitable but its consequences vary.  Some computations 
yield nearly exact results, while others give answers that are 
completely wrong.  This is true even for problems that involve no
other approximations.  For example, solving systems of linear equations 
using gaussian elimination would give the exact answer
if all the computations were performed exactly.  When these computations
are done in finite precision floating point arithmetic, we might or
might not get something close to the right answer.  

A single floating point operation almost
always produces high relative accuracy.  For given $B$ and $C$, let 
$A=B\bigcirc C$, with $\bigcirc$ standing for one of the arithmetic operations:
addition, subtraction, multiplication, or division.  Then, {\em with the
same $B$ and $C$}, the computer will produce $\widehat{A}$ which
satisfies (\ref{error*se}) with 
$\left|\epsilon\right| \leq \epsilon_{\mbox{mach}}$, where 
$\epsilon_{\mbox{mach}}$ is the {\em machine precision}.  Normally
$\epsilon_{\mbox{mach}}$ is $2^{-24} \approx 6\cdot 10^{-8}$ for single
precision and $2^{-53} \approx 10^{-16}$ for double precision.  The 
rough approximation $2^{10} = 1024 \approx 1000 = 10^3$ gives
$2^{24} = 2^4 \cdot 2^{20} = 16 \cdot (2^{10})^2 
\approx 16 \cdot (10^3)^2 = 16 \cdot 10^6$, so 
$2^{-24} \approx \frac{1}{16} 10^{-6} = 6.25\cdot 10^{-8}$.  Also,
$2^{53} \approx 8 \cdot 10^{15}$ so the double precision machine precision
is about $\frac{1}{8} \cdot 10^{-15} = 1.25\cdot 10^{-16}$.

With uniformly high relative accuracy, how does computer arithmetic ever lead
to the wrong answer?  The answer is that $\widehat{A}$ may not be computed 
from the exact $B$ and $C$, but from computed approximations $\widehat{B}$ 
and $\widehat{C}$.  The relative accuracy of 
$\widehat{A}$ can be worse than the relative accuracies of $\widehat{B}$
and $\widehat{C}$.  High relative accuracy can lost, quickly or slowly,
during a multi stage computation.  See Section 6, below.

\subsection{Introducing the standard}    \label{IEEE*se}

The IEEE floating point standard is a set of conventions  
on computer representation and processing of floating point numbers.  
Modern computers follow these standards for the most part.  
The standard has four main goals:
\begin{enumerate}
\item
To make floating point arithmetic as accurate as possible.
\item
To produce sensible outcomes in exceptional situations.
\item
To standardize floating point operations across computers.
\item
To give the programmer control over exception handling.
\end{enumerate}

The standard specifies exactly how numbers are represented in hardware.  
The most basic unit of information that a computer stores is a {\em bit}, 
a variable whose value may be either 0 or 1.  Bits are organized into $32$ 
bit or $64$ bit {\em words} words.  A $32$ bit word is a {\em string} of $32$ 
bits.  The number of $32$ bit words, or {\em bit strings}, is approximately 
(recall $2^{10} \approx 10^3$) 
$2^{32} = 2^2 \cdot 2^{30} \approx 4\times (10^3)^3 = 4$ billion.  A typical 
computer should take well under a minute to list all $32$ bit words.  A computer
running at $1\mbox{GHz}$ in theory can perform one billion operations
per second, though that is rarely achieved in practice.  The number
of $64$ bit words is about $1.6\cdot 10^{19}$ $64$, which is too many 
to be listed in a year.

There are two basic computer data types that represent numbers.
A fixed point number, or integer, has type {\tt int} or {\tt longint} 
in C/C++, depending on the number of bits.  Floating point 
numbers have type {\tt float} for {\em single precision} 32 bit words, 
and {\tt double} for {\em double precision}, or 64 bit words.  In early 
C compilers, a {\tt float} by default had 64 bits instead of 32.

Integer arithmetic, also called fixed point, is very simple.  For example,
with $32$ bit integers, the $4\cdot 10^9$ distinct words represent that 
many consequtive integers, filling the range from about $-2 \cdot 10^9$ 
to about $2 \cdot 10^9$.  Addition, subtraction, and multiplication are done 
exactly whenever the answer is within this range.  The result is 
unpredictable when the answer is out of range ({\em overflow}).  
Results of integer division are rounded down to the nearest integer 
below the answer.

\subsection{Representation of numbers, arithmetic operations}

For scientific computing, integer arithmetic has two drawbacks.  One is that 
there is no representation for numbers that are not integers.  Equally 
important is the small range of values.  The number of dollars in the US 
national debt, several trillion ($10^{12}$), cannot be represented as a 
32 bit integer but is easy to approximate in 32 bit floating point.

A floating point number is a binary version of the exponential (``scientific'')
notation used on calculator displays.  Consider the example expression:
\begin{displaymath}
-.2491{\mbox E}-5
\end{displaymath} 
which is one way a calculator could display the number $-2.491\cdot 10^{-6}$.
This expression consists of a sign bit, $s = -$, a mantissa, $m = 2491$ and an
exponent, $e = -5$.  The expression $s.m{\mbox E}e$ corresponds to the number
$s\cdot .m \cdot 10^e$.  Scientists like to put the first digit of the 
mantissa on the left of the decimal point ($-2.491\cdot 10^{-6}$) while
calculators put the whole thing on the right ($-.2491\cdot 10^{-5}$).
In base 2 (binary) arithmetic, the scientists' convention saves a bit, 
see below.

The IEEE format for floating point numbers replaces the base 10 with base 2,
and makes a few other changes.  When a 32 bit word (bit string) is interpreted 
as a floating point number, the first bit is the sign bit, $s = \pm$.  The next
8 bits form the ``exponent'', $e$, and the remaining 23 bits determine the 
``fraction'', $f$.  There are two possible signs, $2^8 = 256$ possible 
exponents (ranging from 0 to 255),  and $2^{23} \approx 8$ million possible 
fractions.  Normally a floating point number has the value
\begin{equation}
A = \pm 2 ^{e-127} \cdot \left( 1.f \right) _2 \;\; ,
\label{ieee*se}  \end{equation}
where $f$ is base 2 and the notation $\left( 1.f \right) _2$ means that the 
expression $1.f$
is interpreted in base 2.  Note that the mantissa is $1.f$ rather than just
the fractional part, $f$. Any number (except 0) can be normalized 
so that its base 2 mantissa has the form $1.f$.  There is no need to store the 
``$1.$'' explicitly, which saves one bit.

For example, the number $2.752 \cdot 10^3 = 2572$ can be written 
\begin{eqnarray*}
2752 & = & 2^{11} + 2^9 + 2^7 +2^6 \\
     & = & 2^{11} \cdot \left( 1 + 2^{-2} + 2^{-4} + 2^{-5} \right) \\
     & = & 2^{11} \cdot \left( 1 + \left( .01    \right)_2  
                                 + \left( .0001  \right)_2
                                 + \left( .00001 \right)_2 \right) \\
     & = & 2^{11} \cdot \left( 1.01011  \right)_2 \;\; .
\end{eqnarray*}
Altogether, we have, using $11 = (1011)_2$,
$$
2752 = + (1.01011)_2^{(1011)_2} \; .
$$
Thus, we have sign $s = +$.  The exponent is $e - 127 = 11$ 
so that $e = 138 = \left( 10001010 \right)_2$.
The fraction is $f = \left( 01011000000000000000000 \right)_2$.  The entire
32 bit string corresponding to $2.752 \cdot 10^3$ then is:
\begin{displaymath}
\underbrace{1}_{\displaystyle s} \underbrace{10001010}_{\displaystyle e}  
          \underbrace{01011000000000000000000}_{\displaystyle f} \;\;\;\; .
\end{displaymath}

For arithmetic, the standard mandates the rule: {\em the exact answer, 
correctly rounded}.  For example, suppose {\tt x}, {\tt y}, and {\tt z} 
are computer variables of type {\tt float}, and the computer executes the 
statement {\tt x = y / z;}.  Let $B$ and $C$ be the numbers that correspond 
to the $32$ bit strings {\tt y} and {\tt z} using the standard (\ref{ieee*se}).
Let $A$ be the {\em exact mathematical} quotient, $A=B/C$.  Let $\widehat{A}$
be the number closest to $A$ that can be represented exactly in the single
floating point format.  The computer is supposed to set the bit string
{\tt x} equal to the bit string representing that $\widehat{A}$.
For exceptions to this rule, see below.  

The {\em exact answer correctly rounded} rule implies that the only error 
in floating point arithmetic comes from rounding the exact answer, $A$, 
to the nearest floating point number, $\widehat{A}$.  Clearly, this rounding 
error is  determined by the distance between floating point numbers.  The 
worst case would be to put $A$ in the middle of the largest gap between
neighboring floating point numbers, $B$ and $B^{\prime}$.  For a floating 
point number of the form $B=(1.f)_2 \cdot 2^p$, the next larger floating 
point number is usually $B^{\prime} = (1.f^{\prime})_2 \cdot 2^p$, where 
we get $f^{\prime}$ from $f$ by adding the smallest possible fraction, which 
is $2^{-23}$ for $23$ bit single precision fractions.  The {\em relative} 
size of the gap between $B$ and $B^{\prime}$ is, after some algebra, 
$$
\epsilon = \frac{ B^{\prime} - B } { B } = 
\frac{ (1.f^{\prime})_2 - (1.f)_2 } { (1.f)_2 } = 
\frac{ 2^{-23} } { (1.f)_2 }\; .
$$
The largest $\epsilon$ is given by the smallest denominator, which is 
$(1.0\cdots 0)_2 = 1$, which gives $\epsilon_{\mbox{\em max}} = 2^{-23}$.
The largest rounding error is half the gap size, which gives the single
precision machine precision $\epsilon_{\mbox{\em mach}} = 2^{-24}$ stated above.

The $64$ bit double precision floating point format allocates one bit for 
the sign, $11$ bits for the exponent, and the remaining $52$ bits for the 
fraction.  Therefore its floating point precision is given by 
$\epsilon_{\mbox{\scriptsize mach}} = 2^{-53}$.
Double precision arithmetic gives roughly 16 decimal digits of
accuracy instead of 7 for single precision.  There are $2^{11}$
possible exponents in double precision, ranging from $1023$ to $-1022$.
The largest double precision number is of the order of 
$2^{1023} \approx 10^{307}$.  The largest single precision number is 
about $2^{126} \approx 10^{38}$.  Not only is double precision arithmetic
more accurate than single precision, but the range of numbers is far
greater.

\subsection{Exceptions}

The extreme exponents, $e=0$ and $e=255$ ($e=2^{11}-1 = 2047$ for double 
precision) do not correspond to numbers.  Instead, they have 
carefully engineered interpretations that make the IEEE standard distinctive.  
If $e=0$, the value is 
$$
A = \pm 0.f \cdot 2^{-126}  \;\; \mbox{(single precision),} \;\;\;\;\;\;
A = \pm 0.f \cdot 2^{-2046} \;\; \mbox{(double precision).}
$$
This feature is called {\em gradual underflow}.  {\em Underflow} is the 
situation in which the result of an operation is not zero but is closer 
to zero than any floating point number.  The corresponding numbers are
called {\em denormalized}.  Gradual underflow has the consequence that two
floating point numbers are equal, $x = y$, if and only if subtracting one
from the other gives exactly zero.

Introducing denormalized numbers makes sense when you
consider the spacing between floating point numbers.  If we exclude
denormalized numbers then the smallest positive floating pont number (in
single precision) is $A = 2^{-126}$ (corresponding to $e=1$ and 
$f = 00\cdots 00$ (23 zeros)) but the next positive floating point 
number larger than $A$ is $B$, which also has $e=1$ but now has 
$f=00\cdots 01$ (22 zeros and a 1).  Because of the implicit leading 1
in $1.f$, the difference between $B$ and $A$ is $2^{22}$ times smaller 
than the difference between $A$ and zero.  That is, without gradual 
underflow, there is a large and unnecessary gap between
0 and the nearest nonzero floating point number.

The other extreme case, $e = 255$, has two subcases, {\tt inf} (for infinity)
if $f = 0$ and {\tt NaN} (for Not a Number) if $f \neq 0$.  C/C++ 
prints\footnote{In keeping with the Microsoft pattern of maximizing 
incompatibility, the its compiler prints something different.}
 ``{\tt inf}'' and ``{\tt NaN}'' when you print out a variable in
floating point format that has one of these values.  An arithmetic 
operation produces {\tt inf} if the exact answer is larger than the 
largest floating point number, or {\tt 1/x} if $x = \pm 0$.  (Actually $1/+0 =$ 
+{\tt inf} and $1/-0 =$ -{\tt inf}).  Invalid operations such as 
{\tt sqrt(-1.)}, {\tt log(-4.)}, produce {\tt NaN}.  Any operation involving 
a {\tt NaN} produces another {\tt NaN}.  It is planned that $f$ will contain 
information about how or where in the program the {\tt NaN} was created but 
this is not standardized yet.  Operations with {\tt inf} are
common sense: ${\mbox {\tt inf}} + {\mbox \it finite} = {\mbox {\tt inf}}$,
${\mbox{\tt inf}} / {\mbox {\tt inf}} = {\mbox {\tt NaN}}$, 
${\mbox \it finite} / {\mbox {\tt inf}} = 0.$, 
${\mbox {\tt inf}} - {\mbox {\tt inf}} = {\mbox {\tt NaN}}$.

A floating point arithmetic operation is an {\em exeption} if the result is not
a normalized number in floating point format.  The standard mandates that
a hardware {\em flag} (a binary bit of memory in the processor) should
be {\em set} (given the value 1) when an exception occurs.  There should
be a separate flag for the underflow, {\tt inf}, and {\tt NaN} exceptions.
The programmer should be able to specify what happens when an exception
flag is set.  Either the program exection continues without interruption
or an {\em exeption handler} procedure is called.  The programmer should
be able to write procedures that interface with the exception handler to
find out what happened and take appropriate action.  Only the most advanced
and determined programmer will be able to do this.  The rest of us have the
worst of both: the exception handler is called, which slows the program
execution but does nothing useful.

Many features of IEEE arithmetic are illustrated in 
Figure \ref{SourcesOfError.C*se}.  Note that $e^{204}$ gives {\tt inf} 
in single precision but not in double precision because the range of 
values is larger in double precision.  We see that {\tt inf} and {\tt NaN}
work as promised.  The main rule, ``exact answer correctly
rounded'', explains why adding pairs of floating point numbers is
commutative: the mathematical sums are equal so they round to the
same floating point number.  This does not force addition to be 
associative, which it is not.  Multiplication is also commutitative
but not associative.  The division operator gives integer or floating 
point division depending on the types of the operands.  Integer arithmetic
{\em truncates} the result to the next lower integer rather than rounding
it to the nearest integer.

\rule{12cm}{1mm}
\begin{verbatim}
// A program that explores floating point arithmetic in the IEEE 
// floating point standard.  The source code is SourcesOfError.C.

#include <iostream.h>
#include <math.h> 

int main() {

   float  xs, ys, zs, ws;  // Some single precision variables.
   double yd;              // A double precision variable.

   
   xs = 204.;         // Take an exponential that is out of range.
   ys = exp(xs);
   cout << "The single precision exponential of " << xs << 
           " is " << ys << endl;
   yd = exp ( xs );   // In double precision, it is in range.
   cout << "The double precision exponential of " << xs << 
           " is " << yd << endl;
 
   
   zs = xs / ys;      // Divide a normal number by infinity.
   cout << xs << " divided by " << ys  << " gives " << zs << endl;
   
   ws = ys;           // Divide infinity by infinity.
   zs = ws / ys;
   cout << ws << " divided by " << ys << " gives " << zs << endl;
   
   zs = sqrt( -1.) ;  // sqrt(-1) should be NaN.
   cout << "sqrt(-1.) is " << zs << endl;
   
   ws = xs + zs;      // Add NaN to a normal number.
   cout << xs << " + " << zs << " gives " << ws << endl;
   
   xs =        sin(1.); // Some generic single precision numbers.
   ys = 100.  *sin(2.);
   zs = 10000.*sin(3.);
   float xsPys, ysPxs, xsPzs, zsPxs;  // xsPzx holds xs + zs, etc.
   xsPys = xs + ys;   
   ysPxs = ys + xs;   // Try commuting pairs.
   xsPzs = xs + zs;
   zsPxs = zs + xs;  
   if ( ( xsPys == ysPxs ) && ( xsPzs == zsPxs ) ) 
      cout << "Adding " << xs << "  " << ys << " and "<< zs <<
              " in pairs commutes." << endl;
     else
       cout << "Adding " << xs << "  " << ys << " and "<< zs <<
               " in pairs does not commute." << endl; 
              
   float xsPysPzs, ysPzsPxs;     // Test for associativity.
   xsPysPzs = ( xs + ys ) + zs;
   ysPzsPxs = ( ys + zs ) + xs;
   if ( xsPysPzs == ysPzsPxs ) 
      cout << "Adding " << xs << "  " << ys << " and "<< zs <<
              " is associative." << endl;
     else 
       cout << "Adding " << xs << "  " << ys << " and "<< zs <<
               " is not associative." << endl;
               
   int xi, yi; // Some integer variables.
   xi =  9;    // Compute the quotient using integer
   yi = 10;    // and floating point arithmetic.
   zs = xi/yi;  
   ws = ( (float) xi ) / ( (float) yi ); // Convert, then divide.
   cout << "Integer division of " << xi << " by " << yi << " gives " <<
           zs << ".  Floating point gives " << ws << endl;
        
   return(0);
   
   }
\end{verbatim}
\begin{figure}[h]
\rule{12cm}{1mm}
\caption{\it A program that illustrates some of the features of arithmetic
using the IEEE floating point standard.}
\label{SourcesOfError.C*se}
\end{figure}

\section{Truncation error}             \label{TruncErr*se}
Truncation error is the error in analytical approximations such as 
\begin{equation}
f^{\prime}(x) \approx \frac{f(x+h) - f(x)}{h}\;.
\label{fp*se}  \end{equation}
This is not an exact formula, but it can be a useful approximation.  We 
often think of truncation error as arising from truncating a Taylor series.  
In this case, the Taylor series formula,
$$
f(x+h) = f(x) + hf^{\prime}(x) + \frac{1}{2} h^2 f^{\prime\prime}(x) + \cdots
\; ,
$$
is truncated by neglecting all the terms after the first two on the right.
This leaves the approximation
$$
f(x+h) \approx f(x) + hf^{\prime}(x)  \; ,
$$
which can be rearranged to give (\ref{fp*se}).  Usually, truncation is the 
main source of error in numerical integration or solution of differential 
equations.  The analysis of truncation error using Taylor series will
occupy the next two chapters.  

As an example, we take $f(x) = xe^x$, $x=1$, and a variety of $h$ values.
The results are in Figure \ref{TrunkTable*se}.  In Chapter \ref{na*na}
we will see that in exact arithmetic (i.e.\ without rounding), the error
would be roughly proportional th $h$ for small $h$. These numbers were computed
in double precision arithmetic and the effect of roundoff error is 
appearent in for the smallest $h$ values.  It is rare in a practical 
application that $h$ would be so small that roundoff error would be a
significant factor in the overall error.  
\begin{figure}
\begin{tabular}{|c|c|c|c|c|c|c|} \hline
        $h$            & .3   &       .01           &       $10^{-5}$     &    $10^{-8}$         &   $10^{-10}$         \\  \hline
$\widehat{f^{\prime}}$ & 6.84 &      5.48           &        5.4366       &      5.436564        &     5.436562         \\  \hline
$  e_{\mbox{\em tot}}$ & 1.40 & $4.10\cdot 10^{-2}$ & $4.08\cdot 10^{-5}$ & $-5.76\cdot 10^{-8}$ & $-1.35\cdot 10^{-6}$ \\  \hline
\end{tabular}
\label{TrunkTable*se}
\caption{Estimates of $f^{\prime}(x)$ using (\ref{fp*se}).  The error
is $e_{tot}$, which is a combination of truncation and roundoff error.
Roundoff error is apparent only in the last two estimates.}
\end{figure}

\section{Iterative Methods}              \label{IterErr*se}
Suppose we want to find a number, $A$, by solving an equation.  Often it is
impossible to find a formula for the solution.  Instead, {\em iterative
methods} construct a sequence of approximate solutions, $A_n$, for
$n=1, 2, \ldots\;$.  Hopefully, the approximations {\em converge} to the 
right answer: $A_n \to A$ as $n \to \infty$.  In practice, we must stop
the iteration process for some large but finite $n$ and accept that $A_n$
as the approximate answer.  

For example, suppose we have a $y>0$ and we want to find $x$ with 
$xe^x=y$.  There is not a formula for $x$, but we can write a program
to carry out the iteration: $x_1 = 1$, $x_{n+1} = \ln(y) - \ln(x_n)$.
The numbers $x_n$ are the {\em iterates}.
If the limit of the iterates exists, $x_n \to x$ exists as $n \to \infty$, 
then that $x$ should be a {\em fixed point} of the iteration, i.e.\ 
$x = \ln(y) - \ln(x)$.  Figure \ref{ItTable*se} demonstrates the
convergence of the iterates in this case.  If we stop after $10$
iterations, the $20^{\mbox{\em th}}$ iterate has an error 
$e_{20} \approx 2.3\cdot 10^{-5}$, which might be small enough,
depending on the application.
\begin{figure}
\begin{tabular}{|c|c|c|c|c|c|c|} \hline
$n$   &     1   &     3   &        6           &   10               & 20                &    50              \\  \hline
$x_n$ &     1   &  $1.46$ &    $1.80$          & $1.751$            & $1.74555$         & $1.745528$         \\  \hline
$e_n$ & $-.745$ & $-.277$ & $5.5\cdot 10^{-2}$ & $5.9\cdot 10^{-3}$ & $2.3\cdot 10^{-5}$& $1.2\cdot 10^{-12}$\\  \hline
\end{tabular}
\label{ItTable*se}
\caption{Iterates of $x_{n+1} = \ln(y) - \ln(x_n)$ illustrating convergence to a limit 
that satisfies the equation $xe^x = y$.  The error is $e_n = x_n - x$.  Here, $y=10$.}
\end{figure}

\section{Statistical error in Monte Carlo}    \label{StatErr*se}
Monte Carlo is a computational technique that uses random numbers as a
computational tool.  For example, we may be able to use the computer
random number generator to create a sequence of random variables
$X_1$, $X_2$, $\ldots$, with the same distribution.  We could use these
samples to estimate the average, $A$, by the computer generated sample
mean
$$
A \approx \widehat{A} = \frac{1}{n} \sum_{k=1}^n X_k \; .
$$
Statistical error is the difference between the estimate $\widehat{A}$ and 
the correct answer, $A$.  Monte Carlo statistical errors typically would be 
larger than roundoff or truncation errors.  This makes Monte Carlo a method 
of last resort, to be used only when other methods are not available.

We illustrate the behavior of statistical error with a simple example.
Here the $X_k$ are independent random variables with the same probability
distribution that have average value $A = .5$.  Figure \ref{StatTable*se}
gives the approximations and statistical errors for several values of
$n$.  The value $n=10^6$ is repeated several times to illustrate the fact
that statistical error is random (see Chapter {mc*mc} for a clarification
of this).  Note the size of the errors even with a million samples and compare
these to the errors in Figures \ref{} and \ref{}.

\begin{figure}
\begin{tabular}{|c|c|c|c|c|c|c|} \hline
        $n$   &  10   &    100             &  $10^4$            &    $10^6$          &   $10^6$            &    $10^6$          \\  \hline
$\widehat{A}$ & .603  &   .518             &  .511              &    .5004           &    .4996            &    .4991           \\  \hline
     error    & .103  & $1.8\cdot 10^{-2}$ & $1.1\cdot 10^{-2}$ & $4.4\cdot 10^{-4}$ & $-4.0\cdot 10^{-4}$ & $-8.7\cdot 10^{-4}$\\  \hline
\end{tabular}
\label{StatTable*se}
\caption{Statistical errors in a demonstration Monte Carlo computation.}
\end{figure}

\section{Error amplification and unstable algorithms}     \label{Stab*se}

Errors can be amplified during sequence of computational steps.  For example,
suppose we are using the divided difference (\ref{fp*se}) to approximate 
$f^{\prime}(x)$ using approximate values $\widehat{f}_1 = f(x) + e_1$ and 
$\widehat{f}_2 = f(x+h) + e_2$.  This is practically inevitable, since we
cannot evaluate $f$ exactly.  Using these values gives
\begin{eqnarray*}
\widehat{f^{\prime}} 
& = & \frac{ f(x+h) - f(x) }{h} + \frac{e_2 - e_1}{h} + e_r \\
  & = & f^{\prime}(x) + e_{\mbox{\em tr}} + e_{\mbox{\em am}} + e_r\; .
\end{eqnarray*}
Here $e_{\mbox{\em tr}}$ is the truncation error 
$$
e_{\mbox{\em tr}} = \frac{f(x+h) - f(x)}{h} - f^{\prime}(x) \; ,
$$
$e_r$ is the roundoff error in evaluating $(\widehat{f_2} - \widehat{f_1})/h$, 
in floating point, given the floating point values  
$\widehat{f_1}$ and $\widehat{f_1}$.  The remaining error, and often the
largest, is the amplification error
$$
e_{\mbox{\em am}} = \frac{e_2 - e_1}{h} \;\; .
$$
We see that the errors $e_1$ and $e_2$ are {\em amplified} by the factor
$1/h$, which is quite large if $h$ is small.  For example, if $h=.01$,
then the error in evaluating $f^{\prime}$ is likely to be a hundred times
larger than the error in approximating $f$.  

This point may be put in a more technical form: addition or subtraction 
with cancellation can amplify relative errors in the operands.  This is 
only true for relative error, which is a major reason relative error is
important.  Consider the operations $A=B+C$ and 
$\widehat{A} = \widehat{B} + \widehat{C}$.  We suppose for simplicity that 
the sum is done exactly, so that all the error in $\widehat{A}$ is due to error 
in the operands, $\widehat{B}$ and $\widehat{C}$.  The absolute and relative 
errors, $e_A$, $\epsilon_A$, etc.\, are defined by (see (\ref{error*se}):
$$
\widehat{A} = A + e_A = (1+\epsilon_A)A \;\; , \;\;\;\;
\widehat{B} = B + e_B = (1+\epsilon_B)B \;\; , \;\;\;\;
\widehat{C} = C + e_C = (1+\epsilon_C)C \;\; .
$$
With absolute errors there are no surprises: $e_A = e_B + e_C$.  We do not
know the signs of $e_B$ and $e_C$, but we can write
$|e_A| \leq |e_B| + |e_C|$.  This says that the magnitude of the 
absolute error of the result is at most the sum of the absolute errors
in the operands.  If floating point arithmetic were to produce high absolute
accuracy, which it does not, we could perform tens of thousands of single 
precision operations always come out with high absolute accuracy.  

With relative error, {\em amplification} through {\em cancellation} is possible.
The relative error of $A$ may be larger or much larger than sum of the 
relative errors in the operands.  Some algebra leads to the formula
\begin{equation}
\left| \epsilon_A \right| 
= \left| \epsilon_B W_B + \epsilon_C W_C \right| 
\leq \max(\left| \epsilon_B \right|,\left| \epsilon_C \right|) 
(\left|W_B\right| + \left|W_C\right|) \;\; ,
\label{error_amp*se}  \end{equation}
where
$$
W_B = \frac{B}{B+C} \;\;, \;\;\;\;
W_C = \frac{C}{B+C} \;\;, \;\;\;\;
W_B + W_C = 1 \; .
$$
If $B$ and $C$ have the same sign, then 
$\left|W_B\right| = \left|W_C\right| = 1$, and (\ref{error_amp*se}) shows
that $\left|\epsilon_A\right|$ is less than both $\left|\epsilon_A\right|$
and $\left|\epsilon_A\right|$: the relative error has not been amplified.
On the other hand, if $B=1.01$ and $C=-1$, then 
$\left|W_B\right| \approx \left|W_C\right| = 100$: the relative error 
may be amplified by a factor of $100$.  

{\em Cancellation} is the situation $A=B+C$ (or $A=B-C$) and 
$|A| < |B| + |C|$.  Addition has cancellation if the operands $B$ and 
$C$ have opposite signs.  The previous paragraph shows that addition 
and subtraction amplify relative errors only if there is cancellation.
A similar argument shows that multiplication and division produce high
relative accuracy in (almost) all cases, no cancellation there.  
Performing $1,000$ arithmetic operations will generate a relative
error not larger than $1000\epsilon_{\mbox{\em mach}}$, and probably 
much smaller, if there is no cancellation in the additions and subtractions.

It is possible to lose many decimal places of accuracy in a single subtraction.
This is called {\em catistrophic cancellation}.  More insideous is the 
possiblilty of small error amplification at each stage of a long calculation.
For example, if we start with a relative error of $10^{-16}$ characteristic
of double precision arithmetic and amplify by $5\%$ in each of $1000$ steps,
then we have well over $100\%$ relative error at the end.  A computational
algorithm is called {\em unstable} if it loses accuracy because of 
cancellation, either catastrophically or, which is worse, gradually.
One of the main uses of mathematical analysis in scientific
computing is in understanding the stability of multistage computational
algorithms.
 
Inaccuracy because of an unstable algorithm very
is hard to discover by standard debugging techniques.  In a large calculation,
the error may grow a seemingly negligible amount at each step but
grow to swamp the correct answer by the time the computation is
finished. 

\section{Condition number and ill conditioned problems}   \label{Cond*se}

The {\em condition number} of a computational problem measures the sensitivity
of the answer to small changes in the data, in relative terms.  The 
condition number limits the accuracy we can expect to achieve in a computation
even with the best algorithm.  If a problem is sufficiently 
{\em ill conditioned}, the condition number is too large, then even 
inevativle rounding errors in the input may lead to unacceptable errors in $A$.

We write $A(x)$ to indicate that the answer, $A$, depends on the data, $x$.  
In the simplest case, both $A$ and  $x$ are single numbers, and $A(x)$ is a
differentiable function of $x$.  A perturbation in the data, $\Delta x$, 
will lead to a change in the answer, $\Delta A = A(x+\Delta x) - A(x)$.  
If $\Delta x$ is small enough, $\Delta A/\Delta x \approx a^{\prime}(x)$, 
so that
$$
\Delta A \approx A^{\prime}(x)\Delta x\; .
$$
The absolute change in the answer is related to the absolute change in the
data by the derivative, $A^{\prime}$.

The relative change of the answer is related to the relative change in the
data by the condition number, $\kappa$:
\begin{equation}
\frac{\Delta A}{A} \approx \kappa \frac{\Delta x}{x} \;\; .
\label{oneDk1*se}  \end{equation}
For small $\Delta x$, we may use the derivative approximation to get, after
a bit of algebra,
$$
\kappa = \frac{\Delta A/A}{\Delta x / x} 
\approx A^{\prime} \cdot \frac{x}{A}  \;\;\; .
$$
Everybody's definition of $\kappa$ differs from this one in that nobody allows
the condition number to be negative; they use the absolute value of this 
$\kappa$:
\begin{equation}
\kappa = \left| A^{\prime}(x) \frac{x}{A(x))} \right| \;\; .
\label{kapDef*se}    \end{equation}
This has the property that the ratio of the relative size of $|\Delta A|$ 
to the relative size of $\Delta x$ is not very sensitive to $\Delta x$,
provided $\Delta x$ is small enough:
\begin{equation}
\frac{|\Delta A|/|A|}{|\Delta x|/|x|} \approx \kappa \;\;, \;\;\;\;
\mbox{for small $\Delta x$.}
\label{oneDk2*se}    \end{equation}

In most real problems, both the answer and the data are more complex than 
a single number.  We still write $\Delta x$ and $\Delta A$ for perturbations
in the data and answer, but now norms $\left\| \Delta x \right\|$ and
$\left\|\Delta A\right\|$ measure the size of the perturbation.  We will 
discuss norms in Chapter \ref{lc*lc}.  A measure of the relative size of 
$\Delta x$ is $\left\|\Delta x \right\| / \left\|  x \right\|$.  Unlike the
univariate case (\ref{oneDk2*se}), here $\left\| \Delta A \right\|$ depends
on the direction of $\Delta x$, not just its size.  The condition number
is taken to be the worst case sensitivity:
$$
\kappa = \max_{\Delta x \;\mbox{\scriptsize small}} 
\frac{ \left\| \Delta A\right\|/\left\|A\right\| }
     { \left\| \Delta x\right\|/\left\|x\right\| } \;\; .
$$
We rewrite this to resemble (\ref{oneDk1*se}) by saying that $\kappa$ is 
the smallest number so that
\begin{equation}
\frac{\left\|\Delta A\right\|}{\left\|A\right\|} \leq \approx \kappa
\frac{\left\|\Delta x\right\|}{\left\|x\right\|} \;\;,\;\;\;\;
\mbox{for all small enough $\Delta x$.}
\label{k2*se}   \end{equation}

Floating point rounding makes it is inevatible that $|\Delta x/x|$ is of 
the order of $\epsilon_{\mbox{\em mach}}$.  Therefore, if 
$\kappa > 1/\epsilon_{\mbox{\em mach}}$, than (\ref{oneDk1*se}) or 
(\ref{k2*se}) show that $\Delta A$ will be as large as $A$, total error.  
Unfortunately, condition numbers as large as $10^7$ or $10^{16}$ occur in 
practice.  

Computational errors caused by large condition number cannot be cured by
changing the computational algorithm.  None of our discussion of condition
number referred to any particular computational algorithm.  If a problem
has 

Anyone engaged in scientific computing should keep conditioning in mind.
He or she should be aware that a problem may be ill conditioned, and try
to estimate the condition number, at least roughly.  High quality numerical
software often includes estimators of condition number with the results,
which should be checked.  Even in complex situations, it may be possible to
find an analogue of (\ref{kapDef*se}) by differentiation.  Chapter 
\ref{lc*lc} has some examples.

\section{Software tips}       \label{Tips*se}

\subsection{Floating point numbers are (almost) never equal}

Because of inexact floating point arithmetic, two numbers that should be
equal in exact arithmetic often are not equal in the computer.  For example,
the code in Figure \ref{FloatingPointInequality*se} prints {\tt Not equal}.
\begin{figure}
\rule{12cm}{1mm}
\begin{verbatim}
   double x = ( double(1) ) / 7, y1, y2;
   y1 = x + x + x + x + x;
   y2 = 5*x;
   if ( y1 == y2 ) cout << "Equal"     << endl;
   else            cout << "Not equal" << endl;
\end{verbatim}
\rule{12cm}{1mm}
\caption{\it A code fragment in which roundoff error can lead to numbers that
should be equal in exact arithmetic not being equal in floating point.}
\label{FloatingPointInequality*se}
\end{figure}
In general, an equality test between two variables of type {\tt float}
or {\tt double} is a mistake.  Even worse, if we have {\tt y1 = y2;} and
then do not reassign either {\tt y1} or {\tt y2}, many lines later 
the test {\tt ( y1 == y2 )} may evaluate to {\tt false}.  For example,
on a Pentium 4 chip there are $80$ bit registers for holding variables
of type {\tt double}.  The computer often holds a variable in a register
rather than returning it to memory to save time.  If the variable is
returned to memory, only the $64$ bits of the IEEE standard are stored,
the remaining $16$ bits being lost.  If {\tt y1} but not {\tt y2} was
returned to memory, then they may no longer be equal. 

A common mistake in this regard is to use floating point comparisons to
terminate a loop.  Figure \ref{floatLoopTest*se} illustrates this.
In exact aritmetic this would give the desired $n$ iterations.  Because
of floating point arithmetic, after the $n^{th}$ iteration, the variable
{\tt t} may be equal to {\tt tFinal} but is much more likely to be above
or below roundoff error.  Unfortunately, it is impossible to predict which
way the roundoff error will go.  Therefore, we do not know whether this 
code will execute the {\tt while} loop body $n$ or $n+1$ times.
\begin{figure}
\rule{12cm}{1mm}
\begin{verbatim}
   double tStart, tFinal, t, dt;
   int    n;
   tStart = . . .  ;   // Some code that determines the start
   tFinal = . . .  ;   // and ending time and the number of
   n      = . . .  ;   // equal size steps.
   dt     = ( tFinal - tStart ) / n;  // The size of each step.
   while ( t = tStart, t < tFinal, t+= dt ) 
      { . . . }        // Body of the loop does not assign t.
\end{verbatim}
\rule{12cm}{1mm}
\caption{\it A code fragment illustrating a pitfall of using a floating point
variable to regulate a {\tt while} loop.}
\label{floatLoopTest*se}
\end{figure}

To guarantee $n$ executions of the {\tt while} loop body, use integer (fixed
point) variables, as illustrated in Figure \ref{intLoopTest*se}
\begin{figure}
\rule{12cm}{1mm}
\begin{verbatim}
   double tStart, tFinal, t, dt;
   int    n,    , i;
   tStart = . . .  ;   // Some code that determines the start
   tFinal = . . .  ;   // and ending time and the number of
   n      = . . .  ;   // equal size steps.
   dt     = ( tFinal - tStart ) / n;  // The size of each step.
   while ( i = 0, i < n, i++ ) 
      { t = tStart + i*dt;  // In case the value of t is needed 
             . . . }        // in the loop body.   
\end{verbatim}
\rule{12cm}{1mm}
\caption{\it A code fragment using an integer
variable to regulate the loop of Figure \ref{floatLoopTest*se}.}
\label{intLoopTest*se}
\end{figure}

\subsection{Plotting data curves}
Careful visualization is a key step in understanding any data, particularly
the results of a numerical computation.  Time spent making plots carefully
and thoughtfully is usually rewarded.  The first rule is that a plot, like
a piece of code, should be self documenting.  In coding, we take the time to
type comments, format carefully, etc.  In plots, we take the time to label
axes and print titles with relevant plot parameters.  Otherwise, a week later
nobody will know which plot is which.

Careful scaling can make the data more clear.  If you have a function $f(x)$
with values in the range from $1.2$ to $1.22$, the Matlab {\tt plot} function
will lable the vertical axis from $0$ to $2$ and plot the data as a nearly
horizontal line.  This is illustrated in Figure \ref{PlotScaling*se}, where
the first $70$ Fibonacci numbers, are plotted on a linear scale and a log 
scale.  The Fibonacci numbers, $f_i$, are defined by $f_0 = f_1 = 1$, 
and $f_{i+1} = f_i + f_{i-1}$, for $i\geq 1$.  On the linear scale, $f_1$ 
through $f_{57}$ sit on the horizontal axis.
The log plot lets us see how big each of the $70$ numbers is.  It also makes
it clear that $\log(f_i)$ is nearly proportional to $i$.  If 
$\log(f_i) \approx a + bi$, then $f_i \approx cd^i$.  It is common to understand
the data by plotting it in various ways.  
\begin{figure}
\scalebox{.4}[.4]{\includegraphics[width=\linewidth]
          {MatlabScripts/FibLinear_se} }
\hspace{1cm}
\scalebox{.4}[.4]{\includegraphics[width=\linewidth]
          {MatlabScripts/FibLog_se} }
\caption{\it Plots of the first $n$ Fibonacci numbers, linear scale on the 
left, log scale on the right}
\label{PlotScaling*se}
\end{figure}

The Matlab script that made the plots of Figure \ref{PlotScaling*se} is
in Figure \ref{PlotScalingScript*se}.  Normally, one need not comment up
throwaway scripts like this.  The only real parameters are $n$, the largest
$i$ value, and whether the plot is on a linear or log scale.  Both of those
are recorded in the plot.  Note the convenience and clarity of not hard 
wiring $n=70$.  It would take just a moment to make plots up to $n=100$.  
   
\begin{figure}
\rule{12cm}{1mm}
\begin{verbatim}
     % Matlab code to generate and plot Fibonacci numbers.

   clear f     % If you decrease the value of n, the plots still work.
   n    = 70;  % The number of Fibonacci numbers to compute.
   fi   = 1;   % Start with f0 = f1 = 1, as usual.
   fim1 = 1;
   f(1) = fi;  % Record f(1) = f1.
   for i = 2:n
      fip1 = fi + fim1;  % f(i+1) = f(i) + f(i-1) is the recurrence
      fim1 = fi;         % relation that defines the Fibonacci numbers.
      fi   = fip1;
      f(i) = fi;         % Record f(i) for plotting.
   end

   plot(f)
   xlabel('i')           % The horizontal and vertical axes are 
   ylabel('f')           % i and f respectively.
   topTitle = sprintf('Fibonacci up to n = %d',n); % Put n into the title.
   title(topTitle)
   text(n/10, .9*f(n), 'Linear scale');
   grid                  % Make it easier to read values in the plot.
   set ( gcf, 'PaperPosition', [.25 2.5 3.2 2.5]); 
                      % Print a tiny image of the plot for the book.
   print -dps FibLinear_se
\end{verbatim}
\rule{12cm}{1mm}
\caption{\it A code fragment using an integer
variable to regulate the loop of Figure \ref{floatLoopTest*se}.}
\label{PlotScalingScript*se}
\end{figure}

\section{Further reading}
The idea for starting a book on computing with a discussion of sources
of error comes from the book {\em Numerical Methods and Software} by
David Hahaner, Cleve Moler, and Steve Nash.  Another interesting version
is in {\em Scientific Computing} by Michael Heath.  My colleague, Michael
Overton, has written a nice short book {\em IEEE floating point arithmetic}.


\section{Exercises}

\begin{enumerate}

\item  \label{TF*se}
For each of the following, indicate whether the statement is true
or false and explain your answer.
\begin{description}
\item[a.] The {\tt for} loop in line 4 below is will be executed 20 times 
no matter what $x$ is generated by {\tt rand()}.  Here {\tt float rand()} 
generates random numbers uniformly distributed in the interval $[0,1]$.
\begin{verbatim}
     float  x = 100*rand() + 2;                // line 1
     int    n = 20;                            // line 2
     float dy = x/n;                           // line 3
     for ( float y = 0; y < x; y += dy; ) {    // line 4
          body does not change x, y, or dy }  
\end{verbatim}
\item[b.] The {\tt for} loop in line 4 always will always be executed at least
20 times.
\item[c.] The {\tt for} loop in line 4 never will be executed more than
20 times.
\item[d.] There is a function, $f(x)$ that we wish to compute numerically.
We know that for $x$ values around $10^{-3}$, $f$ is about $10^5$
and $f^{\prime}$ is about $10^{10}$.  This function is too ill conditioned
to compute in single precision.
\end{description}

\item Show that in the IEEE floating point standard with any number
of fraction bits, $\epsilon_{\mbox{\em mach}}$ is the largest floating 
point number, $\epsilon$, so that $1+\epsilon$ gives $1$ in floating 
point arithmetic.  Whether this is mathematically equal to the definition
in the text depends on how ties are broken in rounding, but the difference
between the two definitions is irrelevent (show this).

\item  Starting with the declarations
\begin{verbatim}
           float x, y, z, w;
     const float oneThird = 1/ (float) 3;  // const means these can never
     const float oneHalf  = 1/ (float) 2;  // be reassigned
\end{verbatim}
we do lots of arithmetic on the variables {\tt x, y, z, w}.  In each case
below, determine whether the two arithmetic expressions result in the same
floating point number (down to the last bit) as long as no {\tt NaN} or
{\tt inf} values or denormalized numbers are produced.
\begin{description}
\item[a.] $  $
\begin{verbatim}
      ( x * y ) + ( z - w )
      ( z - w ) + ( y * x )
\end{verbatim}
\item[b.] $  $
\begin{verbatim}
      ( x +   y ) + z 
        x + ( y +   z )
\end{verbatim}
\item[c.] $  $
\begin{verbatim}
        x * oneHalf + y * oneHalf
      ( x + y ) * oneHalf
\end{verbatim}
\item[d.]
\begin{verbatim}
        x * oneThird + y * oneThird
      ( x + y ) * oneThird
\end{verbatim}
\end{description}

\item \label{Fib*se}
The {\em fibonacci numbers}, $f_k$, are defined by $f_0 = 1$, $f_1 = 1$, and 
\begin{equation}
f_{k+1} = f_k + f_{k-1}
\label{fibEqn*se}   \end{equation}
for any integer $k>1$.  A small perturbation of them, the {\em pib numbers} 
(``p'' instead of ``f'' to indicate a perturbation), $p_k$, are defined by 
$p_0 = 1$, $p_1 = 1$, and 
$$
p_{k+1} = c \cdot p_k + p_{k-1}
$$
for any integer $k>1$, where $c = 1+\sqrt{3}/100$.  
\begin{description}
\item[a.]
Make one plot of $\log(f_n)$ and $\log(p_n)$, or plot the $f_n$ and $p_n$ 
together on a log scale, as a function of $n$.  On the 
plot, mark $1/\epsilon_{\mbox{\em mach}}$ for single and double 
precision arithmetic.  This can be useful in answering the questions below.
\item[b.]
Rewrite (\ref{fibEqn*se}) to express $f_{k-1}$ in 
terms of $f_k$ and $f_{k+1}$.  Use the computed $f_n$ and $f_{n-1}$ to 
recompute $f_k$ for $k = n-2, n-3, \ldots, 0$.  Make a plot of the difference
between the original $f_0 = 1$ and the recomputed $\widehat{f}_0$ as a 
function of $n$.  What $n$ values result in no accuracy for the recomputed 
$f_0$?  How do the results in single and double precision differ?
\item[c.]
Repeat b.\ for the pib numbers.  Comment on the striking difference in
the way precision is lost in these two cases.  Which is more typical?
{\em Extra credit}: predict the order of magnitude of the error in recomputing
$p_0$ using what you may know about recurrence relations and what you
should know about computer arithmetic.
\end{description}

\item  \label{Binom*se}
The binomial coefficients, $a_{n,k}$, are defined by
$$
a_{n,k} = \left( \begin{array}{c} n \\ k \end{array} \right) 
        = \frac{n! } {k! (n-k)!}
$$
To compute the $a_{n,k}$, for a given $n$, start with $a_{n,0} = 1$
and then use the recurrence relation $a_{n,k+1} = \frac{n-k}{k+1}a_{n,k}$.
\begin{description}
\item[a.]
For a fixed of $n$, compute the $a_{n,k}$ this way, noting the largest 
$a_{n,k}$ and the accuracy with which $a_{n,n}=1$ is computed.  Do this
in single and double precision.  Why is it that roundoff is not a 
problem here as it was in problem \ref{Fib*se}?
\item[b.]
Use the algorithm of part (a) to compute
\begin{equation}
E(k) = \frac{1}{2^n} \sum_{k=0}^n k a_{n,k} = \frac{n}{2}\;\; .
\label{E(k)*se}    \end{equation}
In this equation, we think of $k$ as a random variable, the number
of heads in $n$ tosses of a fair coin, with $E(k)$ 
being the expected value of $k$.  This depends on $n$.  
Write a program without any safeguards against overflow or zero divide 
({\em this time only!})\footnote{One of the purposes of the IEEE
floating point standard was to {\em allow} a program with overflow or
zero divide to run and print results}.  
Show (both in single and double precision)
that the computed answer has high accuracy as long as the intermediate 
results are within the range of floating point numbers.  As with (a), 
explain how the computer gets an accurate, small, answer when the 
intermediate numbers have such a wide range of values.  Why is cancellation
not a problem?  Note the advantage of a wider range of values: we can
compute $E(k)$ for much larger $n$ in double precision.  Print $E(k)$
as computed by (\ref{E(k)*se}) and $M_n = \max_k a_{n,k}$.  For large $n$, one 
should be {\tt inf} and the other {\tt NaN}.  Why?
\item[c.]
For fairly large $n$, plot $a_{n,k}$ as a function of $k$ for a range of
$k$ chosen to illuminate the interesting ``bell shaped'' behavior 
of the $a_{n,k}$ near their maximum.
\item[d.] ({\em Extra credit, and lots of it!}) Find a way to compute
\begin{displaymath}
S(k) =  \sum_{k=0}^n \left(-1\right)^k\sin(2\pi\sin(k/n)) a_{n,k} 
\end{displaymath}
with good relative accuracy for large $n$.  This is very hard, so don't
spend too much time on it.
\end{description}
\item What day is it
\label{day}
\end{enumerate}


