If you choose any option that is not compliant with the standard, that means, exactly as you have noticed, that you claim that you are an expert in FP computations and you know how to write FP programs that will give the desired results even with FP operations that can behave in strange ways.

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.

I think what people really want is a sane floating point mode that throws exceptions instead of these Inf and NaN things, so if their program isn't debugged, they get a clear signal by way of its abnormal termination.

Trapping on math exceptions or otherwise pulling an emergency break doesn't really fit with modern hardware design if you want performance (for one thing, what do you do with the rest of your vector? do you just live with a pipeline flush?).

Adding targeted checks used sparingly can work, but probably requires a deeper understanding of the underlying numerical analysis of your problem. Generally just checking for NaN or Inf at the end is a better solution as these are in most cases absorbing states.

Regarding performance, the cost of a real division is about 3-4 orders worse performance than an if statement that is very consistent, but the usual way is to have fast/safe versions of functions, where you need performance and can deduce if something is always/never true, create/use the fast function but by default everything uses the slower/safer function.

    l <= x <= u

    l' <= x + b <= u'

    l = l' - b
    u = u' - b

That’s a bold claim.

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.

I can't think of any, haven't experienced any, not heard of any examples of it, so you're welcome to break my ignorance on the subject.

NaN and Inf are extremely valuable and useful encodings to represent real things. Particularly signed infinity is wonderful, as functions like exp behave correctly for inf (eg. exp(-inf)=0).

But there are many ways such a calculation could result in NaN besides a division. For example, intermediate sums could result in something like inf - inf. I find it rare that you can definitely conclude that a NaN means some specific thing.

    var jsonish = "...";
    var number = parseFloat(jsonish);
    if (!isNaN(number)) {
        // deserialized a number
    } else {
        // ...try deserializing as a boolean or string next...
    }

I was porting a C++ project to a certain platform - and that platform enabled a -ffast-math equivalent by default in Release (but not Debug) builds! This broke duktape, a JS engine said project embedded, in some nasty and subtle ways. Instead of storing a number/pointer/??? (8 bytes) + type tag (4? bytes) for each dynamically typed JS value, duktape can bit-pack values into a single 8 byte "double" value by storing object/string handles as NaN values - this isn't an uncommon trick for dynamically typed scripting stuff:

https://github.com/svaarala/duktape/blob/c3722054ea4a4e50f48...

Naturally, the -ffast-math equivalent broke isNaN checks, which caused random object/string handles to be mistakenly reinterpreted as "numbers" - but only in Release builds, for this one particular platform, in one rarely taken branch, so neither QA nor CI caught it, leading to hours of manufacturing a repro case, stepping through an absurd amount of code, and then finally looking at the default build rules and facepalming.

Cursing the platform vendor under my breath, I overrode the defaults to align with the defaults of every other config x platform combination we already had: no fast math. If you want those optimizations, use SSE-friendly NaN-avoiding intrinsics - or, if you must use the compiler flags, ensure you do so consistently across build configs and platforms, perhaps limited to a few TUs or modules if possible. This allows you to have a chance at using your Debug builds to debug the resulting "optimizations".

I'd suspect the same is true for NaN. These things were introduced for reasons. If you don't have a deep knowledge of why, don't assume everyone else is the same.

Edit: I found examples using a NaN and Inf intermediates to get a large speed increase [1,2].

[1] https://tavianator.com/2015/ray_box_nan.html

[2] https://tavianator.com/2011/ray_box.html

So they are useful. They are more useful than most people realize that simply assume they only mark errors.

And without them in your case, you would not be alerted that a calculation went wrong. So even there they are very useful. Not using them correctly is like ignoring file API errors in your programs - sure they are rare, but you need to understand and handle them if you want to make good software.

±Inf is a very natural output for many kinds of numerical codes.

Years ago i had some interpolation to do, and used the barycentric rational from Boost [1]. The resulting curve was all over the place. A natural cubic spline gave a much more subjectively sensible-looking fit.

[1] https://www.boost.org/doc/libs/1_77_0/libs/math/doc/html/mat...

> What is the rationale for liking rational functions?

They are cheap to evaluate, (relatively) easy to prove things about symbolically, and will efficiently approximate any reasonably smooth function.

That's against one of the rules of well-behaved floating point programming: Never test for equality.

You should never test for equality numbers which are known only approximately, which is the case for most FP numbers.

There are nonetheless cases when certain FP values are known exactly and it is OK to test them for equality.

In general the rule does not depend on number representation, it applies equally to floating-point numbers, fixed-point numbers, rational numbers and even large integers in some cases.

What counts is whether a number is known exactly or only approximately.

But the type of barycentric interpolation formula that leads to Inf/Inf at an interpolation point isn't one of those exceptions - it's an example which shows why the rule should be borne in mind.

When the input is very close to one of the interpolation points, the calculation becomes ill-conditioned, unsurprisingly as two values are diverging towards Inf, so the result can become incorrect as the input approaches the interpolation point, until it becomes Inf/Inf = Nan at some even closer, but non-zero, distance. Before reaching Nan it can also have an interval where the result is +/-Inf.

It depends on the interpolation points. But under some circumstances, it is necessary to recognise when the input is close to an interpolation point, and adjust the formula appropriately.

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".)

You are correct, it does not. Quite the opposite: the results are excellent when x doesn’t quite equal one of the interpolation nodes, because you end up with a value of S · y / S, where S = «very big but not close to overflowing float» and y = «expected output value».

The only potentially problematic case is if S · y overflows but S does not. But floating point numbers have a ton of dynamic range, so unless you are looking at points really, really close to zero, with very large outputs (say, x − x_j = 1e-200, y = 1e+200), there is no problem. And the failure case there is that you end up with Inf instead of 1e200 as an output.

It isn’t a problem in practice, and if you are worried about it as a problem in theory, you can take care in choice of interpolation nodes or do some scaling of the output values.

I just tried it to be sure, in C with "double", and organized the calculation the way it's described in the paper you cited, formula 2.1.

It failed due to catastrophic rational cancellation close to an interpolation point, depending on the values involved. An equality test was not enough.

Approaching the first interpolation point, before it transitioned from the correct result, there were some incorrect Inf results prior to it switching to Nan for an input not equal to any interpolation point.

This was trivial 1d interpolation between two points, where the points were (0.0 => 10.0, 1.0 => 20.0).

I wonder what we're doing differently, or whether we happened to test at different points. Here's my code:

    def interpolate1(f,x,xs,ws):
        p,s = 1,0
        for (xj,wj) in zip(xs,ws):
            p *= x-xj
            s += wj*f(xj)/(x-xj)
        return p*s

For your example we want f(x) = 10+10x, xs=[0,1], ws=[-1,1].

As I mentioned before, if you implement the calculation of the polynomial called l(x) in the paper -- the thing my code calls p -- in an unfortunate way, you can get cancellation problems, but for this case I'm not sure I see any way for that to happen.

I tried calculating l(x) in a different way that could conceivably cause cancellation issues, aimed at this particular interpolant:

    def foo(x):
        p = x*x-x
        s = -1*10/x + 1*20/(x-1)
        return p*s

What exact values of x are giving you what bad results? What language are you implementing this in, with what data types?*

  #include <stdio.h>
  #include <math.h>

  double interp(double x, double x1, double y1, double x2, double y2)
  {
      double node_poly = (x - x1) * (x - x2);
      double w1 = 1.0 / (x1 - x2);
      double w2 = 1.0 / (x2 - x1);
      double y = node_poly * ((w1 * y1) / (x - x1) + (w2 * y2) / (x - x2));
      return y;
  }

  int main()
  {
      for (int p = -300; p >= -330; p -= 1) {
          double x = pow(10, p);
          double y = interp(x, 0, 10, 1, 20);
          printf("x = %.6g, y = %.6g\n", x, y);
      }
      return 0;
  }

The formula you want is the “second barycentric formula”, (4.2) at the top of page 505 of https://people.maths.ox.ac.uk/trefethen/barycentric.pdf

Or formula (2.6) halfway down page 2 of https://people.maths.ox.ac.uk/trefethen/publication/PDF/2011...

Also cf. Higham (2004) “The numerical stability of barycentric Lagrange interpolation” https://www.maths.manchester.ac.uk/~higham/narep/narep440.pd...

If you want to implement the “first barycentric formula” instead you need to take some additional care to avoid overflow/underflow.

As the papers above note, the choice of interpolation nodes makes a big difference. If you have an arbitrary function to approximate on an arbitrary complex domain, take a look at https://arxiv.org/abs/1612.00337

But none of that should make a difference here. What's actually causing the dubious results jlokier reports isn't exactly proximity to the interpolation points -- I don't think it could happen if one of them weren't exactly zero. What's happening is that the very smallest values of x in jlokier's code are denormal numbers; that is, ones whose floating-point representation has less precision than usual.

Usually a floating-point number looks like 1.<stuff> x 2^n; for IEEE double precision, the <stuff> is 52 bits long. So if you start with 1, 1/2, 1/4, 1/8, and keep dividing by 2, you always have stuff=0 and n decreases by 1 each time. But what should you do once you reach the smallest value of n that will fit in the floating-point format? One answer is just to say that after that you get zero. But you can get a little extra dynamic range if instead you say that when n is as small as it's allowed to be, you can use 0.<stuff> instead of 1.<stuff>. These numbers have reduced precision for the sake of extra dynamic range, and they are called "denormal". Lots of things go subtly (or not so subtly) wrong when dealing with denormal numbers.

One thing that goes wrong is the interpolation formula we're discussing here. The formula looks like f(x) ~= l(x) . sum of wj f(xj) / (x-xj) where xj are the interpolation points (in this case 0 and 1), l(x) is the product of (x-xj), so in this case x(x-1), and wj = product (xj-xi) where xi ranges over all the other interpolation points besides xj, so w0 = -1 and w1 = +1. So we have f(x) ~= x(x-1) . [-10/x + 20/(x-1)]. In fact this is exact provided x is neither 0 or 1, because the thing we're approximating is a polynomial of degree <= (# interpolation points) - 1.

So, what I claimed before was: suppose x = xj + h where h is small, and suppose you compute l(x) as the product of its factors. Then we get l(x) = h times product (xj+h-xi), and the sum is wj f(xj) / h + the sum of the other terms; and if h is tiny then the other terms in the sum are tiny in comparison, and the factors of h and 1/h cancel and everything's fine.

But that doesn't happen when xj=0 and x (and hence h) is denormal -- because you can't pull the denormal trick at the large end of the range, so the reciprocal of a denormal floating-point number is floating-point infinity!

This isn't a cancellation problem. It isn't a matter of the formula going unstable when x is very close to one of the interpolation points. But it is a problem, and it never occurred to me until jlokier demonstrated it, because denormals are kinda obscure and I didn't think to think about what might happen if they get involved. It would be avoided if you made your interpolator check for values very close to one of the xj. But -- I think, but this stuff is treacherous -- it would also be avoided if you checked for infinities as well as NaNs on the output side, which would be more efficient (provided your inputs are usually sensible). Or -- I think, but this stuff is treacherous -- you could just check explicitly for denormal values of x and maybe of the xj, because if none of those is denormal this pathology can't arise.

  double interp(double x, double x1, double y1, double x2, double y2)
  {
      double w1 = 1.0 / (x1 - x2);
      double w2 = 1.0 / (x2 - x1);
      double num = (w1 * y1) / (x - x1) + (w2 * y2) / (x - x2);
      double denom = w1 / (x - x1) + w2 / (x - x2);
      return num / denom;
  }

The interval is small for parameters interp(x, 0, 10, 1, 20). But interp(x, 0, 1e50, 1, 2e50) shows the interval can be much larger, as it depends on the y values.

Both formulas return incorrect results in an input range larger than the denormals (1e50 times larger in the second example above). Plenty of "normal" inputs that can arise from other calculations. And these parameters are not unusual either.

So you can't ensure a correct result from either barycentric formula in floating point by just checking for denormal input.

I think you may be right that checking for +/-Inf result as well as NaN may be an adequate test, but I wouldn't assume it without a proof or some paper demonstrating that for sure.

The weights w1,w2 are -1 and +1, as before.

The numerator is -10^50/x - 10^250. The denominator is -1/x - 1. If x is large enough (which it is in the cases that give spurious infinities) the first term in each of these is enough larger that adding the second doesn't change it; so what you're really computing is (-10^50/x) / (-1/x). And if x is small enough that 10^50/x is infinite, then boom.

Something similar could happen even without numbers as large as 10^50. If the value at 0 were 3 instead, then values of x bigger than 1/3 the largest representable non-infinity would produce the same problem. Using f(0)=10^50 enlarges the range in which this particular overflow happens by a factor of 10^50, though.

The fact that this isn't an exotic catastrophic-cancellation problem doesn't, of course, detract from your original point (with which, for the avoidance of doubt, I agree): a calculation that can blow up for certain exact values is likely to misbehave for values that are merely close to those, so testing for the exact values is likely to be a very bad idea.

If I were actually writing something that does interpolation in this way, I would certainly not want to trust that checking the results for infs and nans was sufficient (without, as you say, either proving it or finding someone else who had -- but I suspect it's actually only almost true).

I think NaN and inf can be really useful in specific cases, but they probably shouldn't be allowed to propagate through your program, and shouldn't occur unless you're specifically making use of them.

The programmer should choose which is most appropriate for an application. Choosing to both ignore the exceptions and not generate a special value that can be examined later, that is seldom acceptable.

R goes even farther, and uses one specific NaN (out of the million-billions) to signal “not available”, while other NaNs are just NaNs.

Its properties, used properly, make code much more simple and sane.

Cumbersome but not that difficult: Range-check all input values.

https://godbolt.org/z/1qjxvW3KK

I found this out the hard way several years ago and still have my own implementations of these functions laying about in case I want to use -ffast-math.

    float
    unit_price (float volume, int count) {
       return volume / count;
    }

    if (count > 0)
        total += unit_price(volume, count);

    total += unit_price(volume, count);
    if (count == 0)
        printf("oops\n");

Unfortunately most compilers have too permissive default options.

C and C++ programs should always be compiled, at least during development, with strict options, e.g.:

-fsanitize=undefined -fsanitize-undefined-trap-on-error

or equivalent options, which eliminate all the horror stories with unexpected behavior caused by integer divide-by-zero or overflow, pointer overflow, out-of-bounds access and many others.

Because unlike integer operations the FP operations are standardized, there is no need for compiler options, the behavior for exceptions can be controlled from the program with the functions defined in <fenv.h> in C or <cfenv> in C++.

The workflow would then be:

    1. prove that the input is not null (if input == 0) yielding a "not null number"
    2. running your calculation (no check needed for multiplication or division, if you add something you have to reprove that the number is not null)

Easy.

"Allow optimizations for floating-point arithmetic that assume that arguments and results are not NaNs or +-Infs."

The obvious implication is that is that arguments and results are always finite and as a programmer you are responsible of guaranteeing the correct preconditions.

Certainly do not use finite-math-only when dealing with external data.

so... never?

I see it as both, as with many such flags¹n though as you say it does put the onus on the developer to avoid inputs the optimisation may cause problems for. You are saying “I promise that I'm being careful not to do things that will be affected by these parts of the standard, or won't blame you for undefined behaviour if I do, so you can skip those rules if it helps you speed things up”.

[1] like those that enable optimisations which can cause significant breakage if pointer aliasing is happening

The picking and choosing specific flags enable might be better in general. The flag `-funsafe-math-optimizations` seems to be useful as general-purpose, although as this article mentions, some of those optimizations need the more unsafe ones to become active.

I would say that `-ffast-math` is a dangerously safe-sounding flag, `-funsafe-fast-math` would be better IMO.

You technically make that promise by writing in C (except a narrow set of allowed conversions) but most compilers on standard settings do let you get away with it.

Do you mean `if-else if` or `if-else`? Because the second would work and would always short-circuit to else.

  if (x < 7)
    doSomething();
  else
    doSomethingElse();