Rearranging equation to prevent denominator blowup in C

Question

I know the title mentions C, because this is a programming related problem, but I think this specific issue is more pure mathematics so I figured here would make more sense!

For a homework assignment, we have to implement an equation in C and Mathematica and evaluate it over [1:500]. The Mathematica code works, and the C code is expected to fail around input value 150.

The relevant parts of the equation are:

Where q = [1:500] and all other variables are in the range -1:1 except epsilon_W which is 10.

The problem with the equation in C is that the denominator of the first equation blows up very quickly, eventually coming out to "inf" in C. Later on in the evaluation, the numerator also becomes infinity, leading to C coming out with "NaN" (not a number).

Is there a way to rearrange this equation in order to prevent or help to tone down this issue? I have been working on it with some math buddies, and we suspect it has something to do with the ln(e) relationship and rules, but we cannot seem to find anything that works. It could also be, instead of a rearranging, a series that approximates the values? If anyone has any suggestions, we would love to hear!

Answer 1

You can use $$ \log \frac{e^a+e^b}2 = \log\left(e^a\frac{1+e^{b-a}}2\right) = a + \log\left(1+\frac{e^{b-a}-1}2\right) $$ which will avoid overflow if you choose $a$ (at runtime) to be the larger of the two exponents.

To protect against precision loss when $a$ and $b$ are close to each other, be sure to use the expm1() and log1p() functions instead of exp() and log() to evaluate the rightmost expression.

When $a$ is small and close to $-b$, the addition of $a$ in the above formula carries some additional risk of cancellation, and so it may be better to switch to $$ \log \frac{e^a+e^b}2 = \log\left(1+\frac{(e^{a+b}-1)-(e^a-1)(e^b-1)}2\right) $$ (again computed using expm1 and log1p) when $|a|$ and $|b|$ are both less than, say, $1/4$. This can be simplified further when $a=-b$ exactly, as is the case in your denominator.