Aritmética de punto flotante

Ejecutar este documento en forma dinámica: Binder

Vamos a hacer una cuenta simple:

In [1]:
0.1 + 0.1 + 0.1 - 0.3
Out[1]:
5.551115123125783e-17

En la matemática que aprendimos en la escuela, $a + a + a - 3 a = 0$, pero cuando ejecutamos la instrucción de arriba el valor obtenido no es cero (aunque es muy pequeño).

Esto no es un error de Python, sino que ocurre debido a la representación interna que hacen las computadoras de los números que hace que no haya forma de representar exactamente valores como 0.1. Esta representación interna es de punto flotante binario. Vamos a ver primero estos conceptos, y luego los problemas que ocurren debido a esta representación finita de los números reales.

Punto flotante decimal

Dado que la memoria de las computadoras es limitada, no es posible almacenar números con precisión infinita. Por lo tanto, de algún modo es necesario realizar un corte en la secuencia de dígitos que representa un números real. Esto es algo que hacemos rutinariamente en todos nuestros cálculos. Al número $1/2$ lo podemos representar exactamente como 0.5, pero al número $1/3$ podemos representarlo, con precisión creciente, como:

0.3
0.33
0.333
0.3333
0.33333
...

La representación de punto flotante es una forma de notación científica usada en las computadoras con la cual se pueden representar números reales extremadamente grandes y pequeños de una manera muy eficiente y compacta, y con la que se pueden realizar operaciones aritméticas. El estándar actual para la representación en coma flotante es el IEEE 754.

Esta representación consiste en descomponer un número en tres partes, de la siguiente forma:

$$r = (-1)^s c \; b^e$$

donde $r$ es el número que queremos representar, $s$ es el signo, $c$ es la mantisa, $b$ es la base de representación (10: decimal, 2: binario, 16: hexadecimal) y $e$ es el exponente. Por ejemplo:

\begin{align*} s &= 0 \\ c &= 31415926 \\ b &= 10 \\ e &= -7 \\ r &= (-1)^0 31415926 \times 10^{(-7)} = 3.1415926 \end{align*}

Este formato satisface los siguientes requisitos:

  • Puede representar números de órdenes de magnitud enormemente dispares (limitado por la longitud del exponente).
  • Proporciona la misma precisión relativa para todos los órdenes (limitado por la longitud de la mantisa).
  • Permite cálculos entre magnitudes: multiplicar un número muy grande y uno muy pequeño conserva la precisión de ambos en el resultado.

Los números de punto flotante decimales se expresan normalmente en notación científica con un punto explícito siempre entre el primer y el segundo dígitos. El exponente o bien se escribe explícitamente incluyendo la base, o se usa una e (o E) para separarlo de la mantisa.

Punto flotante binario

Las computadoras utilizan el mismo formato pero binario, es decir, $b = 2$. Los formatos más utilizados son de 32 o 64 bits en total (4 u 8 bytes):

Precisión Bits Bits significativos Bits del exponente Número más pequeño Número más grande
Simple 32 23 + 1 signo 8 ~1.2e-38 ~3.4e38
Doble 64 52 + 1 signo 11 ~5.0e-324 ~1.8e308

Algunas características:

  • En la secuencia de bits, el primero es el del signo, luego el exponente y finalmente la mantisa.
  • El exponente no tiene signo, en su lugar se le resta un desplazamiento (bias) (127 para simple y 1023 para doble precisión). Esto, junto con la secuencia de bits, permite que los números de punto flotante se puedan comparar y ordenar correctamente incluso cuando se interpretan como enteros: $r = (-1)^s \, c \, b^{(e - E)}$, con $E$ denotando el bias.
  • Se asume que el bit más significativo de la mantisa es 1 y se omite, excepto para casos especiales.
  • Hay valores diferentes para cero positivo y cero negativo. Estos difieren en el bit del signo, mientras que todos los demás son 0. Deben ser considerados iguales aunque sus secuencias de bits sean diferentes.
  • Hay valores especiales no numéricos (NaN, not a number en inglés) en los que los bits del exponente son todos 1 y los de la mantisa no son todos 0. Estos valores representan el resultado de algunas operaciones indefinidas (como multiplicar 0 por infinito, operaciones que involucren NaN, o casos específicos). Incluso valores NaN con idéntica secuencia de bits no deben ser considerados iguales.

Repaso (?): conversión de decimal a binario

Para convertir un número decimal n, por ejemplo 2.71728 a binario con una cierta cantidad k de dígitos (por ejemplo 4), se deben seguir los pasos a continuación:

1- Convertir la parte entera de n en el binario equivalente:

1. Dividir el número decimal por 2 y guadar el resto en una lista.
1. Dividir el cociente por 2.
1. Repetir el paso B hasta que el cociente sea 0.
1. El número binario equivalente resulta de la lista del paso A en orden inverso.

Divisor Cociente Resto
2 1 0
1 0 1
$$2_{10} = 10_{2}$$

2- Convertir la parte fraccionaria de n en el binario equivalente:

1. Multiplicar la parte fraccionaria decimal por 2.
1. La parte entera del resultado de A es el primer dígito de la parte fraccionaria binaria.
1. Repetir el paso 1 usando solo la parte fraccional del número y luego el paso B, hasta obtener los `k` dígitos.

Fracción decimal x 2 Parte entera del resultado
0.71728 1.43456 1
0.43456 0.86912 0
0.86912 1.73824 1
0.73824 1.47648 1
$$0.71728_{10} = 1011_{2} $$

3- Concatenar la parte entera binaria con la fraccionaria binaria:

$$2.71728_{10} = 10.1011_{2}$$

Podemos hacer un script que realice el cálculo por nosotros:

In [ ]:
def dec2bin(n, k):
    binario = ''
    entero = int(n)
    fraccional = n - entero
    # Convertimos la parte entera a su formato binario
    while(entero):
        resto = entero % 2
        binario += str(resto)
        entero //= 2
    binario = binario[::-1]
    binario += '.'
    # Convertimos la parte fraccional a su formato binario
    while(k):
        fraccional *= 2
        fraccional_bit = int(fraccional)
        if fraccional_bit == 1:
            fraccional -= fraccional_bit
            binario += '1'
        else:
            binario += '0'
        k -= 1
    return(binario)
In [ ]:
dec2bin(2.71728,4)

Para volver de binario a decimal, la operación es simple. Por ejemplo, si queremos llevar el número $110.001_b$ a decimal:

\begin{align*} 110.001_b &= 1 \, 2^2 + 1 \, 2^1 + 0 \, 2^0 + 0 \, 2^{-1} + 0 \, 2^{-2} + 1 \, 2^{-3} \\ &= 1 \, 4 + 1 \, 2 + 0 \, 1 + \frac{0}{2} + \frac{0}{4} + \frac{1}{8} \\ &= 6 + \frac{1}{8} = 6.125 \end{align*}

Podemos también hacer un script que transforme de binario a decimal, dando como entrada una string que representa el número binario:

In [ ]:
def bin2dec(s):
    if not '.' in s:
        be, bd = s, ''
    else:
        be, bd = s.split('.')        
    d = 0
    print(be,bd)
    for e in range(len(be)):
        d += int(be[-(e + 1)]) * 2**(e)
    for e in range(len(bd)):
        d += int(bd[e]) * 2**(-(e + 1))
    return(d)
In [ ]:
bin2dec('110.001')

¿Qué sucede cuando queremos volver a decimal el número $10.1011_b$? ¿Por qué ocurre esto?

Número decimal a representación binaria IEEE

Para convertir un número decimal a punto flotante IEEE hay que seguir los siguientes pasos:

  1. Escribir el número en binario
  2. Normalizar el número
  3. Encontrar la mantisa y el exponente
  4. Escribir el número como punto flotante binario IEEE

Ejemplo: convertir $2.25_{10}$ a punto flotante binario IEEE de precisión simple (32 bits).

In [ ]:
dec2bin(2.25,8)

Para normalizar $10.01000000_2$ hay que multiplicar por $2^1$. La mantisa entonces es 001 0000 0000 0000 0000 0000 (7 + 8 + 8 = 23), recordando que se asume que hay un 1 a la izquierda del punto decimal ya que el valor está normalizado.

Para hallar el exponente debemos recordar que en representación de 32 bits, el bias es 127, por lo que el exponente en la representación IEEE es $e - 127 = 1 \Rightarrow e = 128$, que convertido a binario es:

In [ ]:
dec2bin(128,0)

Finalmente, la representación en punto flotante binario IEEE resulta:

s e         c
0 10000000 (1)001 0000 0000 0000 0000

o directamente:

0100000000010000000000000000

Representación binaria IEEE a decimal

Para realizar la conversión inversa hay que recorrer el camino inverso al de la sección anterior, esto es, reconocer la mantisa y el exponente, convertir a decimal las expresiones binarias y realizar la aritmética necesaria. Por ejemplo, dado el binario siguiente (de 32 bits):

0 10000001 10000000000000000000000

obtener su representación decimal.

Dado que el primer bit es 0, determinamos que el número es positivo. Los siguientes ocho bits representan el exponente en formato binario, por lo que debemos convertirlo a entero:

In [ ]:
bin2dec('10000001')

Dado que el bias en 32 bits es 127, el correspondiente exponente es $129 - 127 = 2$. La mantisa es $0.1 = 2^{-1} = 0.5$, con lo que el decimal equivalente es:

$$ +1.5 \; 2^{2} = +6.0 $$

Errores

En el transcurso del tratamiento matemático de datos debemos tratar con el problema de los errores. Existen cuatro formas en que se introducen errores en los cálculos.

  1. Error inherente. Estos errores provienen desde el principio en los datos originales y están fuera del alcance del control de cálculo. Pueden deberse, por ejemplo, de las incertezas inherentes al proceso de medición.
  2. Error de truncamiento. Se producen como resultado de reemplazar un proceso infinito por uno finito, tal como sucede cuando consideramos solo los primeros términos de una expansión en serie de Taylor.
  3. Error de redondeo. Estos errores provienen de la representación con precisión finita de los números en una computadora (de este error nos ocuparemos aquí).
  4. Error por equivocación. Este error es causado por realizar una operación aritmética incorrectamente.

Los primeros tres errores son inevitables, por lo que el "problema de los errores" no consiste en evitarlos, sino en controlar su magnitud con el propósito de obtener un resultado final lo más preciso posible. Este proceso se denomina control de error y consiste en estimar la propagación de los errores a través de un cómputo. Por ejemplo, queremos asegurarnos que los resultados de realizar una operación aritmética entre dos números (que intrínsecamente ya tienen un error), estén dentro de límites tolerables. Además, este error en cálculos subsecuentes también debe mantenerse bajo control.

El error por equivocación es evitable, por lo que no debe ocurrir y por lo tanto no lo consideraremos en este curso.

Error de redondeo

Un número en representación entera es exacto, y por lo tanto la aritmética entre números enteros también es exacta con dos salvedades: (1) el resultado no está por afuera del rango de representación de enteros (usualmente con signo), y (2) la división se interpreta generando también un entero, y descartando todo resto (también entero).

Diferentes patrones de bits pueden representar el mismo número en punto flotante. Si usamos una representación binaria, por ejemplo, una mantisa con los bits más significativos nulos puede ser desplazada a izquierda, es decir, multiplicado por una potencia de 2, siempre que el exponente sea disminuido por la misma cantidad. Los patrones de bits que están tan desplazados a izquierda como se pueda, se dicen que están normalizados. Casi todas las computadoras siempre producen resultados normalizados, dado que no se quiere desperdiciar bits para las mantisas y por lo tanto ganar en precisión en la representación. Dado que el bit de mayor orden de una matisa normalizada (en representación binaria) es siempre 1, algunas computadoras ni siquiera almacenan este bit, ofreciendo entonces un bit extra de significancia.

La aritmética entre números en representación de punto flotante no es exacta, aún si los operandos lo son. Por ejemplo, la suma de dos números en punto flotante se realiza primero desplazando a derecha (dividiendo por dos) la mantisa del número menor (en magnitud), y incrementando simultáneamente el exponente, hasta que los dos operandos tengan el mismo exponente. Los bits menos significativos se pierden debido a este desplazamiento. Si los dos operandos difieren mucho en sus magnitudes, el operando menor es reemplazado efectivamente por cero, dado que es completamente desplazado a derecha.

Por ejemplo:

In [ ]:
a = 80000
b = 1.43E-12
a + b

Veamos cómo sucede esto analizando las representaciones de números, y la suma, en punto flotante decimal IEEE (no usamos los del ejemplo anterior ya que la suma fue calculada con representación de 64 bits, y por comodidad vamos a usar 32).

Haremos $4 + 10^{-7}$. Primero encontraremos su representación binaria IEEE:

In [ ]:
print(dec2bin(4, 0))
print(dec2bin(1.e-7,32))

Estos números se representan (con el bit más significativo de normalización explícito):

   4 = 0 10000001 1000 0000 0000 0000 0000 0000
1e-7 = 0 10010111 1101 0110 1000 0000 0000 0000
     = 0 10000001 0000 0000 0000 0000 0000 0000 (con el mismo exponente que 4)

El desplazamiento hacia un exponente común es necesario para que dos números puedan sumarse. Vemos entonces que todas las cifras significativas se pierden en la representación de $10^{-7}$, por lo que al sumar $4 + 10^{-7}$ obtenemos $4$, aún cuando podemos representar $10^{-7}$ en forma precisa.

El número más chico (en magnitud) en representación de punto flotante que, cuando se suma al número $1.0$ en punto flotante, produce un punto flotante diferente de $1.0$ se denomina precisión de la máquina $\epsilon_m$. En precisión simple tiene un valor cercano a $1.19 \times 10^{-7}$, mientras que para saber cuál es el valor de la máquina que ejecuta este código (usalmente en doble precisión o 64 bits), podemos usar:

In [ ]:
import sys
print(sys.float_info)

Vemos aquí que $\epsilon_m = 2.220446... \times 10^{-16}$. En cierto modo, $\epsilon_m$ es la precisión relativa con la cual se pueden representar números en punto flotante, correspondiente a un cambio de $1$ en el bit menos significativo de la mantisa. Se puede interpretar esto como que cualquier operación aritmética entre números de punto flotante introducen un error fraccional adicional de al meno $\epsilon_m$. Éste es el error de redondeo.

Es importante notar que $\epsilon_m$ no es el número de punto flotante más chico que puede representar la computadora. Ése número depende de cuántos bits hay en el exponente, mientras que $\epsilon_m$ depende de cuántos bits hay en la mantisa.

Al realizar cálculos repetidos, el error de redondeo se acumula. Si en un determinado cómputo realizamos $N$ de tales operaciones aritméticas, un cálculo optimista nos daría un error de redondeo total del orden de $\sqrt{N} \epsilon_m$, si los errores de redondeo se distribuyen aleatoriamente hacia arriba y hacia abajo (la raíz cuadrada proviene de un caminante aleatorio). Sin embargo, esta estimación puede empeorar significativamente por dos razones:

  1. Ocurre frecuentemente que por las regularidades del cálculo, o por las características de la computadora, los errores de redondeo se acumulan preferentemente en una dirección, y en este caso el error total será del orden de $\sqrt{N} \epsilon_m$.

  2. Algunas operaciones pueden incrementar significativamente el error de redondeo en operaciones únicas. Esto puede, generalmente, suceder cuando se realiza la resta de dos números muy cercanos, dando un resultado cuyos únicos bits significativos son unos pocos de orden bajo, en los que los operandos difieren. Podría pensarse que tal substracción casual es poco probable, pero no siempre es así. Algunas expresiones matemáticas magnifican la probabilidad de ocurrencia en forma notable. Por ejemplo, en la solución de una ecuación cuadrática: $$ x = \frac{-b + \sqrt{b^2 -4 a c}}{2 a}$$ la suma se convierte susceptible de que haya un error de redondeo importante siembre que $b > 0$ y $|ac| \ll b^2$.

Propagación y error de comparación

La aritmética de punto flotante no es exacta. Valores como $0.1$ no pueden representarse en forma precisa usando puntos flotantes binarios IEEE, y la precisión limitada de los números de punto flotante significa que pequeños cambios en el orden de las operaciones o en la precisión de los resultados intermedios puede cambiar el resultado. Esto significa que comparar dos floats para ver si son iguales (o distintos) no es generalmente lo que buscamos. Por ejemplo:

In [ ]:
f = 0.1
suma = 0.0
for i in range(10):
    suma += f
prod = f * 10
print(suma, prod)
print(suma == prod)

Es sencillo ver lo que sucede aquí. Comenzamos con f que es una aproximación a $0.1$. Luego sumamos f diez veces, con un error de redondeo en cada paso, por lo que hay mucho espacio para que el error se incremente. Por otra parte, prod realiza una sola operación sobre f, por lo que hay menos oportunidades de incrementar el error de redondeo. Además ocurre que la conversión de $0.1$ a f redondea hacia arriba, mientras que la multiplicación por diez redondea hacia abajo, y algunas veces dos redondeos producen un acierto. Por lo tanto en prod tenemos la respuesta correcta por las razones equivocadas. O tal vez la respuesta equivocada, ya que no es realmente diez multiplicado por 0.1.

Dado entonces que las representaciones de los puntos flotantes son solo aproximadas, al comparar dos valores es necesario verificar si están suficientemente cerca.

Una estrategia muy utilizada son las comparaciones con épsilon. Por ejemplo:

In [ ]:
epsilon = 1E-12
abs(suma - prod) < epsilon

Podríamos estar tentados de usar $\epsilon_m$, lo cual funcionaría perfectamente para números de magnitud del orden de $1$, pero esto no sería suficiente al comparar números que sean órdenes de magnitud menores.

Esto se puede subsanar comparando con un épsilon relativo. Para ello hay que evaluar la diferencia entre dos números, y comparar esta difrencia con sus magnitudes. Con el propósito de obtener resultados consistentes, siempre se debe comparar la diferencia con el mayor de los dos números:

In [ ]:
def casiIgual(a, b, maxDifRel = 2.22E-16): 
    A = abs(a)
    B = abs(b)
    dif = abs(a - b)
    mayor = max(A, B)
    if dif < mayor * maxDifRel:
        return True
    else:
        return False
casiIgual(1.00001E-18, 1.000001E-18)

Aún comparando con un épsilon relativo, hay casos donde este criterio falla, y aparecen entonces criterios más sofiticados realizando comparaciones en representación binaria de punto flotante IEEE. Se puede leer más sobre esto aquí.

Consideraciones finales

Por lo general los errores de redondeo son muy pequeños cuando se trabaja con doble precisión y operaciones aisladas. Sin embargo hay situaciones donde pueden aparecer problemas con operaciones simples o al acumularse en cálculos repetidos.

Hay que tener en cuenta las siguientes consideraciones:

  1. La multiplicación y la división son operaciones "seguras"
  2. La suma y la resta son peligrosas, ya que cuando los operandos son de órdenes de magnitud diferentes, los dígitos del más pequeño se pierden.
  3. Esta pérdida de dígitos menos significativos puede ser inocua (cuando los dígitos perdidos también son insignificantes para el resultado final), o catastrófica (cuando el error se acumula y afecta al resultado en forma importante).
  4. Cuantas más operaciones se realice en un cálculo (especialmente con los métodos iterativos), más atención hay que darle a este problema.
  5. Un método de cálculo puede ser estable (tiende a reducir los errores de redondeo) o inestable (los errores de redondeo se amplifican). Suele haber soluciones estables o inestables a un problema.

El problema de la estabilidad numérica de los algoritmos se estudia en el área de la matemática denominada análisis numérico.

In [ ]: