# M13. Machine Epsilon

One does not so much make a case *for* floating point as require what it so generously offers.

Recall that an $n$-bit machine is capable of representing $2^n$ discrete states, each of which may be assigned a unique integer $z\in \mathbb{Z}$, usually via unsigned binary or two's complement. Whatever the specifics, it is possible to construct a bijective map between the set of possible machine states and an interval $I \subset \mathbb{Z}$ of cardinality $|I|=2^n$.

For sufficiently large $n$, we obtain exact representations of most of the integers we might encounter in everyday life. The range of *uint32*, for instance, is 0 to $2^{32} =$ 4,294,967,296, while the range of *int32* is $-2^{31}=$ 2,147,483,648 to $+2^{31} -1 =$ +2,147,483,647.

We note that this approach is by no means the last word.

Most inconveniently, the Universe is under no compulsion to actually respect these ranges. Some numbers are simply Very Big:

- The observable universe is roughly $10^{27}$ meters across.
- There are around $10^{24}$ particles in a mole.
- Our species consumes roughly $10^{13}$ watts of power *per annum*. 

To compound our troubles, many interesting computations occur not even over $\mathbb{Z}$ but over $\mathbb{R}$ (the entirety of calculus, for instance, as well as continuous-time signal processing). The incorporation of the reals poses a challenge that is in a certain sense insurmountable. We note that while $\mathbb{Z}$ and $\mathbb{R}$ are both infinite sets, there *do* exist finite intervals of $\mathbb{Z}$ which may be equinumerous with our machine's state space, whatever size it is. On the other hand, in general, real intervals are not even countable, much less finite.

A classic discrete-math exercise is to show that every open interval on $\mathbb{R}$ has the same cardinality as $\mathbb{R}$ itself. One approach is to prove that the map $f(x): \mathbb{R}\to (0,1)$, $x \mapsto \frac{1}{1+e^{-x}}$, also known as the *logistic sigmoid*, is a bijection (in turn, one can easily construct a bijection from $(0,1)$ to any interval $(a,b)$ as a matter of scaling and shifting). 

The uncountability of $\mathbb{R}$ (and therefore any open interval of $\mathbb{R}$) is then shown by deploying Cantor's *diagonal argument* on $(0,1)$:

For sake of contradiction, suppose that $(0,1)$ *is* enumerable. By assumption, we may list *all* of the elements in this set by some scheme or other

|      LIST      |
|:-------------:|
| 0.48936032... |
| 0.60012377... |
| 0.23158392... |
| 0.78832125... |
| 0.34829373... |
| 0.11239586... |
|      ...      |

Now consider a number constructed by joining each of the digits in the main diagonal of such a list

|      LIST      |
|:-------------:|
| 0.**4**8936032... |
| 0.6**0**012377... |
| 0.23**1**58392... |
| 0.788**3**2125... |
| 0.3482**9**373... |
| 0.11239**5**86... |
|      ...      |

Then, change each digit of this number to *another* digit (strictly speaking, another digit that isn't 9). We have constructed a number that *cannot* have been enumerated (every listed number differs from the newly-constructed number by at least one digit). Contradiction&mdash;the reals must be uncountable.  

It is thus fundamentally impossible to maintain exact representations of even a 'tiny' segment of the real line on a digital computer. This is good to know. Aspiring to achieve the impossible is, in my experience, demoralizing.

What is instead attempted is the *approximation* of $\mathbb{R}$, a sort of illusion which we have all occasionally broken: 


```python
print(1.2-1.)
```

    0.19999999999999996


By *floating point*, we refer to an approximate representation of the real numbers that also solves, at a stroke, much of the scaling problem discussed prior (as in the representation extends to both very large and very small numbers). The scheme is sketched as follows:

- Define an even base $\beta$, also known as a *radix*. 
- Define a *precision* $p$.
- A *floating point number* is then represented by $$\pm (d_0 + d_1\beta^{-1} + d_2\beta^{-2}+\ldots+d_{p-1}\beta^{-(p-1)})\beta^e$$
where $e$, the *exponent*, is bounded by the machine, and $0 \leq d_i < \beta$. Additionally, one may impose the constraint $d_0 \not = 0$, also called *normalization*, to ensure the uniqueness of all representations. The number $d_0.d_1 d_2 \ldots d_{p-1}$ is often called the *mantissa* of the number (sometimes, *significand*).

We will not discuss the actual encodings by which floating point is usually implemented (namely, the standards set forth in IEEE 754). What is more important is the bounding of error (and thus the bounding of our ignorance).

To start, note that for a fixed exponent $e$, discrete floating point representations are available every $\beta^{e-(p-1)}$. This is clear from decrementing the unit in the last place (ulp) of the mantissa by 1: 

\begin{align}
(d_0 + d_1\beta^{-1} + d_2\beta^{-2}+\ldots+d_{p-1}\beta^{-(p-1)})\beta^e-(d_0 + d_1\beta^{-1} + d_2\beta^{-2}+\ldots+(d_{p-1}-1)\beta^{-(p-1)})\beta^e &= d_{p-1}\beta^{-(p-1)}\beta^{e}-(d_{p-1}-1)\beta^{-(p-1)}\beta^{e} \\
&=(d_{p-1}-(d_{p-1}-1))\beta^{-(p-1)}\beta^{e}\\
&=\beta^{e-(p-1)}.
\end{align}

Thus, the floating point representation of a real number may differ from that real number by as much as $\frac{\beta^{e-(p-1)}}{2}$ assuming that the representation is obtained through conventional rounding. Normalizing this difference by $\beta^{e}$ gives a notion of *maximum relative error*:

$$\frac{\frac{\beta^{e-(p-1)}}{2}}{\beta^e} = \frac{\beta^{-(p-1)}}{2}.$$

This number, subject to the parameters $\beta$ and $p$, is known as the *machine epsilon* and denoted as $\varepsilon$. Simply put, it is an upper bound on the relative error of a floating point representation $r'$ of a real number $r$

$$\frac{r'-r}{r} \leq \varepsilon.$$

Somewhat frustratingly (actually, extremely frustratingly), we often see the machine epsilon defined alternatively as the smallest number $\varepsilon$ such that $1+\varepsilon \not = 1$ on a given machine (this is equivalent to the distance to the next floating point representation after 1.).

If this should be desired, it can be easily approximated by something like this: 


```python
def machineEpsilon():
    epsilon = 1.
    while(1. + epsilon != 1.):
        epsilon = epsilon/2.
    return epsilon*2 # roll back one iteration

machineEpsilon()
```




    2.220446049250313e-16



The result is in agreement with the parameter provided by the system module:


```python
import sys

sys.float_info.epsilon
```




    2.220446049250313e-16



266244f533ea44152d5b0366eb09ab5090f30bd8
ipynb
Jupyter Notebook
M13_Machine_Epsilon.ipynb
brekekekex/computer_organization_memoranda
d16dc251075d3da49aaf01f148676f857d05dc4b
Unlicense
2020-01-17T16:34:17.000Z
2020-02-23T22:06:07.000Z
M13_Machine_Epsilon.ipynb
brekekekex/computer_organization_memoranda
d16dc251075d3da49aaf01f148676f857d05dc4b
Unlicense
M13_Machine_Epsilon.ipynb
brekekekex/computer_organization_memoranda
d16dc251075d3da49aaf01f148676f857d05dc4b
Unlicense
Qwen/Qwen-72B
1. YES
2. YES