How to test for overflow when multiplying floats

Question

I am trying to implement a 3-term recurrence relation: $$ p_{n+1} = ap_n + bp_{n-1} $$ This can be implemented as

  p0 = P0 // starting value p_0
  p1 = P1 // starting value p_1
  loop n = 2,...,N
    pn = a*p1 + b*p0
    p0 = p1
    p1 = pn

Depending on $(a,b)$ and also the stopping point $N$, this can sometimes overflow, and I would like to detect such a failure. I believe it is bad form to simply let the loop run to completion and then use something like isfinite to check the final value.

A more elegant method would be to somehow test whether the term a*p1 + b*p0 will exceed the floating point maximum (DBL_MAX = 1.7976931348623157e+308) prior to computing it. I could not find much information about this online, so does anyone have suggestions?

One idea I had, was suppose we are just dealing with the first term $a p_n$. Then we could test: $$ \log{|a p_n|} > \log{\textrm{DBL_MAX}} $$ or $$ \log{|a|} + \log{|p_n}| > \log{\text{DBL_MAX}} $$ It seems that this test would work even if $a p_n$ is not representable as a floating point number. I'm not sure how to extend this test to the general case $a p_n + b p_{n-1}$, or even if this test will account for all possible failure cases in double precision.

Answer 1

I want to further emphasize and underscore gnasher729's comment which is completely correct. NaN and Infinity were added exactly so you wouldn't need to do such checks all the time. That is, they were added to you wouldn't need to do exactly the thing you are suggesting doing. One of the consequences is that it often is the right thing to do to just blindly continue a calculation, checking whether it has overflowed or produced NaN only when you need to make a decision based on the result.

There are at least a few reasons it is this way. First, as you are starting to see, it can complicate your code quite a bit to attempt to preemptively detect potential errors. Second, such checking code tends to significantly slow down computation¹. This is especially bad as the vast majority of the time you will not overflow, and so these checks are just wastes of time. Adding early abort code is then optimizing for a very rare case at the significant expense of the vastly more common case. Third, many aspects of IEEE754 have been carefully designed to behave sensibly even in "error" cases. Often this allows correct results despite overflows in intermediate calculations and allows even more error checking code to be eliminated as well as expanding the domain of operations. For example, dividing a finite number by Infinity doesn't cause an error but instead produces a (signed) 0, which is likely the desired result.

So, to reiterate gnasher729's comment, the way you test for overflow when multiplying IEEE754 floating point numbers is by multiplying them and then checking whether they've overflowed. And you shouldn't even do that that often. You should only check for overflow when you need to make a code branch based on the decision², and otherwise you should just let the error values propagate.

¹ This especially true for SIMD architectures such as GPUs. GPUs are much happier doing the same thing to a bunch of data in parallel than doing different things in data dependent ways.

² These code branches may be hidden, e.g. while calling a print routine may not require a branch in your code, there is one in the print routine to know whether to display "Infinity" or not. Therefore, you should check for overflow before printing assuming you want a behavior other than printing "Infinity", et cetera.

Answer 2

I found a solution, though I don't know if it is optimal. One could apply the following test:

if (fabs(p1) > DBL_MAX / fabs(a) || fabs(p0) > DBL_MAX / fabs(b)) { // overflow }

This would catch the error before it actually occurs. Though as gnasher says in his comments, one could do the operation, and then check afterward if the result is inf or nan. I don't know which method is considered better.

One consideration could be checking how long it takes to apply the above test vs calling isfinite. I don't know if one method is more optimal than the other.