Мы знаем log_add, но как сделать log_subtract?

Умножение двух чисел в лог-пространстве означает их добавление:

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

Добавление двух чисел в лог-пространстве означает, что вы выполняете специальную операцию добавления журнала:

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

который реализован в следующем коде, таким образом, что нам не требуется брать две экспоненты (и потерять скорость и точность выполнения):

  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) ));
  }

(Здесь - еще один.)

Но вот вопрос:

Есть ли уловка для его вычитания?

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

без необходимости брать экспоненты и потерять точность?

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

Ответы

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));
}