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

Question

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

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) {
    // ?
  }

Answer 1

Как насчет

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

Это просто на основе какой-то быстрой математики, которую я сделал ...

Answer 2

Библиотечные функции для exp и log теряют точность для экстремальных значений. log1p поможет вам на полпути, но вам нужна функция, которая обрабатывает ошибку как для журнала, так и для части exp.

См. Эту статью: http://cran.r -project.org / web / packages / Rmpfr / виньетки / log1mexp-note.pdf

Заголовок «Журнал точных вычислений (1 - exp (- | a |))».

В статье обсуждается, как безошибочно объединять различные алгоритмы для создания хороших границ ошибок для большего диапазона входных данных.