Aritmética en punto flotante

La idea

Un ordenador no trabaja con números reales sino con un subconjunto finito: los números en punto flotante. El estándar IEEE 754 de doble precisión (el double de C, Java o JavaScript) usa 64 bits: 1 de signo, 11 de exponente (almacenado con un sesgo de 1023) y 52 de mantisa, a los que se añade un 1 implícito. El valor es (−1)ˢ · 1,m · 2ᵉ⁻¹⁰²³. Los exponentes 0 y 2047 se reservan para los subnormales y el cero, y para los infinitos y NaN.

Con 53 bits significativos solo se representan exactamente las fracciones cuyo denominador es potencia de 2. El famoso 0,1 es en realidad 0,1000000000000000055511151…, y por eso 0,1 + 0,2 da 0,30000000000000004 y no es igual a 0,3. El épsilon de máquina, ε = 2⁻⁵² ≈ 2,22 · 10⁻¹⁶, es la distancia entre 1 y el siguiente double: mide la precisión relativa (unas 16 cifras decimales). Los doubles no están espaciados uniformemente: entre 1 y 2 hay 2⁵² números, y la misma cantidad entre 2⁵² y 2⁵³, donde ya no caben decimales y 2⁵³ + 1 se redondea a 2⁵³.

El peligro principal es la cancelación catastrófica: al restar dos números casi iguales, los dígitos comunes se anulan y quedan solo los errores de redondeo. Calcular √(x + 1) − √x para x = 10¹² pierde la mitad de las cifras; reescribirlo como 1/(√(x + 1) + √x) las recupera todas. También importa el orden de las sumas: sumar una serie de mayor a menor va perdiendo los términos pequeños, mientras que sumarla de menor a mayor, o con el algoritmo de Kahan, es mucho más exacto. La regla práctica: nunca comparar flotantes con ==, sino con una tolerancia relativa.

Fórmulas clave

• Valor de un double normalx = (−1)ˢ · (1 + m/2⁵²) · 2ᵉ⁻¹⁰²³, s: 1 bit, e: 11 bits, m: 52 bits
• Épsilon de máquinaε = 2⁻⁵² ≈ 2,220446 · 10⁻¹⁶
• Espaciado (ulp) cerca de xulp(x) = 2^(⌊log₂|x|⌋ − 52)
• Error de redondeofl(x ∘ y) = (x ∘ y)(1 + δ), |δ| ≤ ε/2
• Cancelación evitada√(x+1) − √x = 1 / (√(x+1) + √x)
• Suma de Kahany = xᵢ − c; t = s + y; c = (t − s) − y; s = t

Dónde se aplica

Los sistemas bancarios nunca guardan dinero en doubles: usan enteros de céntimos o tipos decimales, porque acumular 0,1 mil veces no da 100. El fallo del misil Patriot en 1991 se debió a que 0,1 s en binario, acumulado durante 100 horas, desviaba el reloj un tercio de segundo. En ciencia de datos, calcular la varianza con la fórmula Σx² − n·x̄² sufre cancelación catastrófica con datos grandes de poca dispersión (por eso se usa el algoritmo de Welford), y en aprendizaje automático el truco log-sum-exp evita desbordamientos al calcular softmax. Conocer los bits explica por qué JSON no puede representar 2⁵³ + 1 y por qué los identificadores largos se envían como texto.

Con programación (Python)

import math

def campos(x):
    """Signo, exponente y mantisa (en bits) de un double a partir de float.hex()"""
    h = x.hex()                       # p. ej. '0x1.999999999999ap-4'
    signo = '1' if h.startswith('-') else '0'
    mant_hex, exp = h.lstrip('-')[4:].split('p')
    mantisa = bin(int(mant_hex, 16))[2:].zfill(52)
    exponente = bin(int(exp) + 1023)[2:].zfill(11)
    return signo, exponente, mantisa

s, e, m = campos(0.1)
print("0.1 -> signo", s, "exponente", e, "(=", int(e, 2) - 1023, ") mantisa", m)
print("0.1 con 25 decimales:", f"{0.1:.25f}")
print("0.1 + 0.2 == 0.3 ->", 0.1 + 0.2 == 0.3, "  0.1 + 0.2 =", repr(0.1 + 0.2))

eps = 1.0
while 1.0 + eps / 2 > 1.0:
    eps /= 2
print("épsilon de máquina =", eps, "= 2^", round(math.log2(eps)))

x = 1e12
print("sqrt(x+1) - sqrt(x): ingenuo", math.sqrt(x + 1) - math.sqrt(x), " estable", 1 / (math.sqrt(x + 1) + math.sqrt(x)))

n = 20000
adelante = sum(1 / (k * k) for k in range(1, n + 1))
atras = sum(1 / (k * k) for k in range(n, 0, -1))
print("suma adelante - suma atrás =", adelante - atras)
print("2**53 + 1 ==", float(2**53) + 1)

La versión en JavaScript, editable y ejecutable en el navegador, está en la página interactiva del tema.

Temas relacionados

Método de bisección · Complejidad y notación O grande · Polinomios de Taylor