📢 Updates

Lecture 9: Hermite interpolation available Assignment 3 available Lecture 8 online Lecture 7 online Schedule now available in "Admin" Problem class exercises now online Lecture 6 online Quiz 1 solutions online Lectures 2-3 exercises have solutions! Lecture 1 exercises have solutions!

2  How computers add

Lecture 2: Floating point arithmetic

Overview

1. Review of last week,
   
2. Whiteboard I (Standard/scientific form; general floating point numbers; underflow, overflow, machine epsilon; Examples: half, single, double precision)
3. Half precison = Float16 (Julia)
4. Whiteboard II (round-off error, relative error in floating point representations and the distribution of floating-point numbers; Conditioning, Stability)

Warning

These lecture notes are mainly a record of what we discussed in class. They should not be used as a replacement for attending lectures or reading books. Please let me know if you find any errors or have any questions about the material.

TipHypothesis

You can use the Hypothesis commenting system available on the right side panel (you can login with your U/Google account). This allows you annotate the notes (private) and leave comments or reply to other people’s comments (public).

Note

I encourage you to download the lecture notes and play around with the Julia code.

Download (.ipynb) Ask a Question ChimeIn (opens in a new tab)

Warning

Notes for this lecture have not been updated since we met (stay tuned for updates)

Note

This lecture is based on Chapter 1 of [Fundamentals of Numerical Computation by Driscoll, Braun] which is avaliable online (opens in a new tab)

2.1 Recap

Last time: Admin, Notebooks (Markdown, \(\LaTeX\), Julia), Errors, Significant digits, rate of convergence, “seeing” the rate on graphs.

Exercise 2.1

In England, when we need to convert between lb and kg some people use the phrase “2.2 pounds of jam is the same as a kilogram”.

• What is the error, absolute error, and relative error of approximating 1 kg with 2.2 lb?
• To how many significant digits is this approximation correct to?

You may use 1 kg = 2.20462262185 lb.

The error is simply “exact minus approx.” \(= 2.20462262185 - 2.2 \approx 0.00462262185\). This is also the absolute error (because it is positive).

For the relative error, we divide by the exact value:

x = 2.20462262185 
y = 2.2 

println("relative error = $( (x-y)/x )")

relative error = 0.0020967860005540483

The error is \(0.002096... \leq 5 \times 10^{-3}\) and so \(2.2\) is correct to 3 significant digits. (This is sufficient for all recipes you will likely use).

Exercise 2.2

Most people say a marathon is 26 miles when it is actually 26 miles 385 yards. In Europe, people say a marathon is 42 km when it is actually 42.194988 km.

• Which approximation has a smaller error?
• What is the advantage of computing the relative errors in this case?

Hint: You may use the fact that 1 mile \(=\) 1760 yards.

A marathon is \(26 \times 1760 + 385\) yards. The error and relative error in this approximation is

x = 26*1760 + 385
y = 26*1760

err, r_err = x-y, (x-y)/x

(385, 0.00834326579261025)

The same calculation for the approximation in km is

x′ = 42.194988
y′ = 42 

err′, r_err′ = x′-y′, (x′-y′)/x′

(0.19498800000000216, 0.004621117560218341)

r_err / r_err′

1.8054649516893149

The approximation in miles has a \(1.8\times\) larger relative error than the approximation in km.

Why is the error in the approximation in miles much larger than \(1.8\times\) bigger than the error in the approximation in km?

2.2 Standard/scientific form

Example 2.1 (Planck’s constant) Planck’s constant is \(6.62607015 × 10^{-34} \text{m}^2 \text{kg s}^{-1}\).

Exercise 2.3 (Astronomical Unit)

The astronomical unit is defined as \(1 \text{AU} = 149,597,870,700 \text{m}\). Write this in standard form.

\(1.495978707 \times 10^{11} \text{m}\)

2.3 Floating point numbers

There are many different versions of floating point numbers; we will pick one to give general ideas rather than specific details.

For fixed \(\ell, k\), we define

\[ \begin{align} \mathbb F_{k\ell} := &\Bigg\{ (-1)^s (1 + f) 2^n : f = (0.b_1b_2....b_k)_2, \nonumber\\ &\qquad\qquad n = c - (2^{\ell-1}-1), \quad c \in \{1,\dots,2^{\ell}-2\}\nonumber \Bigg\} \nonumber\\ % &\cup \{\pm0, \pm\infty, \text{NaN}\} \cup \{ \text{subnormal numbers} \} \nonumber \end{align} \] {eq-floats}

Here,

• \(k\) is the binary precision,
  
• \(s\in \{0,1\}\) is the sign bit, so \((-1)^s\) is the sign,
  
• \(1 + f\) is the significand,
  
• \(f\) is the mantissa, also called the fraction,
  
• \(n\) is the exponent, belongs to \(\{ 2 - 2^{\ell-1}, ...., 2^{\ell-1}-1 \}\),
  
• \(c\) is the characteristic,
  
• \(2^{\ell-1}-1\) is the bias

This is stored in a computer as the following \(\ell + k + 1\) bits:

\[\begin{align} \texttt{bitstring} \quad = \quad s \quad | \quad c \text{ written in binary} \quad | \quad b_1b_2....b_k \nonumber \end{align}\]

Special values.

When \(c = 0\), \(\texttt{bitstring}\) represents the number \((-1)^s f \cdot 2^{-(2^{\ell-1} - 1)}\). When \(f = 0\), this encodes \(\pm 0\) and when \(f \not=0\) these are called subnormal numbers,

When \(c = 2^{\ell}-1\), \(\texttt{bitstring}\) represents \(\pm \infty\) when \(f = 0\) and \(\text{NaN}\) when \(f\not=0\).

Exercise 2.4

Write down all the floating point numbers with \(k = 2\) in the interval \([1,2)\). What about the interval \([2,4)\)? \([2^n,2^{n+1})\)?

There are \(2^k\) equally spaced floating point numbers in any interval \([2^n,2^{n+1})\) (as long as \(\ell\) is sufficiently large).

Floats between \(1\) and \(2\) are \((1.b_1b_2)_2 = 1 + b_1 2^{-1} + b_2 2^{-2}\) and so are \(1\) (\(b_1 = b_2 = 0\)), \(1.25\) (\(b_1 = 0, b_2 = 1\)), \(1.5\) (\(b_1 = 1, b_2 = 0\)), \(1.75\) (\(b_1 = b_2 = 1\)).

Similarly, floats in the interval \([2^n, 2^{n+1})\) are of the form \((1 + b_1 2^{-1} + b_2 2^{-2} ) 2^n\). That is, \(2^n\) (\(b_1 = b_2 = 0\)), \(1.25 \times 2^n\) (\(b_1 = 0, b_2 = 1\)), \(1.5 \times 2^n\) (\(b_1 = 1, b_2 = 0\)), \(1.75 \times 2^n\) (\(b_1 = b_2 = 1\)).

Maximum number is \(2^{2^{\ell-1} -1}( 2 - 2^{-k} ) \approx 2^{2^{\ell-1}}\). If the result of a computation exceeds this value, it is stored as \(+\infty\) and we have an OVERFLOW error,

Minimum positive (normal) number is \(2^{2-2^{\ell-1}}\). If the result of a computation is between this value and zero, it is stored as \(0\) and we have an UNDERFLOW error,

Machine precision is \(\epsilon_{\text{mach}} = 2^{-k}\). The next largest float after \(1\) is \(1 + \epsilon_{\text{mach}}\).

Exercise 2.5

Take \(k=\ell=2\).

• What number does \(0 0 1 1 0\) encode?
• Write the bitstring of \(-3.5\).
• What are the maximum and minimum positive (normal) numbers and the machine precision in this system?

Recall that floats are stored as \(s | c | f\) where \(s\) is the sign bit, \(c\) the characteristic, and \(f\) the fraction or mantissa. Recall that the bias is \(b = 2^{\ell-1}-1 = 2 - 1 = 1\) and \(c \in \{1, \dots, 2^\ell-2 \} = \{1,2\}\).

In this case, we have \(s = 0\), \(c = (01)_2 = 1\) and \(f = (0.10)_2 = \tfrac12\) and so this number represents \((-1)^s (1 + f) 2^{c - b} = (-1)^0 (1 + \tfrac12) 2^{1-1} = \tfrac32\).

\(-3.5 = (-1)^1 (1 + 0.75) 2 = (-1)^1 (1 + (0.11)_2 ) 2^{2-1}\) and so the bitstring is \(1 1 0 1 1\).

The maximum is \((-1)^0 (1 + (0.11)_2 ) 2^{ 2 - 1} = \tfrac72\),

The minimum positive float is \((-1)^0 (1 + (0.00)_2) 2^{ 1 - 1 } = 1\),

Machine precision is \(2^{-k} = \tfrac14\).

In the following, we will write \(1|\ell|k\) for the number of bits used to represent the sign (\(1\)), the exponent (\(\ell\)), and the significand (\(k\)).

In practice, we have:

• Float16 = half precision = \(1|5|10\),
  
• Float32 = single precision = \(1|8|23\), and
  
• Float64 = double precision = \(1|11|52\),

Exercise 2.6

What are the maximum and minimum positive numbers and the machine precision for each of Float16, Float32, and Float64?

Recall, the bitstring is stored as \(s | c | f\) with \(1 | \ell | k\) bits for each where

• Float16 = half precision = \(1|5|10\),
  
• Float32 = single precision = \(1|8|23\), and
  
• Float64 = double precision = \(1|11|52\).

In general, the maximum number is \(2^{2^{\ell-1} -1}( 2 - 2^{-k} ) \approx 2^{2^{\ell-1}}\), the minimum positive (normal) number is \(2^{2-2^{\ell-1}}\), and the machine precision is \(\epsilon_{\text{mach}} = 2^{-k}\). Therefore:

• Float16 (\(\ell = 5\), \(k = 10\)): the maximum is \(2^{15}(2 - 2^{-10}) = 65504\), the minimum positive (normal) number is \(2^{-14} \approx 6.10 \times 10^{-5}\), and the machine precision is \(2^{-10} \approx 9.77 \times 10^{-4}\).
• Float32 (\(\ell = 8\), \(k = 23\)): the maximum is \(2^{127}(2 - 2^{-23}) \approx 3.40 \times 10^{38}\), the minimum positive (normal) number is \(2^{-126} \approx 1.18 \times 10^{-38}\), and the machine precision is \(2^{-23} \approx 1.19 \times 10^{-7}\).
• Float64 (\(\ell = 11\), \(k = 52\)): the maximum is \(2^{1023}(2 - 2^{-52}) \approx 1.80 \times 10^{308}\), the minimum positive (normal) number is \(2^{-1022} \approx 2.23 \times 10^{-308}\), and the machine precision is \(2^{-52} \approx 2.22 \times 10^{-16}\).

You can check these in Julia with floatmax, floatmin, and eps:

for T ∈ (Float16, Float32, Float64)
    label = "$T:  "
    pad = " "^length(label)
    println(label, "max = ", floatmax(T), "\n", pad, "min = ", floatmin(T), "\n", pad, "eps = ", eps(T))
end

Float16:  max = 6.55e4
          min = 6.104e-5
          eps = 0.000977
Float32:  max = 3.4028235e38
          min = 1.1754944e-38
          eps = 1.1920929e-7
Float64:  max = 1.7976931348623157e308
          min = 2.2250738585072014e-308
          eps = 2.220446049250313e-16

2.4 Floating point arithmetic

We now consider a set of floating point numbers \(\mathbb F := \mathbb F_{k,\ell}\) and define a function \(\mathrm{fl} : \mathbb R \to \mathbb F\) that maps real numbers to (one of) the closest points in \(\mathbb F\).

Theorem 2.1 If \(x \not= 0\) and \(|x| \leq (2-2^{-k}) 2^{2^{\ell-1}-1}\), then

\[\begin{align} \frac {| x - \mathrm{fl}(x) |} {|x|} \leq \frac12 \epsilon_{\rm mach}. \end{align}\]

Proof. Suppose that \(x \in [2^n, 2^{n+1})\). The floats in the interval \([2^n, 2^{n+1})\) are of the form \((1 + z \cdot 2^{-k}) \cdot 2^n\) for \(z = 0,\dots,2^{k}-1\). Therefore, consectutive floats are spaced apart by \(2^{n-k}\) and since \(|x| \geq 2^n\), we have

\[\begin{align} \frac {| x - \mathrm{fl}(x) |} {|x|} \leq \frac12 \frac {2^{n-k}} {|x|} \leq \frac12 2^{-k} = \frac12 \epsilon_{\rm mach} \end{align}\]

2.5 Float16

In Julia you can specify the precision by using FloatN() where \(N = 16, 32,\) or \(64\). e.g.

x = Float16( π )
x == 3.140625

true

Exercise 2.7

Write down the bitstring for \(\pi\) in Float16.

\(\pi\) in Float16 rounds to

x = Float16( π )
x == 3.140625

true

and, in Float16, we have \(1|\ell|k = 1|5|10\). Therefore, the bias is \(2^{\ell-1}-1 = 2^{4}-1 = 15\)

Since \(3.140625 = (-1)^0 (1 + 1/2 + 0/4 + 0/8 + 1/16 + 0/32 + 0/64 + 1/128)\times 2^{16-15}\), the bitstring is \(0 | 1 0 0 0 0 | 1 0 0 1 0 0 1 0 0 0\)

bitstring(x)

"0100001001001000"

You can see how this is represented in binary using bitstring():

bs = bitstring( x )

"0100001001001000"

Here, the first bit is the sign and the next 5 bits correspond to the characteristic:

s = parse(Int, bs[1:1], base=2);
c = parse(Int, bs[2:6], base=2);
@show ( s, c );

(s, c) = (0, 16)

The final bits bs[7:16] \(:= b_1b_2 ... b_{10}\) represent the significand/mantissa in the following way

\[\begin{align} f &= \big(0 . b_1 b_2 b_3 b_4 b_5 b_6 b_7 b_8 b_9 b_{10}\big)_2 \nonumber\\ % &= \sum_{i=1}^{10} b_i 2^{-i} \end{align}\]

That is, bs[7:16] is \(2^{10} \times\) larger than \(f\):

f = 2^(-10) * parse(Int, bs[7:16], base=2)

0.5703125

# exponent is c - 15 (characteristic minus bias)
# significand is 1 + f
[ exponent(x), significand(x) ]

2-element Vector{Float16}:
 1.0
 1.57

y = (-1)^s * 2.0^( c - 15 ) * ( 1 + f )

3.140625

Recall that x = Float16( pi )

x == y

true

2.6 Errors in Float16

using Plots 
using LaTeXStrings 

ϵ = eps( Float16 );
x = range( -2.0^(16), 2.0^(16), 10000 ); 
y = @. ( ( x - Float16(x) ) );

plot( x , y,  seriestype=:scatter, ms=1, title=L"$x - \mathrm{fl}(x)$ in Float16", legend=false)
vline!( [2.0^n for n in 12:1:16] )

ϵ = eps( Float16 );
x = range( -2^(16), 2^(16), 10000 ); 
y = @. ( ( x - Float16(x) ) /x );

scatter( x , y; 
    ms=1, title=L"$\frac{x - \mathrm{fl}(x)}{x}$ in Float16", legend=false)
vline!( [2.0^n for n in 13:1:16] )
hline!( [-ϵ/2, ϵ/2] )

Special values:

a = Float16( 0 );
(a, bitstring(a) )

(Float16(0.0), "0000000000000000")

b = -Float16( 0 );
(b, bitstring(b) )

(Float16(-0.0), "1000000000000000")

c = Float16( 1/0 );
(c, bitstring(c) )

(Inf16, "0111110000000000")

c = Float16( -1/0 );
(c, bitstring(c) )

(-Inf16, "1111110000000000")

c = Float16( 0/0 );
(c, bitstring(c) )

(NaN16, "1111111000000000")

By default, Julia uses Float64. We will use this datatype for the remainder of this course.

typeof( 1.0 + π )

Float64

2.7 Subtractive cancellation and conditioning

when adding or subtracting numbers so that the result is much smaller in magnitude. For example, we know \(x = -1.0000000000000000000147\) to \(23\) significant digits but \(x + 1 = 1.47 \times 10^{-22}\) to only \(3\) significant digits.

This concept can be illustrated using (relative) condition number: suppose we are interested in a quantity \(f(x)\) for \(x\in \mathbb R\). In practice, we can only store \(\widetilde{x} := \mathrm{fl}(x)\) in a computer (recall that \(\widetilde{x} = x(1 + \epsilon)\) for some \(|\epsilon|\leq \frac12 \epsilon_{\mathrm{mach}}\)). Then the ratio of the relative error in the output to the relative error in the input in the limit as the precision goes to zero is known as the (relative) condition number \(\kappa_f(x)\):

\[\begin{align} \frac {\left|\frac {f(x) - f(\widetilde{x})} {f(x)} \right|} {\left|\frac {x-\widetilde{x}} {x}\right|} % &= \left| \frac {f(x) - f(x + \epsilon x)} {\epsilon f(x)} \right| \nonumber\\ % &\to \kappa_f(x) \quad \text{as } \epsilon \to 0 \end{align}\]

That is, the (relative) condition number is approximately the constant of proportionality between the relative error in the output and the relative error of the input:

\[\begin{align} \left|\frac {f(x) - f(\widetilde{x})} {f(x)} \right| \approx \kappa_f(x) |\epsilon| \end{align}\]

Example 2.2 The condition number of \(f(x) = x + 1\) is \(\kappa_f(x) = \left|\frac{x}{x+1}\right|\). This blows up when \(x\approx -1\) which explains subtractive cancellation that we mentioned above.

If the condition number is large, we say the problem is ill-conditioned or the problem is badly conditioned.

Notice that

\[ \begin{align} \kappa_f(x) &= \lim_{\epsilon\to0}\left| \frac {f(x) - f(x + \epsilon x)} {\epsilon f(x)} \right| \nonumber\\ % &= \lim_{\epsilon\to0}\left| \frac {f(x) - f(x + \epsilon x)} {\epsilon x} \right| \left| \frac{x}{f(x)} \right| \nonumber\\ % &= \left| \frac{x f'(x)}{f(x)} \right| \end{align} \tag{2.1}\]

TipTO-DO

• Assignment 1: Due Wednesday (Install Julia),
• Exercises throughout the notes and below,
• (Sign-up to office hours if you like)

2.8 Extra exercises

Exercise 2.8

Compute the relative condition numbers of (i) \(x \mapsto cx\) for \(c \not= 0\), (ii) \(\log(x)\), (iii) \(\exp\), (iv) \(x^p\), and (v) \(f\circ g\) in terms of \(\kappa_f\) and \(\kappa_g\)

Recall that \(\kappa_f(x) = \left| \frac{x f'(x)}{f(x)} \right|\).

Since \((x \mapsto cx)' = c\), \(\log'(x) = \tfrac1x\), \(\exp'(x) = \exp(x)\), \((x^p)' = p x^{p-1}\), and \((f \circ g)'(x) = f'(g(x)) g'(x)\), we have

\[\begin{align} \kappa_{x \mapsto c x}(x) &= \left| \frac{x c}{cx} \right| = 1 \nonumber\\ % \kappa_{\log}(x) &= \left| \frac{x \cdot \tfrac1x}{\log(x)} \right| = \frac{1}{|\log x|} \qquad \text{for } x \not= 1, \nonumber\\ % \kappa_{\exp}(x) &= \left| \frac{x e^x}{e^x} \right| = |x| \nonumber\\ % \kappa_{x\mapsto x^p}(x) &= \left| \frac{x \cdot p x^{p-1}}{x^p} \right| = p \nonumber \end{align}\]

and

\[\begin{align} \kappa_{f\circ g}(x) &= \left| \frac{x f'(g(x)) g'(x)}{f(g(x))} \right| \nonumber\\ % &=\left| \frac{g(x) f'(g(x))}{ f(g(x))} \right| \left| \frac{x g'(x)}{ g(x) } \right| \nonumber\\ % &= \kappa_{f}(g(x)) \kappa_g(x) \nonumber \end{align}\]

Exercise 2.9

Midterm, Fall 2025, Compulsory Section, 3 pts. Let us define

\[\begin{align*} f(x) = \frac{\sqrt{1+x}-1}{x}, \qquad g(x) = \frac{1}{\sqrt{1+x} + 1}. \end{align*}\]

One may verify that \(f\) and \(g\) are mathematically identical (but you don’t have to). The following table is the output of evaluating \(f\) and \(g\) on \(\{10^{-14}, 10^{-15}, 10^{-16}\}\) in double precision floating point. Why are the function values different?

f = x -> ( sqrt(1 + x) - 1 ) / x
g = x-> 1/( 1 + sqrt(1 + x) )

t = [1e-14; 1e-15; 1e-16];
   
println( "                         f(x)                  g(x)" )
println( "x = 10^(-14):     ", f(t[1]), "    ",    g(t[1]) )
println( "x = 10^(-15):     ", f(t[2]), "   ",    g(t[2]) )
println( "x = 10^(-16):     ", f(t[3]), "                   ",    g(t[3]) )

                         f(x)                  g(x)
x = 10^(-14):     0.4884981308350689    0.4999999999999988
x = 10^(-15):     0.44408920985006256   0.4999999999999999
x = 10^(-16):     0.0                   0.5

Subtractive cancellation! The formula \(f\) suffers from subtractive cancellation, but the formula \(g\) does not.

(The formulae are the same since

\[\begin{align} f(x) &= \frac{\sqrt{1+x}-1}{x} \frac{\sqrt{1+x}+1}{\sqrt{1+x}+1} \nonumber\\ % &= \frac{1+x-1}{\sqrt{1+x}+1} = g(x). \nonumber \end{align}\]

Moreover, we can see that \(\lim_{x\to 0} f(x) = \lim_{x\to 0} g(x) = \frac12\).)