Предотвращение переполнения при добавлении двух логарифмов

Question

Я использую добавление в уравнении пространства журналов, описанное в статье Wikipedia log probability, но я получаю недостаточное значение при вычислении exp очень больших отрицательных логарифмов. В результате моя программа выйдет из строя.

Примеры ввода: a = -2 и b = -1033.4391885529124.

Мой код, реализованный прямо из статьи Википедии, выглядит следующим образом:

double log_sum(double a, double b) 
{ 
    double min_ab = std::min(a, b); 
    a = std::max(a, b); 
    b = min_ab; 
    if (isinf(a) && isinf(b)) { 
    return -std::numeric_limits<double>::infinity(); 
    } else if (isinf(a)) { 
    return b; 
    } else if (isinf(b)) { 
    return a; 
    } else { 
    return a + log2(1 + exp2(b - a)); 
    } 
} 

Я придумал следующие идеи, но не может решить, что лучше:

• Перед оценкой проверьте наличие вне диапазона.
• Отключить (как-то) исключение и вымыть или закрепить результат после оценки
• Реализация пользовательских функций журнала и экспозамера, которые не генерируют исключений и автоматически очищают или зажимают результаты.
• Некоторые другие способы?

Кроме того, мне было бы интересно узнать, какой эффект имеет выбор базы логарифмов при вычислении. Я выбрал базовую пару, потому что считал, что другие базы журналов будут вычисляться от log_n(x) = log_2(x)/log_2(n) и будут страдать от потери точности из-за деления. Это верно?

Answer 1

По http://en.cppreference.com/w/cpp/numeric/math/exp:

Для IEEE-совместимый тип двойной, переполнение гарантируется, если 709,8 < Arg и сгущенного гарантируется, если аргумент < -708,4

Таким образом, вы не можете предотвратить нижний поток. Однако:

Если ошибка диапазона возникает из-за недостаточного потока, возвращается правильный результат (после округления).

Таким образом, не должно быть никаких программных сбоев - «просто» потеря точности.

Однако обратите внимание, что

1 + exp(n) 

потеряет точность намного раньше, то есть уже при п = -53. Это связано с тем, что следующее представимое число после 1.0 равно 1.0 + 2^-52.

Так потери точности из-за exp гораздо меньше, чем точность теряется при добавлении 1.0 + exp(...)

Answer 2

Проблема здесь точно вычисления выражения log(1+exp(x)) без промежуточных в/переполнения. К счастью, Мартин Маэчлер (один из основных разработчиков R) рассказывает, как это сделать в section 3 of this vignette.

Он использует естественные базовые функции: он должен быть способен перевести его на base-2 путем соответствующего масштабирования функций, но он использует функцию log1p в одной части, и я не знаю какой-либо математической библиотеки, которая снабжает вариант base-2.

Выбор базы вряд ли повлияет на точность (или производительность), и наиболее разумные математические библиотеки могут предоставить гарантии на 1-ульп для обеих функций (т. Е. У вас будет одно из двух значений с плавающей запятой ближе всего к точному ответу). Довольно распространенный подход состоит в том, чтобы разбить число с плавающей запятой на его базовый уровень 2 k и значение 1+f, такое как 1/sqrt(2) < 1+f < sqrt(2), а затем использовать полиномиальное приближение для вычисления log(1+f): из-за некоторых математических причуд (в основном, тот факт, что второй член ряда Тейлора можно представить точно) она оказывается более точным, чтобы сделать это в природной основе, а не основывают-2, так что типичная реализация будет выглядеть следующим образом:

log(x) = k*log2 + p(f) 
log2(x) = k + p(f)*invlog2 

(например, см log и log2 в openlibm), поэтому нет реальной выгоды от использования одного над другим.