Aritmética finita y fuentes del error¶
Contenidos¶
Sistemas decimal y binario¶
En el sistema decimal el número $105.8125$ significa:
Es decir, es una suma cuyos términos son potencias de $10$ y utiliza $10$ dígitos $0,1,\ldots,9$.
Los ordenadores usan el sistema binario porque así los únicos dígitos que se almacenan son el $0$ y el $1$. En sistema binario los números se representan como sumas de potencias de $2$:
Conversión decimal a binario¶
Para convertir la parte entera, dividimos por $2$ de forma repetida. Los restos de estas divisiones son los dígitos en base $2$, de menos a más significativos.
El número $(105)_{10}$ en sistema binario es $(1101001)_2$.
Nota La función de python bin realiza esta conversión.
bin(105)
Escribir una función de2bi_a que transforme la parte entera
de un número decimal a binario. Aplicarla al número x = 105.8125.
NOTA: pueden ser útiles las siguientes funciones
np.fix(x): redondeaxhacia0.a//bcalcula el cociente de la división entera.a%bcalcula el resto de la división entera.
%run Ejercicio1.py
x = 105.8125
binario1 = de2bi_a(x)
print ''.join(binario1) # si binario1 está almacenado como una lista de cadenas
Para convertir a binario la parte decimal, multiplicamos por $2$, sustraemos la parte entera, que será nuestro dígito binario y repetimos el proceso.
%run Ejercicio2.py
x = 105.8125
binario2 = de2bi_b(x)
print ''.join(binario2) # si binario2 está almacenado como una lista de cadenas
Y el número 105.8125 en binario sería
print ''.join(binario1 + ['.'] + binario2)
Conversión binario a decimal¶
De acuerdo con $(1)$ para convertir $(1101001.1101)_2$ a base $10$
1*(2**6)+1*(2**5)+0*(2**4)+1*(2**3)+0*(2**2)+0*(2**1)+1*(2**0)+1*(2**-1)+1*(2**-2)+0*(2**-3)+1*(2**-4)
Representación de enteros¶
Si tenemos $m$ dígitos o bits de memoria podemos almacenar $2^m$ números diferentes en binario.
Enteros positivos¶
Si sólo tomamos enteros positivos podríamos representar enteros con valores entre
Por ejemplo, para $m=3$ bits podríamos representar los enteros del $0$ al $7$.
Enteros con signo¶
Si queremos representar enteros con signo reservaremos el primer bit para el signo y entonces tendremos espacio para $m-1$ cifras significativas. Es decir $2^{m-1}-1$ números positivos, otros tantos negativos y dos ceros, uno con signo positivo y otro negativo. El rango de números representados sería de $[-2^{m-1}+1, 2^{m-1}-1]$. Si $m=3$ podremos representar números de $-3$ a $3$.
El primer bit es el signo, $0$ para los positivos y $1$ para los negativos.
Ejemplo
Calcular cómo se almacenará $(-80)_{10}$ si usamos representación binaria con una precisión de 8-bits con signo.
Solución: como es un número negativo, la representación es la misma que en el caso anterior pero el primer bit es el signo y como es negativo, es $1$. Es decir, $(-80)_{10}=(11010000)_2$
Representación al complemento de dos¶
Una forma de evitar la doble representación del cero es usar la representación al complemento de dos: los números negativos se construyen tomando los dígitos de los positivos y cambiando los $0$ a $1$ y los $1$ a $0$ y sumándole $1$. De esta manera la suma de un número y su opuesto es siempre $0$.
Ejemplo
Para representar $(-2)_{10}$ empezamos convirtiendo $(2)_{10}$ a binario $(010)_2$
- invertimos sus dígitos $\rightarrow (101)_2$
- le sumamos $001$ (tener en cuenta que en binario $0+0=0 \quad 0+1=1 \quad 1+1=10$) $\rightarrow (110)_2$
Si sumamos $(010)_2+(110)_2=(000)_2$
En este caso, el primer bit de los negativos es $1$ y el de los positivos es $0$. Hemos representado enteros en el rango de $[-2^{m-1},2^{m-1}-1]$.
Ejemplo
Calcular cómo se almacenará $(-80)_{10}$ si usamos representación binaria representación al complemento de dos con una precisión de 8-bits.
Solución: como $(80)_{10}=(01010000)_2$
- en el número anterior cambiamos los $0$ por $1$ y viceversa $\rightarrow 10101111 $
- y sumándole $1$ ya tenemos la representación $(-80)_{10}=(1011000)_2$
Representación sesgada¶
Otra forma es la representación sesgada. Los números negativos se representan con valores positivos consecutivos, empezando por el menor negativo, y los positivos toman los valores siguientes. La representación de los números la obtenemos añadiendo el sesgo $2^{m-1}$ al número $x_r=x+2^{m-1}\in[0,2^m-1]$.
De nuevo vemos que hemos representado enteros en el rango de $[-2^{m-1},2^{m-1}-1]$.
¿De cual de las formas descritas se almacenan los enteros? La mayor parte de las máquinas almacenan los enteros con la representación del complemento al dos y los exponentes de los números en punto flotante, que son enteros, con representación sesgada. ¿Por qué? Hay varios motivos. Uno de ellos es que los números con representación sesgada son más fáciles de comparar. Los códigos comparan frecuentemente números reales, y se empieza comparando el exponente y sólo si son iguales se compara la mantisa.
Representación de enteros en Python¶
Las variables numéricas en Python pueden ser int, float, long y complex. Las variables booleanas son un subtipo de los enteros.
Los enteros sencillos o int tienen al menos 32 bits de precisión, normalmente 32 o 64 bits. Esta depende del ordenador y sistema operativo que estemos utilizando.
El valor máximo que puede tomar un entero sencillo en nuestro ordenador nos lo da
import sys
sys.maxint
En representación sesgada el mayor entero viene dado por $2^{m-1}-1$. Si $m=64$ es el número de bits el máximo entero es
2**(64-1)-1
Y para los valores negativos por $-2^{m-1}$
-sys.maxint-1
-2**(64-1)
Así que, en este caso, el entero sencillo tiene una precisión de 64 bits.
Si necesitamos representar un entero mayor, Python utiliza enteros long. Su precisión solo está limitada por la memoria disponible. Se reconocen porque acaban en L.
sys.maxint+1
-sys.maxint-2
Representación de los números reales¶
En el caso de los números reales se utiliza una representación en punto flotante en base $\beta=2$.
donde:
- $s$ es el signo. Se almacena un $1$ si el número es negativo y $0$ si es positivo.
- $m$ es la mantisa que si está normalizada tiene un valor $1\leq m < \beta$ y el primer dígito ha de ser distinto de $0$. Este dígito en base dos solo puede ser $1$ y por lo tanto, en una mantisa normalizada siempre es $a_0=1$, no es necesario almacenarlo y ganamos un bit. Esto se llama técnica del bit escondido.
- $e$ es el exponente que es un entero con signo y utiliza representación binaria sesgada.
Los números se almacenan en "palabras" de 32 bits (precisión sencilla), 64
bits (precisión doble), o incluso 128 bits (precisión cuádruple). En casi todas las máquinas, la
precisión en Python para el tipo float es la precisión doble.
En el caso de precisión doble los bits se reparten:
- 1 bit para el signo.
- 11 bits para el exponente.
- 52 bits para la mantisa.
Si tenemos $11$ bits para el exponente con signo, significa que tenemos espacio para $2^{11}=2048$ números en formato binario. Si empezamos en el cero, el $ 0 \le E \le 2047 $.
El primer valor $00000000000$ lo reservamos para cero y números
desnormalizados y el último $11111111111$ para Inf y NaN. Así que el rango de valores que puede tomar el exponente queda reducido a $ 1 \le E \le 2046 $.
Y como el sesgo es $2^{m-1}-1=2^{10}-1=1024-1=1023$ estos valores representan a los exponentes $ -1022 \le E - 1023 \le 1023 $.
Por lo tanto el exponente máximo es $e_{max}=1023$, y el mínimo $e_{min}=-1022$. Si dividimos por el valor mínimo:
y no se supera el valor máximo (overflow).
Escribir un programa Ejercicio3.py que calcule la representación de los números los números 105.8125, 120.875, 7.1 y -1.41 en precisión simple, según la norma IEEE 754, utilizando las funciones de los ejercicios 1 y 2. Considerar únicamente el caso en el que la parte entera es mayor que cero. Redondear por truncamiento.
Nota: Recordar que en precisión simple se utilizan 32 bits, de los cuales, uno es para el signo, 8 para el exponente y el resto para la mantisa. El sesgo del exponente es 127.
%run Ejercicio3.py
Las características del tipo float en nuestra máquina nos las da la orden sys.float_info
sys.float_info
Vemos que el número de dígitos de la mantisa, mant_dig es 53 (52 bits pero 53 dígitos debido al bit escondido). Es decir, el float que usa Python por defecto es el de 64 bits.
Pero si observamos esta salida dice
max_exp=1024 y min_exp=-1021
que no concuerda con lo dicho arriba.
Pero si ejecutamos
help(sys.float_info)
entre otras cosas obtenemos
| max_exp | DBL_MAX_EXP -- maximum int e such that radix**(e-1) is representable
| min_exp | DBL_MIN_EXP -- minimum int e such that radix**(e-1) is a normalized float
Como habíamos dicho, el primer $1$ no hace falta guardarlo y nos quedan 52 bits donde almacenamos los unos de $0.11\ldots 11$. Por lo tanto, en decimal este número será:
Hacer un programa Ejercicio4a.py que calcule esta suma. Como norma general, conviene sumar los elementos de la serie de menores a mayores (más adelante veremos por qué). Su valor tiene que coincidir con el valor obtenido con la orden python sys.float_info.max. Llamar sal al valor de salida.
import sys
%run Ejercicio4a.py
print '%.16e' % sal # forma exponencial con 16 decimales
sys.float_info.max
Comprobamos que, efectivamente, son iguales
print sal == sys.float_info.max
(b) El valor real positivo mínimo representable en punto flotante de forma normalizada es
Calcularlo y llamarlo sal. Su valor tiene que coincidir con orden el valor de salida de python sys.float_info.min.
%run Ejercicio4b.py
print '%.16e' % sal
sys.float_info.min
Mayor entero representable exactamente de forma normalizada tal que todos los enteros positivos menores se pueden almacenar de forma exacta en forma normalizada
El mayor entero consecutivo, utilizando el formato en punto flotante, del que podemos almacenar todos los dígitos significativos y que, por tanto, podemos almacenar de forma exacta en binario podría ser:
Es decir, tenemos $53$ cifras de precisión en binario. El entero siguiente:
también podemos almacenarlo de forma exacta. Pero el siguiente a este ya no porque necesitaríamos un bit más a la derecha en la mantisa. Este entero es
print(2**53)
Todos los enteros menores se pueden almacenar de forma exacta. Si contamos, vemos que tiene 16 dígitos. Por lo tanto, podemos almacenar todos los enteros con 15 dígitos decimales y casi todos los enteros con 16 dígitos de forma exacta.
Pero, podría decirse, el número
es más grande y también es entero. Sin embargo sus $10$ últimas cifras serían forzosamente $0$ porque no las podemos almacenar y desde el punto de vista de la precisión no cuentan porque no les podemos dar el valor que queramos.
Valores desnormalizados
¿Qué sucede con los valores de los exponentes fuera del rango $[-1022,1023]$?
Si el exponente es menor que $-1022$ , el número es desnormalizado o cero y estamos usando para los bits correspondientes al exponente el valor especial $00000000000$. Entonces ya no se asume que el bit escondido es $1$ sino $0$.
x = 0.5*(2**-1023)
print x # desnormalizado
print 1/x # su inverso
print 2**-1080 # Underflow
El menor valor desnormalizado es
es decir, $2^{-1022-52}$
num_min = 2.**(-1022-52)
print num_min
Si el valor es menor, devuelve es Underflow y devuelve cero.
2.**(-1022-53) # Underflow
Si el exponente es mayor que $1023$ devuelve OverflowError.
2.**1024 # Overflow
Precisión de la máquina¶
Vamos a calcular el número más pequeño que se puede representar después de $1$ utilizando la representación binaria en punto flotante con doble precisión:
El número real $1$ en representación en punto flotante de doble precisión normalizada (con $52$ ceros en la mantisa) es
y no podemos representar, de forma exacta en punto flotante, un número real entre los dos. La diferencia entre $1$ y $1+\epsilon$ es
que pasado a decimal es
1*2.**(-52)
o lo que es lo mismo
sys.float_info.epsilon
Este valor, el menor número $\epsilon$ tal que $1+\epsilon>1$ se llama precisión de la máquina y es una cota del error relativo que cometemos al almacenar cuanlquier número real en punto flotante, en este caso, en doble precisión.
Como, en doble precisión, $\epsilon\approx 2.22\cdot 10^{-16}$ se puede decir que corresponde, aproximadamente, a $16$ cifras decimales.
Si queremos el número más pequeño comparable con $x$ hay que sumarle un valor eps
x = 1.
eps = np.spacing(x)
print(eps)
x = 10.
eps = np.spacing(x)
print(eps)
x = 100.
eps = np.spacing(x)
print(eps)
Y vemos que la diferencia entre dos números consecutivos representables de forma exacta en punto flotante, aumenta conforme nos alejamos del cero. Por ello, la densidad de números representados en punto flotante no es uniforme. Es mayor cerca del cero y menor a medida que nos alejamos.
Resultados especiales¶
Algunas operaciones dan resultados especiales
x = 1.e200
y = x*x
print(y)
y/y
nan significa "not a number". Su aparición significa, en general,
que algo ha ido mal con el programa y/o que se ha realizado una operación
no válida.
Overflow. Ocurre cuando el resultado de la operación realizada es finita pero el resultado supera el rango de reales que pueden ser representados.
x = 1.e300
x**2
Aproximación y error de redondeo¶
Redondeo
A menudo, para un número real $x$, no exite una representación exacta en punto flotante y cae entre dos números racionales consecutivos $x^- \lt x \lt x^+$. Como representación de $x$ elegiremos uno de los dos dependiendo del método de redondeo usado. IEEE propone 5 sistemas de redondeo:
- Hacia arriba.
- Hacia abajo.
- Hacia cero ("truncamiento").
- Hacia el más alejado de cero más cercano.
- Hacia el par más cercano.
Este último método de redondeo es el más utilizado.
Pérdida de dígitos¶
Supongamos que estamos usando un sistema de representación decimal que almacena 3 cifras significativas y queremos sumar los números $a=123$ y $b=1.25$. El resultado es $s=124.25$ pero como solo podemos almacenar tres cifras redondeamos, por ejemplo, a $s_r=124$. Hemos perdido cifras y estamos cometiendo un error. En general, cualquier operación aritmética conlleva errores de redondeo.
Ejemplo
Si realizamos la operación $3\times 0.1 \times 100$:
3*0.1*100
El resultado no es $30$ debido al error de redondeo. Estos se pueden propagar y afectar hasta invalidar los resultados. En el caso anterior el error se debe a que, como vimos al principio, el número $x=0.1$ en base dos es periódico y por ello su representación en el ordenador no es exacta.
Ejemplo
Si calculamos $$\sum_{k=1}^{10000} 0.1$$
suma = 0.
for i in range(10000):
suma += 0.1
print '%.16f'% suma # En forma decimal con 16 decimales
y el error absoluto es
import numpy as np
abs(suma-1000.)
Ejemplo
Si sumamos o restamos dos números muy diferentes en magnitud se pierde exactitud debido a los errores de redondeo. Por ejemplo
a = 1.e+9
epsilon = 1.e-8
suma = a
for i in range(10000):
suma += epsilon
print '%.16e'% suma # En forma exponencial con 16 decimales
Sin embargo el resultado debería haber sido
suma = a + 10000*epsilon
print '%.16e'% suma
Ejemplo
Vamos a sumar los $N$ primeros términos de la serie armónica $$\sum_{n=1}^N \frac{1}{n}$$
El resultado teórico para 1000 términos sería
import numpy as np
from scipy import special
N = 10000
suma = special.polygamma(0,N+1)-special.polygamma(0,1)
print '%.20f'% suma
Empezando la suma por el término $n=1$
suma1 = 0.;
for n in range(1,N+1):
suma1 += 1./n
print '%.20f'% suma1
error1 = abs(suma1 - suma)
print error1
Empezando por el término final
suma2 = 0.;
for n in range(N,0,-1):
suma2 += 1./n
print '%.20f'% suma2
La diferencia entre ambos resultados es
error2 = abs(suma2 - suma)
print error2
Es mayor el error en el caso primero porque los términos de la sucesión son cada vez más pequeños mientras que la suma es cada vez más grande. Por lo tanto, conforme avanzamos en la suma se producirá pérdida de dígitos por sumar números de magnitud muy distinta.
Error de cancelación¶
Aparece al restar números de valores muy parecidos. Por ejemplo
import numpy as np
x = 1000000.
ep = 0.01
a = np.sqrt(x**2 + ep**2) - x
print a
En este caso el error de calcelación se puede evitar utilizando una expresión matemática equivalente
b = ep**2 / (np.sqrt(x**2 + ep**2) + x)
print b
O bien
b = ep**2 / (2*x)
print b
El error relativo cometido al calcularlo con la primera expresión es muy grande
Er = abs((a - b)/b)
print Er
Además, este ejemplo muestra que expresiones matemáticamente equivalentes no son, necesariamente, equivalentes desde el punto de vista computacional. Y viceversa, que expresiones matemáticas que no son iguales, pueden tener el mismo resultado en el ordenador y ser, por tanto, computacionalmente equivalentes.
(a) Al resolver la ecuación de segundo grado
utilizando la fórmula
una de las soluciones es:
a = 1.; b = 10.**8; c = 1.
x1 = (-b + np.sqrt(b**2 - 4*a*c)) / (2*a)
print x1
Pero al sustituir en el polinomio para calcular el residuo
residuo1 = a*x1**2 + b*x1 + c
print residuo1
es relativamente grande. Encontrar una expresión matemática equivalente que dé una solución con menor residuo y calcular de nuevo la solución y su residual. Escribir para ello un programa llamado Ejercicio5a.py
%run Ejercicio5a.py
print x2
print residuo2
(b) Al resolver la ecuación de segundo grado
utilizando la fórmula
una de las soluciones es:
a = 1.; b = -10.**8; c = 1.
x3 = (- b - np.sqrt(b**2 - 4*a*c)) / (2*a)
print x3
Pero al sustituir en el polinomio
residuo3 = a*x3**2 + b*x3 + c
print residuo3
es relativamente grande. Encontrar una expresión matemática equivalente que dé una solución con menor residual y calcular de nuevo la solución y su residual. Escribir para ello un programa llamado Ejercicio5b.py.
%run Ejercicio5b.py
print x4
print residuo4
Ejemplo
Tenemos que
Aproximemos el valor de $e$, usando $x=1$, con un número finito de términos de la suma, por ejemplo, $50$
import math
suma = 0.
x = 1.
for n in range(100):
suma += x**n / math.factorial(n)
if n in np.array([5,10,15,20,40,60,80,100])-1:
error = abs(np.exp(1.) - suma)
print u"Número de términos %i" % (n+1)
print "Error %.16e" % error
print ""
Vemos que primero el error disminuye y luego se estanca. ¿Por qué? Si estuviéramos usando aritmética infinita el error disminuiría con las iteraciones. Pero como estamos usando aritmética finita, el error, en este caso, está limitado por el espaciado entre números representables en la zona de la solución $e$ que es
Ea = np.spacing(np.exp(1.))
print Ea
Calculemos el menor número de términos necesarios para obtener el número $e$ con el menor error posible.
import math
x = 1.
Ea = np.spacing(np.exp(x))
n = 0
suma = 1.
itermax = 100
error = abs(np.exp(x) - suma)
while (error > Ea and n < itermax):
n += 1
suma += x**n / math.factorial(n)
error = abs(np.exp(x) - suma)
print u"Número de iteraciones %i" % n
print "Error %.16e" % error
print "Espaciado alrededor de e %.16e" % Ea
Y aunque sabemos que con $17$ iteraciones todavía cometemos error, este ya no es apreciable por la máquina y aunque realicemos más iteraciones este error computado no disminuye.
calcular el número de términos necesarios para calcular $\log 1.5$ con el menor error posible.
%run Ejercicio6.py
print "Número de iteraciones %i" % n
print "Error %.16e" % error
print "Espaciado alrededor de ln(1.5) %.16e" % Ea