Explicar este código para calcular $\log(1+x)$

Question

Es bien sabido que usted necesita tomar el cuidado al escribir una función para calcular un $\log(1+x)$ al $x$ es pequeña. Porque de punto flotante de redondeo, $1+x$ tienen menos precisión que las $x$, lo que puede traducirse en grandes relativa de error en la computación de la $\log(1+x)$. En librerías numéricas que vemos a menudo una función log1p(), que calcula el registro de una, además de su argumento, por esta razón.

Sin embargo, hace poco estuve intrigado por esta implementación de log1p:

def log1p(x):
    y = 1 + x
    z = y - 1
    if z == 0:
        return x
    else:
        return x * (log(y) / z)

¿Por qué es x * (log(y) / z) una mejor aproximación a $\log(1+x)$ que sólo log(y)?

Answer 1

La respuesta implica una sutileza de números de punto flotante. Voy a usar de punto flotante decimal en esta respuesta, pero se aplica igualmente bien a binaria de punto flotante.

Considerar decimal aritmética de punto flotante con cuatro lugares de precisión, y vamos a

$$x = 0.0001234$$

A continuación, en virtud de punto flotante, además de a $\oplus$, tenemos

$$y = 1 \oplus x = 1.0001$$ y $$z = y \ominus 1 = 0.0001000$$

Si ahora nos denotan la perdida de precisión en $x$$s = 0.0000234$, y la parte restante por $\bar{x}$, entonces podemos escribir

$$x = \bar{x} + s$$ $$y = 1 + \bar{x}$$ $$z = \bar{x}$$

Ahora, el valor exacto de $\log(1+x)$ es

$$\log(1+x) = \log(1+\bar{x}+s) = \bar{x}+s + O(\bar{x}^2) = 0.0001234$$

Si calculamos el $\log(y)$ entonces tenemos

$$\log(1+\bar{x}) = \bar{x} + O(\bar{x}^2) = .0001000$$

Por otra parte, si calculamos el $x \times (\log(y)/z)$ entonces tenemos

$$(\bar{x}+s) \otimes (\log(1+ \bar{x}) \div \bar{x}) = (\bar{x}+s)(\bar{x}\div \bar{x}) = \bar{x}+s = 0.0001234$$

así que mantener los dígitos de precisión que se habrían perdido sin esta corrección.