Prevent underflow when adding two logarithms

Question

I am using the addition in log space equation described on the Wikipedia log probability article, but I am getting underflow when computing the exp of very large, negative, logarithms. As a result, my program crashes.

Example inputs are a = -2 and b = -1033.4391885529124.

My code, implemented straight from the Wikipedia article, looks like this:

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

I've come up with the following ideas, but can't decide which is best:

• Check for out-of-range inputs before evaluation.
• Disable (somehow) the exception, and flush or clamp the output after evaluation
• Implement custom log and exp functions that do not throw exceptions and automatically flush or clamp the results.
• Some other ways?

Additionally, I'd be interested to know what effect the choice of the logarithm base has on the computation. I chose base two because I believed that other log bases would be calculated from log_n(x) = log_2(x) / log_2(n), and would suffer from precision loss due to the division. Is that correct?

Answer 1

According to http://en.cppreference.com/w/cpp/numeric/math/exp:

For IEEE-compatible type double, overflow is guaranteed if 709.8 < arg, and underflow is guaranteed if arg < -708.4

So you can't prevent an underflow. However:

If a range error occurs due to underflow, the correct result (after rounding) is returned.

So there shouldn't be any program crash - "just" a loss of precision.

However, notice that

1 + exp(n)

will loose precision much sooner, i.e. already at n = -53. This is because the next representable number after 1.0 is 1.0 + 2^-52.

So loss of precision due to exp is far less than the precision lost when adding 1.0 + exp(...)