4 Numbers and Floating Point Numbers

4.1 Representation of Numbers

We usually represent numbers in the decimal system. For example, the number \(273\) we interpret as \[\begin{aligned} 273 &=& 2\cdot 100 + 7\cdot 10 + 3\cdot 1\\ & = & 2\cdot 10^2 + 7\cdot 10^1 + 3\cdot 10^0 \, . \end{aligned}\] Therefore, we talk about a decimal representation, as numbers are expressed in multiples of powers of \(10.\) However, looking at the above number, we see that we could also write it as \[\begin{aligned} 273 & = & 256 + 16 + 1 \\ & = & 1 \cdot 2^8 + 1 \cdot 2^4 + 1 \cdot 2^0 \\ & = & 1 \cdot 2^8 + 0 \cdot 2^7 + 0 \cdot 2^6 + 0 \cdot 2^5 + 1 \cdot 2^4 + 0 \cdot 2^3 1 + 0 \cdot 2^2 + 0 \cdot 2^1 +1 \cdot 2^0 \end{aligned}\] i.e. using powers of two. Thus, we see that the decimal number \((273)_{10},\) where the subscript \({}_{10}\) indicates that this number is to be interpreted with respect to the basis \(10,\) can equally be represented as the binary number \((100010001)_{2}.\)

We now investigate how natural numbers can be represented in a \(b\)-adic system.

The proof is carried out via induction. For the proof we will need the floor and ceiling functions (Gaussian brackets) which are defined for \(x\in\mathbb{R}\) as \[\begin{aligned} \lfloor x\rfloor &=& \max \{k\in\mathbb{Z}\colon k\leq x\}\\ \lceil x\rceil &=& \min \{k\in\mathbb{Z}\colon k\geq x\} \end{aligned}\] as well as the modulo operation mod defined for \(x\in\mathbb{Z}\) and \(y\in \mathbb{N}\setminus\{0\}\) as \[x\mod y = x-y\cdot \left\lfloor \frac{x}{y}\right\rfloor\, .\] From the above definitions, it follows that for \(z\in\mathbb{N}\) we have \[z=\lfloor z/b\rfloor\cdot b + (z \mod b)\, .\]

4.2 Floating Point Representation and Machine Accuracy

Computers are predominantly based on numbers systems with basis \(b=2.\) We therefore consider the representation of integers in a computer using bits \(d\in \{0,1\}.\)

4.2.1 Integers

For representing an integer \(z\in\mathbb{Z}\) we use \(N\) bits, i.e. \[z \hat{=} d_{N-1}d_{N-2}\cdots d_0\quad\text{ with }d_i\in\{0,1\}.\] That means \[\tag{4.2} z = d_0+d_1 2^1 + d_2 2^2 + \ldots + d_{N-1} 2^{N-1}.\] Natural numbers are represented as unsigned integers; integers additionally have one bit for the sign \(s\in\{0,1\}.\) That means \[z \hat{=} s|d_{N-2}d_{N-3}\cdots d_0 \hat{=} (-1)^s\sum\limits_{i=0}^{N-2}d_i2^i.\]

In the case of the above integer representation, zero is not represented uniquely; moreover, subtraction cannot be treated the same way as addition. Thus, for negative numbers, we use the following representation as two’s complement \[1|d_{N-2}\ldots d_0\hat{=}-\left(1+\sum\limits_{i=0}^{N-2}(1-d_i)2^i\right) \, .\] For subtraction, each \(z<0\) is represented as above. Then all bits are changed and finally \(1\) is added.

For given \(N,\) we can represent integers within the following range \[z_{\text{min}}= -2^{N-1}=-\left(1+\sum\limits_{i=0}^{N-2}1\cdot2^i\right) \leq z \leq \sum\limits_{i=0}^{N-2}1\cdot2^i=2^{N-1}-1=z_{\text{max}}\] Usually \(N = 32\) is called integer and \(N=64\) long integer. Numbers \(z\) with \(z<z_{\text{min}}\) or \(z>z_{\text{max}}\) lead to underflow or overflow, respectively.

4.2.2

The representation of integers within the admissible range is exact. For real numbers, this is not possible.

We use an approximate representation of real numbers. Thereby, we aim for

• a small rounding error

• a huge range of numbers, which can be represented

• stability of the elementary operations \(\{+,-,\cdot,/\}\)

We now have two possibilities: We can approximate a real number \(x\in\mathbb{R}\) with a floating point number \(\text{rd}(x)\) such that either the absolute error is bounded by some \(\varepsilon > 0,\) i.e. \[|x-\text{rd}(x)|\leq\varepsilon,\] or such that the relative error is bounded: \[|x-\text{rd}(x)|\leq\varepsilon|x|\, .\] We choose the relative error, since it is in general the appropriate error measure and since this choice allows for a much larger range of representable numbers.

The proof may be carried out as an exercise.

By means of the mapping \(\text{rd}: [x_{\min},x_{\max}]\longrightarrow {\mathcal G},\) each \(x\in\mathbb{R}\) gets assigned one and only one floating point number. We moreover have \(\text{rd}\circ\text{rd}=\text{rd},\) i.e. \(\text{rd}\) is a projection. We call the number \(\text{rd}(x)\) floating point representation of \(x.\) What does \(1\)% error mean? It means, there exists \(\tilde{x}\) with \[\frac{|x-\tilde{x}|}{|x|}\leq \frac{1}{100}\, .\]

The mapping \(\text{rd}\) is not injective. Several different numbers will be rounded to the same floating point number. Computing in \({\mathcal G}\) is computing with intervals: each floating point number \(\tilde{x}\) represents an interval around \(\tilde{x}\) (the size of the interval can vary with \(x\)). We set \[x\stackrel{{\mathcal G}}{=}y:\iff x\tilde{=}y:\iff \text{rd}(x)=\text{rd}(y).\] As a consequence, within source code statements such as x==y should never be used. Instead, \(|x-y|\leq r\) should be taken, where \(r\) is suitably choosen (e.g. \(r=\mathrm{eps}\cdot |x|\)). Let us finally note that in the denormalized floating point representation the relative error is not bounded by \(\mathrm{eps}\) any more.

The IEEE-standard (IEEE : Institute of Electrical and Electronics Engineers) specifies the data types float and double as follows:

We call the \(32\)-bit number float (single precision), and the \(64\)-bit number double (double precision).

For double, this leads —taking into account the hidden bit— to \(\mathrm{eps}=\frac{1}{2} 2^{-52}=2^{-53}\approx 1.1 \cdot 10^{-16}.\) Often, for the machine accuracy \(\hat{\mathrm{eps}}_{\text{double}}=2.2 \cdot 10^{-16} =2\cdot \mathrm{eps}\) can be found. This does not correspond to the definition given here.

4.2.3 Floating Point Arithmetic and Cancellation of Leading Digits

Ideally, when doing computations in floating point arithmetic, the relative error of the result of one of the elementary operations \(\{+,-,\cdot,/\}\) should be bounded by the machine accuracy. Let us therefore start with the following

We now investigate the influence of relative errors on the elementary operations. We start with the addition. Let \(x,y\) be floating point numbers. Then, \(x+y\) can have a mantissa with length \(>t.\) For example, let \(\beta=10, L=-9, U=9\) and \(t=2.\) Then we have \[0.10 \cdot 10^0 + 0.10 \cdot 10^2 = 0.101 \cdot 10^2\,\hat{=}\, 0.10\cdot 10^2\,.\] Due to \(t=2,\) on the right hand side information is destroyed. In case of \[0.51 \cdot 10^9 + 0.5 \cdot 10^9 = 0.101\cdot 10^{10}\hat{=} 0.1 \cdot 10^{10}\notin{\mathcal G}\] we find an overflow in the exponent. An underflow in the exponent results from \[0.51 \cdot 10^{-9} - 0.5 \cdot 10^{-9} = 0.10 \cdot 10^{-10} \notin{\mathcal G}.\]

4.2.4 Cancellation of leading digits

The subtraction of approximately equal numbers in a floating point system can lead to unwanted consequences in terms of accuracy. We illustrate this with the following two examples:

• Let \[\begin{aligned} a&=&0.123467\!\ast\\ b&=&0.123456\ast \end{aligned}\] be given. Both numbers might be disturbed at the 7th digit. Then, we have \[a-b=0.000011\ast=0.11\!\ast\cdot10^{-4}.\] An error in the data in the 7th digit leads to an error in the result in the 3rd digit. This means \(\kappa_{\text{rel}}\approx10^{4}.\)

• Consider the polynomial \[f(x,y)=(x-2y)^3=x^3-6x^2y+12xy^2-8y^3\] for \(x=10 000 000,\ y=4 999 999.\) We have \[\begin{aligned} a&:=&x^3+12xy^2\approx3.9999000001202585\cdot10^{21}\\ b&:=&8y^3+6x^2y\approx3.9999000001973427\cdot10^{21} \end{aligned}\] and, in floating point arithmetic (double), \[a-b = 524.388.\] For \(\mathrm{eps}\approx 1\cdot 10^{-16}\) the error is increased as follows: \[\kappa_{\text{rel}}=\frac{a+b}{|a-b|}\approx \frac{8\cdot10^{21}}2=4\cdot10^{21} \, .\]

As a consequence we see that the subtraction of approximately equal numbers should be avoided. If it is unavoidable, it should be carried out at the beginning of an algorithm.

4.2.5 Special Numbers and their Representation

Certain operations will not lead to a floating point number but to a so called “exception”. We might encounter overflow or underflow, or division by \(0\) might occur. In these cases, special numbers will be returned to the programmer.

4.2.6 Summation Techniques

We would like to illustrate that the floating point addition \(\hat{+}\) in general is not associative.

We consider infinite series \[1+\sum\limits_{k=2}^\infty \frac {1}{k} = 1+ \frac{1}{2} +\frac{1}{3} +\cdots \, ,\] which is called harmonic series. The harmonic series is said to converge iff the sequence \((S_K)_{K\in \mathbb{N}}\) of the so called partial sums \[S_K=\sum\limits_{k=2}^{K} \frac {1}{k}\] converges. Since the harmonic series is divergent (as you will see in Calculus 1), also the sequence \((S_K)\) of the partial sums diverges. However, in floating point arithmetic we expect the summation to stall at a certain point, as the sequence \(x_k=\frac{1}{k}\) converges to \(0.\) In order to investigate this, we compute the \(2^m\)–th partial sum, \(m>0\) of the above series, i.e. \[S_{2^m}=1+\sum\limits_{k=2}^{2^m} \frac {1}{k}\, .\] For doing so, we use three different summation techniques.

The first technique is the intuitive forward summation:

fwd = 1;
for k=2:n
 fwd = fwd + 1/k;
end

The second technique is backward summation:

n = 2^fix(m);   
bwd = 0;
for k=n:-1:2
 bwd = bwd + 1/k;
end
bwd = bwd + 1;

The third technique is the so called hierarchical summation.

n = 2^fix(m);   
a1 = 0;
a2 = 0;
I=[3:2:n];
n = n/2;
sum = zeros(1,n);
sum(1) = 1 + 1/2;
for k=I
  a1 = 1/k;
  a2 = 1/(k+1);
  sum(1+(k-1)/2) = a1+a2;
end

while(n > 2)
  for i=1:2:n
   sum(i) = sum(k)+sum(k+1);
  end
  sum = sum([1:2:n]);
  n = n/2;
end
hier = sum(1)+sum(2); 

We illustrate the hierarchic summation for \(K=8.\) We have \[\begin{aligned} S_8 &= x_1+ \cdots + x_8 \\ & = (\underbrace{(x_1+x_2)}_{=s_{1,2}}+\underbrace{(x_3+x_4)}_{=s_{3,4}} + \underbrace{(x_5+x_6)}_{=s_{5,6}}+\underbrace{(x_7+x_8)}_{=s_{7,8}}) \\ & = \underbrace{(s_{1,2}+s_{3,4})}_{=s_{1,2,3,4}} +\underbrace{(s_{5,6}+s_{7,8})}_{=s_{5,6,7,8}}\\ & = s_{1,2,3,4} + s_{5,6,7,8} \, . \end{aligned}\] Figure 4.1 shows a comparison of the three different approaches. As can be seen, the summation stalls already very relatively small values of the partial sums, as for the individual contributions \(1/k\) it also holds \(1/k\to 0\) for \(k\to\infty.\) The best results are achieved for the hierarchical summation.

Algorithmics

Home

Chapters

Contents

---

✕