See, my suspicion is of the fact that it is “justified.” There are all kinds of places where exact equality of floating point numbers is the only correct solution.
Right now, I am in the middle of implementing a rather complicated mathematical function full of tricubic and bicubic splines. I.e. the kind of place where you would expect to see exorbitant amounts of numerical roundoff. But even in a place like here, exact equality of floats can be important. The function is chock full of switching functions, which take the form
0 if x <= xmin
f(x) = { g(x) if xmin < x < xmax
1 if xmax <= x
where g(x) is a function that transitions smoothly from 0 to 1 as
x goes from xmin to xmax
For any stable molecule, these numbers are never in the switching region. And as a result, there’s an awful lot of things in the function that almost always take on the value of an integer… including the parameters to those tricubic splines I mentioned earlier.
Tricubic splines are produced by fitting curves to a set of known values at integer points, so when implementing this code, it is extremely useful to have a fast path in the splines for checking if the point is an integer (in which case we don’t need to evaluate a polynomial, and can simply return the known value).
The way that I do this:
impl TricubicSpline {
fn evaluate(&self, point: V3<f64>) -> (EvalKind, (f64, V3<f64>)) {
let indices = point.map(|x| x as usize);
if point == indices.map(|x| x as f64) {
// Fast path (i, j, and k are all integers)
unimplemented!() // look up the known values
} else {
// Slow path (fractional point)
unimplemented!() // evaluate a cubic polynomial in 3 variables
}
}
}
minimizes roundoff error, while the method used by people following dogma:
if (fabs(Nij-floor(Nij)) < TOL && fabs(Nji-floor(Nji)) < TOL
&& fabs(Nijconj-floor(Nijconj)) < TOL) {
// fast path
} else {
// slow path
}
needlessly introduces small discontinuities into the value and derivatives of the computed spline.
