# Finite precision point arithmetic

## Topics in this lab

## Introduction

We all know that when we use the value of pi in
Matlab, or compute cos(3.4), we are not getting
the "exact" value of $\pi$ or the cosine function, but rather an
approximation. In fact, we don't expect to be able to compute any
irrational function to all of its digits, if only because we know that
such values are nonterminating, non-repeating decimals values. In
typical practical situations, our answers will not be significantly
affected by the approximations that are made. What may be suprising,
however, is that even typical "rational" numbers, such as 0.1 are only
approximately represented on the computer.

In this lab, we will explore the number system represented by floating
point arithmetic, and discuss some of the consequences for scientific
computing. The ideas presented here extend to most modern computing
systems, not just Matlab.

clear all format long e

Back to the top

## Examples : Floating point arithmetic

### Example 1

Compute the following.

$$ x = 0.1 + 0.1 + 0.1 + ... + 0.1 \qquad \mbox{(10 times)} $$x = 0; for i = 1:10, x = x + 0.1; end

What is the difference between $x$ and 1?

`fprintf('%g\n',abs(x-1))`

1.11022e-16

### Example 2

Compute $x$.

$$ x = 2 - 3\left(\frac{4}{3} - 1\right) $$```
% Algebra tells us that this should be 1
x = 2 - 3*(4/3 - 1);
```

What is the difference between $x$ and 1?

`fprintf('%g\n',abs(x-1));`

2.22045e-16

### Example 3

Verify the following exact mathematical expression for the value $a = 0.3$.

$$ 1 + a + a^2 + a^3 + a^5 = \frac{1-a^6}{1-a} $$a = 0.3; S_left = sum(a.^(0:5)); S_right = (1-a^6)/(1-a);

How close are the left and right sides of this expression?

`fprintf('%g\n',abs(S_left - S_right));`

2.22045e-16

### Example 4

For the function
$f(x) = x$, we can
compute the derivative $f'(x)$
*exactly* using the formula
$$
\begin{eqnarray*}
f'(x) & = & \frac{f(x+h) - f(x)}{h} = \frac{(x+h) - x}{h} = 1
\end{eqnarray*}
$$

Verify this formula for $x = 1$ and $h = 10^{-3}$.

% Compute the derivative of f(x) = x using the secant method x = 1; h = 1e-3; dfdx = ((x + h) - x)/h; % This should be exactly 1

Do we get $f'(x) = 1$?

`fprintf('%g\n',abs(dfdx-1));`

1.10134e-13

Back to the top

## More examples : Is floating point arithmetic commutative and associative?

You may recall from your first introduction to algebra that we can
often re-arrange the order of operations like addition or multiplication.
This *communative* property of these operations allows us to write
$$
a + b + c = a + c + b = b + c + a
$$

and so on. Unfortunately, this is not always true for floating
point arithmetic.

The associative property that we learned in algebra, i.e. that
$$
a + (b + c) = (a + c) + b
$$

may not always hold either.

### Example 5

The order in which we add numbers can matter. Consider these two expressions $$ x = 10^{16} + 1 - 10^{16} $$ $$ y = 10^{16} - 10^{16} + 1 $$ $$ z = 10^{16} - (10^{16} - 1) $$

```
x = 1e16 + 1 - 1e16;
y = 1e16 - 1e16 + 1;
z = 1e16 - (1e16 - 1); % Test the associative property
```

Are $x$, $y$ and $y$ equal?

fprintf('x = %g\n',x); fprintf('y = %g\n',y); fprintf('z = %g\n',z);

x = 0 y = 1 z = 0

### Example 6

The above can happen for much smaller values as well

x = 1 + 0.1 - 1; y = 1 - 1 + 0.1; z = 1 - (1 - 0.1);

Comparing the differences of these three values, we get

fprintf('x-y = %g\n',x-y); fprintf('y-z = %g\n',y-z); fprintf('x-z = %g\n',x-z);

x-y = 8.32667e-17 y-z = 2.77556e-17 x-z = 1.11022e-16

### Example 7

In this example, we add up 100 random numbers in different orders.

rand('seed',1110); % Get the same random numbers each time x = rand(100,1); xperm = x(randperm(100)); % Permute the values in array x sum_orig = sum(x); sum_perm = sum(xperm);

Compare the sum of the original array and the permuted array.

`fprintf('%g\n',abs(sum_orig-sum_perm));`

1.42109e-14

Back to the top

## Examples (modified)

In some of the examples above, a slight change of the constants involved can fix the problems that we saw above. Below are modified versions of some of the above examples.

### Example 1 (modified)

Compute the following.

$$ x = 0.125 + 0.125 + 0.125 + ... + 0.125 \qquad \mbox{(8 times)} $$x = 0; for i = 1:8, x = x + 0.125; end

What is the difference between $x$ and 1?

`fprintf('%g\n',abs(x-1))`

0

### Example 2 (modified)

Compute $x$.

$$ x = 2 - 2\left(\frac{3}{2} - 1\right) $$```
% Algebra tells us that this should be 1
x = 2 - 2*(3/2 - 1);
```

What is the difference between $x$ and 1?

`fprintf('%g\n',abs(x-1));`

0

### Example 3 (modified)

Verify the following exact mathematical expression for the value $a = 0.0625$.

$$ 1 + a + a^2 + a^3 + a^5 = \frac{1-a^6}{1-a} $$a = 0.0625; S_left = sum(a.^(0:5)); S_right = (1-a^6)/(1-a);

How close are the left and right sides of this expression?

`fprintf('%g\n',abs(S_left - S_right));`

0

### Example 4 (modified)

For the function $f(x) = x$, verify the following formula for $x = 1$ and $h = 0.015625$. $$ \begin{eqnarray*} f'(x) & = & \frac{f(x+h) - f(x)}{h} = \frac{(x+h) - x}{h} = 1 \end{eqnarray*} $$

```
x = 1;
h = 0.015625;
dfdx = ((x + h) - x)/h; % This should be exactly 1
```

Do we get $f'(x) = 1$?

`fprintf('%g\n',abs(dfdx-1));`

0

Back to the top

## Machine epsilon

As as way to gauge how close two numbers are to each other, we can use the Matlab function eps(x). This measures the distance to the next representable number. For example,

### Example 7

In this example, we want to see if the following mathematical expression holds for all $\x > 0$.

Compute $x$.

$$ 1 + x > 1, \qquad x > 0 $$x = 1; while (1 + x > 1) x = x/2; end

Is $x$ equal to 0?

`fprintf('x = %8.4e \n',x);`

x = 1.1102e-16

Back to the top

## Lab exercises

- Can we modify examples 5,6 and 7 above to fix the apparent problems?
- Using a loop similar to the one above to compute the smallest value we could add to 1 and still have something greater than 1, try to find the smallest value we can add to $10^6$ and still get something larger than a million?

Compare your answers with the solutions.

Back to the top

## Get the code

Do you want to try the above code fragments on your own? Download the Matlab script that produces this page here. (lab_19.m)

Published with MATLAB® 8.5