<jupyter_start><jupyter_text>INF-285 - Computación Científica      Floating Point Arithmetic       [S]cientific [C]omputing [T]eam       Version: 1.20   Table of Contents* [Introduction](intro)* [The nature of floating point numbers](nature)* [Visualization of floating point numbers](visualization)* [What is the first integer that is not representable in double precision?](firstinteger)* [Loss of significance](loss)* [Loss of significance in funcion evaluation](func)* [Another analysis (example from textbook)](another)* [Acknowledgements](acknowledgements)<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline<jupyter_output><empty_output><jupyter_text>Introduction[Back to TOC](toc)Hello! This notebook is an introduction to how our computers handle the representation of real numbers using double-precision floating-point format. To understand the contents of this notebook you should have at least a basic notion of how binary numbers work.The double-precision floating-point format occupies 64 bits which are divided as follows:* 1 bit for the sign* 11 bits for the exponent* 52 bits for the mantissaThis means that the very next representable number after $1$ is $1 + 2^{-52}$, and their difference, $2^{-52}$, is the called $\epsilon_{mach}$.Additionally, if you'd like to quickly go from a base-2 integer to a base-10 integer and viceversa, Python has some functions that can help you with that.<jupyter_code>int('0b11', 2)
bin(9)
bin(2**53)<jupyter_output><empty_output><jupyter_text>The nature of floating point numbers[Back to TOC](toc) As we know until now, float representations of real numbers are just a finite and bounded representation of them. But another interesting thing, is that these floating numbers are distributed across the real numbers. To see that, it's really important to keep in mind the following property of floating point numbers:$$\begin{equation} \left|\frac{\text{fl}(x)-x}{x}\right| \leq \frac{1}{2} \epsilon_{\text{mach}} \end{equation},$$where $\text{fl}(x)$ means the float representation of $x \in \mathbb{R}$. What it says is that **the relative error in representing any non-zero real number x, is bounded by a quantity that depends on the precision bign used**, i.e. ($\epsilon_{\text{mach}}$).Maybe now you're thinking: what does this relationship have to do with the distribution of floating point numbers? So, if we rewrite the previous as follows we get:$$\begin{equation} |\text{fl}(x)-x| \leq \frac{1}{2} \epsilon_{\text{mach}} |x| \end{equation}.$$It's clearer: **The absolute error (distance) between a real number and its floating point representation is proportional to the real number's magnitude.**Intuitively speaking, the representation error of a number increases as its magnitude increases, this implies that **the distance between a floating point number and the next representable floating point number will increase as the magnitude of such number increases (and conversely)**. Could you prove that?  For now we will prove it experimentally.We will use a library named **bitstring** to handle different number representations. You can install it with:```pip install bitstring```<jupyter_code>import bitstring as bs<jupyter_output><empty_output><jupyter_text>The next two functions are self-explanatory:1. `next_float(f)` computes the next representable float number right after $f$.2. `gap(f)` computes the difference between $f$ and the next representable float number.Notice that here we are considering that we are using double-precision.<jupyter_code>def next_float(f):
    #packing double-precision foat
    b = bs.pack('>d', f)
    
    #extracting mantisa as unsigned int and adding up 1.
    # There are two cases, (1) if the bits of the mantissa are all 1,
    # (2) all the other cases.
    m = b[12:].uint
    if m==4503599627370495:
        # Case (1)
        m=0
        b[12:] = m
        exp=b[1:12].uint
        exp +=1
        b[1:12] = exp
        
    else:
        # Case (2)
        m += 1
        # putting the result in his place
        b[12:] = m
    
    return b.float

def gap(f):
    next_f = next_float(f)
    return next_f - f<jupyter_output><empty_output><jupyter_text>So if we compute `gap(1)` we should get the epsilon machine number. Let's try it:<jupyter_code>gap(1)
# This is the case (1) mentioned in the definition of next_float.
next_float(2-gap(1))<jupyter_output><empty_output><jupyter_text>In order to prove our hypotesis, we will create an array of values: $[2^{-5},...,2^{60}]$ and compute their corresponding gaps.<jupyter_code>values = np.array([2**i for i in range(-5,60)]).astype(float)

# Corresponding gaps:
# The Numpy function "vectorize" is very useful to be use one wants to apply
# a scalar function to each element of an array.
vgap = np.vectorize(gap)
gaps = vgap(values)<jupyter_output><empty_output><jupyter_text>We include now a comparison between a linear scale plot and a loglog scale plot. Which one is more useful here?<jupyter_code>fig = plt.figure(figsize=(10,5))

plt.subplot(121)
plt.plot(values, gaps,'.',markersize=20)
plt.xlabel('float(x)')
plt.ylabel('Gap between next representable number')
plt.title('Linear scale')
plt.grid(True)

plt.subplot(122)
plt.loglog(values, gaps,'.')
plt.xlabel('float(x)')
plt.ylabel('Gap between next representable number')
plt.title('Log scale')
plt.grid(True)

fig.tight_layout()
plt.show()<jupyter_output><empty_output><jupyter_text>As you can see, the hypotesis was right. In other words: Floating point numbers are not linearly distributed across the real numbers, and the distance between them is proportional to their magnitude. **Tiny numbers (~ 0) are closer between each other than larger numbers are.**Moreover, we can conclude that for large values the gap is larger that $1$, **which means that there will be integers that will not be stored**.  Visualization of floating point numbers[Back to TOC](toc)With the help of `bitstring` library we could write a function to visualize floating point numbers in his binary representation<jupyter_code>def to_binary(f):
    b = bs.pack('>d', f)
    b = b.bin
    #show sign + exponent + mantisa
    print(b[0]+' '+b[1:12]+ ' '+b[12:])<jupyter_output><empty_output><jupyter_text>Let's see some intereseting examples<jupyter_code>to_binary(1.)
int('0b01111111111', 2)
to_binary(1.+gap(1.))
to_binary(+0.)
to_binary(-0.)
to_binary(np.inf)
to_binary(-np.inf)
to_binary(np.nan)
to_binary(-np.nan)
to_binary(2.**-1074)
print(2.**-1074)
to_binary(2.**-1075)
print(2.**-1075)
to_binary(9.4)<jupyter_output>0 10000000010 0010110011001100110011001100110011001100110011001101<jupyter_text>What is the first integer that is not representable in double precision?[Back to TOC](toc) Recall that $\epsilon_{\text{mach}}=2^{-52}$ in double precision.<jupyter_code>to_binary(1)
to_binary(1+2**-52)<jupyter_output>0 01111111111 0000000000000000000000000000000000000000000000000000
0 01111111111 0000000000000000000000000000000000000000000000000001<jupyter_text>This means that if we want to store any number in the interval $[1,1+\epsilon_{\text{mach}}]$, only the numbers $1$ and $1+\epsilon_{\text{mach}}$ will be stored. For example, compare the exponent and the mantissa in the previous cell with the following outputs:<jupyter_code>for i in np.arange(1,11):
    to_binary(1+i*2**-55)<jupyter_output>0 01111111111 0000000000000000000000000000000000000000000000000000
0 01111111111 0000000000000000000000000000000000000000000000000000
0 01111111111 0000000000000000000000000000000000000000000000000000
0 01111111111 0000000000000000000000000000000000000000000000000000
0 01111111111 0000000000000000000000000000000000000000000000000001
0 01111111111 0000000000000000000000000000000000000000000000000001
0 01111111111 0000000000000000000000000000000000000000000000000001
0 01111111111 0000000000000000000000000000000000000000000000000001
0 01111111111 0000000000000000000000000000000000000000000000000001
0 01111111111 0000000000000000000000000000000000000000000000000001<jupyter_text>Now, we can scale this difference such that the scaling factor multiplied but $\epsilon_{\text{mach}}$ is one. The factor will be $2^{52}$. This means $2^{52}\,\epsilon_{\text{mach}}=1$. Repeating the same example as before but with the scaling factor we obtain:<jupyter_code>for i in np.arange(11):
    to_binary((1+i*2**-55)*2**52)<jupyter_output>0 10000110011 0000000000000000000000000000000000000000000000000000
0 10000110011 0000000000000000000000000000000000000000000000000000
0 10000110011 0000000000000000000000000000000000000000000000000000
0 10000110011 0000000000000000000000000000000000000000000000000000
0 10000110011 0000000000000000000000000000000000000000000000000000
0 10000110011 0000000000000000000000000000000000000000000000000001
0 10000110011 0000000000000000000000000000000000000000000000000001
0 10000110011 0000000000000000000000000000000000000000000000000001
0 10000110011 0000000000000000000000000000000000000000000000000001
0 10000110011 0000000000000000000000000000000000000000000000000001
0 10000110011 0000000000000000000000000000000000000000000000000001<jupyter_text>Which means we can only store exactly the numbers:<jupyter_code>to_binary(2**52)
to_binary(2**52+1)<jupyter_output>0 10000110011 0000000000000000000000000000000000000000000000000000
0 10000110011 0000000000000000000000000000000000000000000000000001<jupyter_text>This means, the distance now from $2^{52}$ and the following number representable is $1$ !!!! So, what would happend if I can to store $2^{53}+1$?<jupyter_code>to_binary(2**53)
to_binary(2**53+1)<jupyter_output>0 10000110100 0000000000000000000000000000000000000000000000000000
0 10000110100 0000000000000000000000000000000000000000000000000000<jupyter_text>I can't stored the **Integer** $2^{53}+1$! Thus, the first integer not representable is $2^{53}+1$.  Loss of significance[Back to TOC](toc) As we mentioned, there's a small gap between 1 and the next representable number, which means that if you want to represent a number between those two, you won't be able to do so; what you would need to do is to round it to a representable number before storing it in memory.<jupyter_code>a = 1.
b = 2.**(-52) #emach
result_1 = a + b     # arithmetic result is 1.0000000000000002220446049250313080847263336181640625
result_1b = result_1-1.0
print("{0:.1000}".format(result_1))
print(result_1b)
print(b)
c = 2.**(-53)
result_2 = a + c     # arithmetic result is 1.00000000000000011102230246251565404236316680908203125
np.set_printoptions(precision=16)
print("{0:.1000}".format(result_2))
print(result_2-a)
to_binary(result_2)
to_binary(result_2-a)
d = 2.**(-53) + 2.**(-54)

result_3 = a + d     # arithmetic result is 1.000000000000000166533453693773481063544750213623046875
print("{0:.1000}".format(result_3))
to_binary(result_3)
to_binary(d)<jupyter_output>1.0000000000000002220446049250313080847263336181640625
0 01111111111 0000000000000000000000000000000000000000000000000001
0 01111001010 1000000000000000000000000000000000000000000000000000<jupyter_text>As you can see, if you try to save a number between $1$ and $1 + \epsilon _{mach}$, it will have to be rounded (according to IEEE rounding criteria) to a representable number before being stored, thus creating a difference between the _real_ number and the _stored_ number. This situation is an example of loss of significance.Does that mean that the _gap_ between representable numbers is _always_ going to be $\epsilon _{mach}$? Of course not! Some numbers will have smaller gaps, and some others will require larger gaps, as studied before. In any interval of the form $[2^n,2^{n+1}]$ for representable $n\in \mathbb{Z}$, the gap is constant. For example, all the numbers between $2^{-1}$ and $2^0$ have a distance of $\epsilon _{mach}/2$ between them. All the numbers between $2^0$ and $2^1$ have a distance of $\epsilon _{mach}$ between them. Those between $2^1$ and $2^2$ have a distance of $2\,\epsilon _{mach}$ between them, and so on.<jupyter_code># What does it mean to store 0.5+delta?
e = 2.**(-1)
f = b/2. # emach/2

result_4 = e + f     # 0.50000000000000011102230246251565404236316680908203125
print("{0:.100}".format(result_4))

result_5 = e + b     # 0.5000000000000002220446049250313080847263336181640625
print("{0:.100}".format(result_5))
g = b/4.

result_5 = e + g     # 0.500000000000000055511151231257827021181583404541015625
print("{0:.100}".format(result_5))<jupyter_output>0.5<jupyter_text>We'll let the students find some representable numbers and some non-representable numbers.<jupyter_code>num_1 = a
num_2 = b
result = a + b
print("{0:.100}".format(result))<jupyter_output>1.0000000000000002220446049250313080847263336181640625<jupyter_text>Loss of significance in function evaluation[Back to TOC](toc) Loss of Significance is present too in the representation of **functions**. A classical example (which you can see in the guide book), is the next function: \begin{equation}f(x)= \frac{1 - \cos x}{\sin^{2}x} \end{equation}Applying trigonometric identities, we can obtain the 'equivalent' function:\begin{equation}f(x)= \frac{1}{1 + \cos x} \end{equation}Both of these functions are apparently equals. Nevertheless, its graphics say to us another thing when $x$ is equal to zero.<jupyter_code>x = np.arange(-10,10,0.1)
y = (1-np.cos(x))/(np.sin(x)**2)
plt.figure()
plt.plot(x,y,'.')
plt.grid(True)
plt.show()
x = np.arange(-10,10,0.1)
y = 1/(1+np.cos(x))
plt.figure()
plt.plot(x,y,'.')
plt.grid(True)
plt.show()
x = np.arange(-1,1,0.01)
y = (1-np.cos(x))/np.sin(x)**2
plt.figure()
plt.plot(x,y,'.',markersize=10)
plt.grid(True)
plt.show()
y = 1/(1+np.cos(x))
plt.figure()
plt.plot(x,y,'.',markersize=10)
plt.grid(True)
plt.show()<jupyter_output><empty_output><jupyter_text>This happens because when $x$ is equal to zero, the first function has an indetermination, but previously, the computer makes a subtraction between numbers that are almost equals. This generates a loss of significance, turning the expression close to this point to be zero. However, modifying this expression to make the second function, eliminates this substraction, fixing the error in its calculation when $x=0$.In conclusion, for us, two representations of a function can be equals, but for the computer they can be different!  Another analysis (example from textbook)[Back to TOC](toc)<jupyter_code>f1 = lambda x: (1.-np.cos(x))/(np.sin(x)**2)
f2 = lambda x: 1./(1+np.cos(x))
f3 = lambda x: (1.-np.cos(x))
f4 = lambda x: (np.sin(x)**2)
x = np.logspace(-19,0,20)[-1:0:-1]
o1 = f1(x)
o2 = f2(x)
o3 = f3(x)
o4 = f4(x)

print("x,            f1(x),        f2(x),        f3(x),                       f4(x)")
for i in np.arange(len(x)):
    print("%1.10f, %1.10f, %1.10f, %1.25f, %1.25f" % (x[i],o1[i],o2[i],o3[i],o4[i]))<jupyter_output>x,            f1(x),        f2(x),        f3(x),                       f4(x)
1.0000000000, 0.6492232052, 0.6492232052, 0.4596976941318602349895173, 0.7080734182735711756961905
0.1000000000, 0.5012520863, 0.5012520863, 0.0049958347219742904599116, 0.0099667110793791851425238
0.0100000000, 0.5000125002, 0.5000125002, 0.0000499995833347366414046, 0.0000999966667111107946004
0.0010000000, 0.5000001250, 0.5000001250, 0.0000004999999583255032576, 0.0000009999996666667110814
0.0001000000, 0.4999999986, 0.5000000012, 0.0000000049999999696126451, 0.0000000099999999666666680
0.0000100000, 0.5000000414, 0.5000000000, 0.0000000000500000041370185, 0.0000000000999999999966667
0.0000010000, 0.5000444503, 0.5000000000, 0.0000000000005000444502912, 0.0000000000009999999999997
0.0000001000, 0.4996003611, 0.5000000000, 0.0000000000000049960036108, 0.0000000000000100000000000
0.0000000100, 0.0000000000, 0.5000000000, 0.0000000000000000000000000, 0.0000000000000001000000000
0.0000000010, 0.0000000000, 0.50[...]
{'fraction_non_alphanumeric': 0.10574974670719352, 'fraction_numerical': 0.22403748733535966, 'mean_word_length': 5.798536375376668, 'pattern_counts': {'":': 0, '<': 67, '<?xml version=': 0, '>': 69, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 83, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '19', 'dir': 'jupyter-structured-clean-dedup', 'id': '96151', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 7637, 'n_tokens_neox': 4749, 'n_words': 1575, 'path': 'SC1v2/02_floating_point_arithmetic.ipynb', 'repo_name': 'FrumpleQuest/Scientific-Computing', 'size': '15792'}