Optimizations Enabled by -ffast-Math

> So yes, that means that you become responsible to either guarantee that erroneous results do not matter or that you will take care to always check the ranges of input operands, as "okl" has already posted, to ensure that no overflows, underflows or undefined operations will happen.

This is why I dislike the fact that it removes simple methods for checking if values are NaN, for example. I do not find it ergonomic to disable these checks.

Yeah, I really like that Julia and LLVM allow applying it on a per-operation basis.

Because most of LLVM's backends don't allow for the same level of granularity, they do end up propagating some information more than I would like. For example, marking an operation as fast lets LLVM assume that it does not result in NaNs, letting nan checks get compiled away even though they themselves are not marked fast:

  julia> add_isnan(a,b) = isnan(@fastmath(a+b))
  add_isnan (generic function with 1 method)
  
  julia> @code_llvm debuginfo=:none add_isnan(1.2,2.3)
  define i8 @julia_add_isnan_597(double %0, double %1) #0 {
  top:
    ret i8 0
  }

meaning it remains more dangerous to use than it IMO should be. For this reason, LoopVectorization.jl does not apply "nonans".

It's worth being a bit more explicit about this.

The appearance of NaNs doesn't result from breaking the "never test for equality" rule. It happens when the point where you're evaluating the function just happens to exactly match one of the interpolation points.

But if you decide to deal with the NaN issue by (1) testing whether your input exactly matches any of the interpolation points, or (2) testing the result for NaN, then you're testing for equality or doing something very similar, and then you may have one of the problems the "never test for equality" is meant to prevent: if bad stuff happens when x exactly equals one of the x_j, then perhaps less-obvious bad stuff happens when x almost exactly equals one of the x_j, and an equality test won't catch that.

So, does it? Actually, I think not. Suppose x = xj+h with h very small. Then your interpolant is (h times something of normal size) times (1/h times something of normal size + smaller terms) and although one of those is very small and the other very large, ordinary floating-point arithmetic doesn't care about that. There's no catastrophic cancellation or anything. And floating point has much more dynamic range than precision, so to speak, so until h literally reaches zero you aren't running into trouble. (I guess there might be ways of organizing the calculation that do give you catastrophic cancellation, but if you do it the obvious way then that doesn't happen.)

So in this particular case you really could just test for equality or test for NaNs. I think.

(For some interesting investigation of other aspects of the numerical stability of this kind of interpolation, see https://people.maths.ox.ac.uk/trefethen/publication/PDF/2011.... But note that the authors of this paper are using the formula for extrapolation or, as they prefer to call it, "numerical analytic continuation".)

There are fast float intrinsics: https://doc.rust-lang.org/std/intrinsics/fn.fadd_fast.html

but better support dies in endless bikeshed of:

• People imagine enabling fast float by "scope", but there's no coherent way to specify that when it can mix with closures (even across crates) and math operators expand to std functions defined in a different scope.

• Type-based float config could work, but any proposal of just "fastfloat32" grows into a HomerMobile of "Float"

• Rust doesn't want to allow UB in safe code, and LLVM will cry UB if it sees Inf or Nan, but nobody wants compiler inserting div != 0 checks.

• You could wrap existing fast intrinsics in a newtype, except newtypes don't support literals, and

As a scientist, it depends on what you mean by wrong. It's nice to push an array of numbers through some array operation. If some of the outputs look like NaN or Inf then that tells me something went wrong in the operation and to take a closer look. If some optimizer was told that NaN or Inf couldn't happen, meaning that they wouldn't be generated, then some of this output could be pure nonsense and I wouldn't notice. As NaNs propagate through the calculation they're extremely helpful to indicate something went wrong somewhere in the call chain.

NaN is also very useful itself as a data missing value in an array. If you have some regularly sampled data, with some items missing, it makes a lot of things easier to populate the holes with values which guarantee they won't produce misleading output.