Floating point exceptions in MOAB
MOAB divides by zero sometimes. I discovered this because a CGM dependency enabled floating point exceptions without my knowledge, and that caused my program to abort when a zero division occurred inside MOAB. The particular zero division I ran into occurs in GeomUtil.cpp, in the function segment_box_intersect(). (Only the KD tree code calls this function.) This zero division is not a bug: the result is checked using finite(), and discarded if infinite. Grepping the code for other occurrences of finite() suggests that this idiom exists in a handful of other places in MOAB. (Most of these places are geometry-related: KD trees, BSP trees, oriented boxes.) Code like this isn't inherently buggy, but it does lead to unexpected behavior when floating point exceptions are enabled. The library that silently enabled floating point exceptions was CGM, or rather, one of CGM's dependencies. Specifically, in a CGM built against Cubit 10.2 (with ACIS 16.0), one of the Cubit/ACIS shared libraries that CGM links to enables floating point exceptions. (I'm not sure which library it is specifically.) This does not occur when I build CGM against ACIS 16.7. So the situation is this: if I build CGM against Cubit 10.2, and then build MOAB against that CGM installation, certain MOAB functions may crash unless the calling program first disables floating point exceptions (e.g. by calling feholdexcept()). More generally, any program that enables floating point exceptions may encounter those exceptions when calling into MOAB. Is this mentioned in the MOAB documentation? ~S
On Wed, 2010-04-28 at 13:15 -0500, Steve Jackson wrote:
The library that silently enabled floating point exceptions was CGM, or rather, one of CGM's dependencies. Specifically, in a CGM built against Cubit 10.2 (with ACIS 16.0), one of the Cubit/ACIS shared libraries that CGM links to enables floating point exceptions. (I'm not sure which library it is specifically.) This does not occur when I build CGM against ACIS 16.7.
Incidentally, CGM tried to do this in iGeom_getPntRayIntsct but I removed it in r3716 because it was one layer of what caused crashes in that function. If we fix the configuration issue causing this, we could probably revert that change too. - Jim
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. - tim Steve Jackson wrote:
MOAB divides by zero sometimes. I discovered this because a CGM dependency enabled floating point exceptions without my knowledge, and that caused my program to abort when a zero division occurred inside MOAB.
The particular zero division I ran into occurs in GeomUtil.cpp, in the function segment_box_intersect(). (Only the KD tree code calls this function.) This zero division is not a bug: the result is checked using finite(), and discarded if infinite. Grepping the code for other occurrences of finite() suggests that this idiom exists in a handful of other places in MOAB. (Most of these places are geometry-related: KD trees, BSP trees, oriented boxes.) Code like this isn't inherently buggy, but it does lead to unexpected behavior when floating point exceptions are enabled.
The library that silently enabled floating point exceptions was CGM, or rather, one of CGM's dependencies. Specifically, in a CGM built against Cubit 10.2 (with ACIS 16.0), one of the Cubit/ACIS shared libraries that CGM links to enables floating point exceptions. (I'm not sure which library it is specifically.) This does not occur when I build CGM against ACIS 16.7.
So the situation is this: if I build CGM against Cubit 10.2, and then build MOAB against that CGM installation, certain MOAB functions may crash unless the calling program first disables floating point exceptions (e.g. by calling feholdexcept()). More generally, any program that enables floating point exceptions may encounter those exceptions when calling into MOAB. Is this mentioned in the MOAB documentation? ~S
-- ================================================================ "You will keep in perfect peace him whose mind is steadfast, because he trusts in you." Isaiah 26:3 Tim Tautges Argonne National Laboratory ([email protected]) (telecommuting from UW-Madison) phone: (608) 263-8485 1500 Engineering Dr. fax: (608) 263-4499 Madison, WI 53706
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; } }
participants (4)
-
James Porter -
Jason Kraftcheck -
Steve Jackson -
Tim Tautges