Hacker News Board Hacker News, browsed like a classic message board
Live data from Hacker News

Optimizations Enabled by -ffast-Math

kristerw.github.io

Post reply on HN
Print view View original

111–117 of 117 posts

Re: Optimizations Enabled by -ffast-Math

#111

Post by jlokier »

gjm11 wrote: post #109
jacobolus wrote: post #98

Earlier quoted context omitted.

You implemented the wrong formula. You don’t want to explicitly calculate "node poly". 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 interpolati…

There are advantages to using either version of the formula. If you use the "first" then you do need to take care about overflow/underflow when the number of interpolation points is large, but that isn't an issue here. The "first" behaves better than the "second" for extrapolation, as discussed in the paper by Trefethen et al. And yes, the choice of interpolation points is important. But none of that should make a di…

Let's say we switch to the "second barycentric formula":

  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;
  }
This still gives incorrect +Inf results in an interval around the first interpolation point.

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.

Re: Optimizations Enabled by -ffast-Math

#112

Post by gjm11 »

jlokier wrote: post #111
gjm11 wrote: post #109

Earlier quoted context omitted.

There are advantages to using either version of the formula. If you use the "first" then you do need to take care about overflow/underflow when the number of interpolation points is large, but that isn't an issue here. The "first" behaves better than the "second" for extrapolation, as discussed in the paper by Trefethen et al. And yes, the choice of interpolation points is important. But none of that should make a di…

Let's say we switch to the "second barycentric formula": 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; } This still gives incorrect +Inf results in an interval around the first interpolation point. The interval is…

Nice catch on using enormous values of y, where indeed it isn't just denormal values of x that cause trouble. Still not a catastrophic-cancellation problem, though. Let's look more closely. Suppose x is very very close to 0. (In particular, close enough that x-1 is indistinguishable from -1. This is not a source of trouble in this case.)

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

Re: Optimizations Enabled by -ffast-Math

#113

Post by jacobolus »

jlokier wrote: post #111
gjm11 wrote: post #109

Earlier quoted context omitted.

There are advantages to using either version of the formula. If you use the "first" then you do need to take care about overflow/underflow when the number of interpolation points is large, but that isn't an issue here. The "first" behaves better than the "second" for extrapolation, as discussed in the paper by Trefethen et al. And yes, the choice of interpolation points is important. But none of that should make a di…

Let's say we switch to the "second barycentric formula": 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; } This still gives incorrect +Inf results in an interval around the first interpolation point. The interval is…

If x is very small (1e+100 or something) you should rescale it.

Re: Optimizations Enabled by -ffast-Math

#114

Post by [deleted] »

rocqua wrote: post #11

I found the following note for -ffinite-math-only and -fno-signed-zeros quite worrying: The program may behave in strange ways (such as not evaluating either the true or false part of an if-statement) if calculations produce Inf, NaN, or -0.0 when these flags are used. I always thought that -ffast-math was telling the compiler. "I do not care about floating point standards compliance, and I do not rely on it. So opti…

[deleted]

Re: Optimizations Enabled by -ffast-Math

#115

Post by nwallin »

heavenlyblue wrote: post #35
rocqua wrote: post #11

I found the following note for -ffinite-math-only and -fno-signed-zeros quite worrying: The program may behave in strange ways (such as not evaluating either the true or false part of an if-statement) if calculations produce Inf, NaN, or -0.0 when these flags are used. I always thought that -ffast-math was telling the compiler. "I do not care about floating point standards compliance, and I do not rely on it. So opti…

> such as not evaluating either the true or false part of an if-statement Do you mean `if-else if` or `if-else`? Because the second would work and would always short-circuit to else.

No. Imagine you have this:

  if (x
With -ffast-math, the compiler is allowed to call neither doSomething() nor doSomethingElse() if x is nan. Which IMHO is not great.

Re: Optimizations Enabled by -ffast-Math

#116

Post by marcan_42 »

adrian_b wrote: post #44
marcan_42 wrote: post #21

Earlier quoted context omitted.

Flushing denormals to zero only matters if your calculations are already running into the lower end of floating-point exponents (and even with denormals, if they're doing that, they're going to run into lost precision anyway sooner or later). The useful thing denormals do is make the loss of precision at that point gradual, instead of sudden. But you're still losing precision, and a few orders of magnitude later you'…

Your arguments are correct, but the conclusion does not result from them. If we assume that underflows happen in your program and this, as you say, is a sign that greater problems will be caused by that, then you must not enable flush-to-zero, but you must enable trap-on-underflow, to see where underflows happen and to investigate the reason and maybe rearrange your formulas to avoid the too small results. Flush-to-z…

It is normal to get values close to zero in a correctly implemented floating point algorithm - what matters is that, at that point, you don't rely on their dynamic range or precision. Trap-on-underflow would produce many false positives for algorithms which would be perfectly happy having those numbers be treated as zero.

This is not uncommon - it is completely expected with, say, anything implementing an exponential decay! A slow exponential decay could run into subnormals for large spans of time, yet be fine having those treated as zero.

Re: Optimizations Enabled by -ffast-Math

#117

Post by kazinator »

celrod wrote: post #71
adrian_b wrote: post #15

Earlier quoted context omitted.

The whole rationale for the standard for floating-point operations is to specify the FP operations with such properties that a naive programmer will be able to write programs which will behave as expected. 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 gi…

> 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 t…

If NaN and Inf values can still happen, but you don't have reliable ways of working with them, that is poor.

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.

Post reply on HN