Sabemos log_add, pero ¿cómo hacer log_subtract?

Question

Multiplicando dos números en el espacio de registro de los medios de adición:

log_multiply(x, y) = log(exp(x) * exp(y)) 
        = x + y 

Adición dos números en el espacio de registro significa que haces un registro agregado especial de operación:

log_add(x, y) = log(exp(x) + exp(y)) 

que se implementa en el siguiente código, de una manera que no requiere que tomemos las dos exponenciales (y perdamos la velocidad y la precisión del tiempo de ejecución):

double log_add(double x, double y) { 
    if(x == neginf) 
     return y; 
    if(y == neginf) 
     return x; 
    return max(x, y) + log1p(exp(-fabs(x - y))); 
    } 

(Here es otra.)

Pero aquí es la pregunta:

¿Hay un truco para que lo haga por sustracción así?

log_subtract(x, y) = log(exp(x) - exp(y)) 

sin tener que tomar los exponentes y perder precisión?

double log_subtract(double x, double y) { 
    // ? 
    } 

Answer 1

¿Qué tal

double log_subtract(double x, double y) { 
    if(x <= y) 
    // error!! computing the log of a negative number 
    if(y == neginf) 
    return x; 
    return x + log1p(-exp(y-x)); 
} 

eso es sólo basado en un cálculo rápido que hice ...

Answer 2

Las funciones de biblioteca para exp y log perder precisión para los valores extremos. log1p lo lleva a medio camino hasta allí, pero lo que necesita es una función que trate el error tanto para el registro como para las partes exp.

Ver este artículo: http://cran.r-project.org/web/packages/Rmpfr/vignettes/log1mexp-note.pdf

El título es "con precisión Computing log (1 - exp (- | a |))".

El artículo describe cómo fusionar de manera uniforme diferentes algoritmos para crear buenos límites de error para una mayor gama de entradas.