Короче говоря, используйте следующее выражение:

fmax(la, lb) + log1p(exp(-fabs(lb - la)))

где я использую la и lb в качестве переменных, хранящих log(a) и log(b) соответственно, а имена функций взяты из языка C math.h. Большинство других языков также будут иметь эти функции, но они могут называться по-другому, например. abs и max вместо fabs и fmax. (Примечание: это соглашения, которые я буду использовать в этом ответе.)

Говоря о функциях, У Марка Дикинсона есть хорошая мысль : вы можете проверить, есть ли у вас доступ к той, которая сделает это за вас напрямую. Например, если вы используете Python и NumPy, у вас есть доступ к logaddexp , а в SciPy — logsumexp.

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

Более детально

Нет такого простого правила, как правила умножения и деления, но есть математическое тождество, которое поможет:
log(a + b) = log(a) + log(1 + b/a)

Мы можем немного поиграть с этим идентификатором, чтобы получить следующее:

log(a + b) = log(a) + log(1 + b/a)
           = log(a) + log(1 + exp(log(b) - log(a)))
           = la + log(1 + exp(lb - la))

Здесь еще есть проблема. Если la намного больше, чем lb, у вас будет 1 + 0.000...000something внутри log. В мантиссе с плавающей запятой недостаточно цифр для хранения something, поэтому вы просто получите log(1), полностью потеряв lb.

К счастью, в стандартной библиотеке большинства языков программирования есть функция для решения именно этой задачи, log1p, которая вычисляет логарифм единицы плюс ее аргумент. То есть log1p(x) возвращает log(1 + x), но точно для очень маленьких x.

Итак, теперь у нас есть:

log(a + b) = la + log1p(exp(lb - la))

Мы почти на месте. Есть еще одна вещь, которую следует учитывать. В общем, вы хотите, чтобы la было больше, чем lb. Это не всегда будет иметь значение, но иногда это даст вам дополнительную точность.* Если разница между lb и la действительно велика, это убережет вас от переполнения в exp(lb - la). В самом крайнем случае расчет работает, когда lb равно отрицательной бесконечности (т.е. b равно 0), но не когда la равно.

Иногда вы будете знать, какой из них больше, и вы можете просто использовать его как la. Но когда это может быть любой из них, вы можете использовать максимальное и абсолютное значение, чтобы обойти это:

log(a + b) = fmax(la, lb) + log1p(exp(-fabs(lb - la)))

Сумма коллекции

Если вам нужно взять сумму более чем двух чисел, мы можем получить расширенную версию тождества выше:

log(a[0] + a[1] + a[2] + ...)
    = log(a[0] * (a[0]/a[0] + a[1]/a[0] + a[2]/a[0] + ...))
    = log(a[0]) + log(a[0]/a[0] + a[1]/a[0] + a[2]/a[0] + ...)
    = la[0] + log(1 + exp(la[1]-la[0]) + exp(la[2]-la[0]) + ...)

Мы захотим использовать те же приемы, что и при сложении двух чисел. Таким образом, мы получаем наиболее точный ответ и максимально избегаем переполнения и недополнения. Во-первых, log1p:

log(a[0] + a[1] + a[2] + ...)
    = la[0] + log1p(exp(la[1]-la[0]) + exp(la[2]-la[0]) + ...)

Другое соображение заключается в том, какой операнд вытащить перед log1p. До сих пор я использовал la[0] для демонстрации, но вы хотите использовать тот, который лучше. Это по тем же причинам, по которым мы хотели la > lb при добавлении двух чисел. Например, если бы la[1] было наибольшим, вы бы сделали такой расчет:

log(a[0] + a[1] + a[2] + ...)
    = la[1] + log1p(exp(la[0]-la[1]) + exp(la[2]-la[1]) + ...)

Поместив это в правильный код, он будет выглядеть примерно так (это C, но он должен хорошо переводиться на другие языки):

double log_sum(double la[], int num_elements)
{
    // Assume index_of_max() finds the maximum element
    // in the array and returns its index
    int idx = index_of_max(la, num_elements);

    double sum_exp = 0;
    for (int i = 0; i < num_elements; i++) {
        if (i == idx) {
            continue;
        }
        sum_exp += exp(la[i] - la[idx]);
    }

    return la[idx] + log1p(sum_exp);
}

Вычисление разницы в пространстве журнала

Это не было частью вопроса, но, тем не менее, это может быть полезно: вычитание в пространстве журнала может быть выполнено аналогичным образом. Основная формула такова:

log(a - b) = la + log(1 - exp(lb - la))

Обратите внимание, что это все еще предполагает, что la больше, чем lb, но для вычитания это еще важнее. Если la меньше lb, вы логарифмируете отрицательное число!

Как и в случае сложения, у этого есть проблема с точностью, которую можно исправить с помощью специализированных функций, но, как оказалось, есть два способа. Один использует ту же функцию log1p, что и выше, но другой использует expm1, где expm1(x) возвращает exp(x) - 1. Вот оба способа:

log(a - b) = la + log1p(-exp(lb - la))
log(a - b) = la + log(-expm1(lb - la))

Какой из них вы должны использовать, зависит от значения -(lb - la). Первый более точен, когда -(lb - la) больше примерно 0,693 (т. е. log(2)), а второй более точен, когда он меньше. Подробнее о том, почему это так и откуда взялось log(2), см. в этой заметке из проекта R, в которой оцениваются два метода.

Конечный результат будет выглядеть так:

(lb - la < -0.693) ? la + log1p(-exp(lb - la)) : la + log(-expm1(lb - la))

или в форме функции:

double log_diff(double la, double lb)
{
    if (lb - la < -0.693) {
        return la + log1p(-exp(lb - la));
    } else {
        return la + log(-expm1(lb - la));
    }
}

* Есть что-то вроде приятного места для этого. Когда разница между la и lb невелика, ответ будет точным в любом случае. Когда разница слишком велика, результат всегда будет равен большему из двух, поскольку числа с плавающей запятой не обладают достаточной точностью. Но когда разница в самый раз, вы получите лучшую точность, когда la больше.

Предположим без ограничения общности, что b<a.

log(a+b) = log(a) + log(1 + b/a)
         = log(a) + b/a - 1/2(b/a)^2 + 1/3(b/a)^3 etc.

Члены, включающие b / a, могут быть вычислены в терминах известных величин, поскольку

b/a = exp(log(b) - log(a))

Если при вычислении b/a происходит потеря значимости, то log(a+b) = log(a) с машинной точностью.