4 ms·

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

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.

gjm11 5y ago
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).
jacobolus 5y ago
If x is very small (<1e-100 or something) you can alternatively evaluate a linear approximation accurate near 0. If y is going to be enormous (>1e+100 or something) you should rescale it.

gjm11 5y ago

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

jacobolus 5y ago

If x is very small (<1e-100 or something) you can alternatively evaluate a linear approximation accurate near 0. If y is going to be enormous (>1e+100 or something) you should rescale it.