[MOAB-dev] Floating point exceptions in MOAB

[Jason Kraftcheck]

Tim Tautges wrote:
> I'd like to hear the case for why division by zero be allowed anywhere. 
> This was always considered a no-no in any analysis code I've worked on. 
> The proper way to handle this has always been to test for zero before
> dividing.  The case for calling finite() has more to do with overflow,
> I'd think, though I'm not certain.  I do *not* think that having MOAB
> control exception generation is the right approach.
> 

Proceeding with the division and calling finite() afterwards catches both
division by zero and overflow.  Either of these could trigger a floating
point exception depending on what kind of mucking with the floating point
control register was done by some other library.  In general, changing the
floating point behavior and not restoring it to the system default before
passing control to code outside of the module is a bug.  The IEEE spec
specifies NaN and +/-Inf values precisely for this use.

Anyway, if you really want to catch both division by zero and overflow w/out
ever generating a NaN, +/-Inf or floating point exception, then something
more computationally expensive needs to be done.  For example the following:

/**\brief Perform safe division (no overflow or div by zero).
 *
 * Do result = num/den iff it would not result in an overflow or div by zero
 *\param num  Numerator of division.
 *\param den  Denominator of division
 *\param result The result of the division if valid, zero otherwise.
 *\return false if the result of the division would be invalid (inf or nan),
 *        true otherwise.
 */
inline bool divide( double num, double den, double& result )
{
  const double fden = fabs(den);
    // NOTE: First comparison is necessary to avoid overflow in second.
    // NOTE: Comparison in second half of condition must be '<'
    //       rather than '<=' to correctly handle 0/0.
  if (fden >= 1 || fabs(num) < fden*std::numeric_limits<double>::max()) {
    result = num/den;
    return true;
  }
  else {
    result = 0.0;
    return false;
  }
}