3
$\begingroup$

I want to implement the function $f(x,a,b) = a + x(b-a)$ where all the inputs are floating point (doubles, say), such that (a) $f(0,a,b)=a$ exactly; (b) $f(1,a,b)=b$ exactly; (c) $f(x,a,b) \le f(y,a,b)$ whenever $x \le y$; and preferably (d) it is accurate (correct up to rounding) for $0 \le x \le 1$.

Implementing $f(x,a,b)=a+x(b-a)$ directly does not work because for example $f(1.0,-1.0,\operatorname{prev}(1.0)) = 1.0$ (where $\operatorname{prev}(a)$ is the floating point number before $a$). And $f(x,a,b)=b-(1-x)(b-a)$ has the same issue.

Now $$f(x,a,b) = (1-x)(a+x(b-a))+x(b-(1-x)(b-a))$$ has the first two properties.

  • Does it have property (c)?
  • How accurate is it?
  • Is there a more performant way to do this?
$\endgroup$
  • $\begingroup$ What values can $x$ take? The formula suggests it is used to interpolate between $a$ and $b$, that is, $x$ in $[0,1]$. Is that the case? $\endgroup$ – njuffa May 19 at 14:23
  • $\begingroup$ @njuffa Absolutely. I have updated the question. $\endgroup$ – Doris May 19 at 15:48
1
$\begingroup$

Here is a proposed alternative $$f(x,a,b) = \operatorname{ifelse}(x<0.5, a+x(b-a), b-(1-x)(b-a))$$ which also satisfies (a) and (b), but it is not accurate near $0.5$. It also satisfies (c) except maybe at $\operatorname{prev}(0.5)$ and $0.5$.

For example with double-precision arithmetic we compute $f(0.5,-1.0,\operatorname{prev}(1.0))$ to be -1.11e-16 but it is actually -5.55e-17, so out by a factor of 2.

It's worth noting that the implementation proposed in the question gets this example correct.

$\endgroup$
1
$\begingroup$

In general, you can exclude cases x ≥ 0.5 as Doris shows: Find an implementation g (x, a, b) which is good for 0 ≤ x ≤ 0.5. If x ≥ 0.5 then 1 - x is calculated exactly without rounding error, so let f (x, a, b) = g (x, a, b) if x ≤ 0.5, and f (x, a, b) = g (1-x, b, a) if x ≥ 0.5.

a + x (b - a) is reasonably accurate if b and a have the same sign. Problem is b, a having different signs (because then there is some x where f (x, a, b) should be very close to zero, and there the smallish absolute error turns into a huge relative error, again as Doris showed - except that will happen for any a, b with different sign for some 0 ≤ x ≤ 1).

So assume a, b have opposite sign. Let S = round (b - a), that is the result of floating-point arithmetic. It can be shown that if abs(b) ≥ abs(a) then s = ((S - b) + a) is calculated without rounding error, and therefore b - a = S - s (calculated as real numbers). But if abs(a) ≥ abs(b) then s = ((S + a) - b) is calculated without rounding error, and again b - a = S - s. All this improves the situation, because we can calculate (a + xS) - xs, which is better because we handled the rounding error in b - a.

If you did the calculation in double precision (53 bit mantissa), then you can round both x and S to single precision (24 bit mantissa) and you get two numbers whose predict is exact. Let $x = x_{hi} + x_{lo}$, $S = S_{hi} + S_{lo}$, then $xS = x_{hi}S_{hi} + x_{hi}S_{lo} + x_{lo}S_{hi} + x_{lo}S_{lo}$. The first of the products is calculated exactly without rounding error, and the others are small. So we can calculate a + xS with less rounding error:

$a + xS ≈ (a + x_{hi}S_{hi}) + x_{hi}S_{lo} + x_{lo}S_{hi} + x_{lo}S_{lo}$

$\endgroup$
  • $\begingroup$ Thanks for the excellent answer! I didn't get notified so missed it at the time. How do you prove that s is computed exactly? $\endgroup$ – Doris May 19 at 15:46
  • $\begingroup$ If you add or subtract two numbers, the result is always exact if the exponent of the result is not larger than the exponent of either number. Because of that, x - y is always calculated correctly if x/2 ≤ y ≤ 2x. $\endgroup$ – gnasher729 May 19 at 17:57
0
$\begingroup$

Linear interpolation is a deceptively simple problem, for which, based on a literature search, no comprehensive nor definitive (i.e. with proven properties and error bounds) treatment seems to have been published to date. My previous attempt to contact a French researcher who has worked on formal proofs for similar "trivial" problems did not yield a response. The following discussion is based on running various arithmetic arrangements with a large number (billions) of test cases. In terms of the basic operations used in such rewrites, one should consider fused multiply-add (FMA) along with additions and multiplications.

It should be noted that due to the issue of subtractive cancellation, accurate results cannot be achieve when $a$ and $b$ are of opposite sign without the use of some form of extended precision. The initial portion of this answer therefore assumes that $a$ and $b$ have like signs. Likewise it is assumed that none of the intermediate computation leads to underflow, which in my test framework is guaranteed by limiting the range of $a$ and $b$ appropriately. All variants discussed in the following handle the end cases $x = 0$ and $x = 1$ correctly.

(1 - x) * (a + x * (b - a)) + x * (b - (1 - x) * (b - a))

is the formula suggested in the question above. It is non-monotonic, that is, property (c) does not hold, as simple random test cases show. My tests suggest an error bound of $2.75$ ulp for this formula.

A question on accurate and monotonic linear interpolation very similar to the present one was previously asked on the companion site Math Stack Exchange. The accepted answer there proposes the following arrangement:

ifelse (x < 0.5, a + x * (b - a), b - (1 - x) * (b - a))

My tests indicate (but cannot prove) that this variant is monotonic, including in the vicinity of $0.5$. Error seems to be bounded by $2$ ulps. Personally, I have subsequently adopted a modified form of this, taking advantage of FMA:

ifelse (x <= 0.5, fma (b - a, x, a), fma (b - a, x - 1, b))

Note the slightly different guard condition. Again, this variant appears to be monotonic according to my tests, and the error bound seems to be $1.5$ ulps. A very fast version for situations where monotonicity doesn't matter is the following:

fma (b, x, fma (-a, x, a))

This provides an apparent error bound of $1$ ulp, but is easily shown to be non-monotonic via random test cases.

One way of addressing the issue of subtractive cancellation when the signs of $a$ and $b$ differ (which in the worst case leads to results that are off by a factor of $2$ as noted by Doris), is to perform the intermediate computations in double native precision, using techniques originally developed by Knuth and Dekker.

My own FMA-enhanced attempt at this is shown as C-family code below. Tests indicate it provides monotonicity and an error bound of $0.501$ ulp, that is, results are nearly correctly rounded. It might be more convenient and faster to use built-in higher precision for this approach if available, as in gnasher729's answer.

typedef struct {
    T lo, hi;
} T2;

/* compute x+y, with |x| >= |y| */
T2 add_x_ge_y (T2 x, T2 y)
{
    T2 z;
    T r, s, e;
    r = x.hi + y.hi;
    e = x.hi - r;
    s = ((e + y.hi) + y.lo) + x.lo;
    z.hi = e = r + s;
    z.lo = (r - e) + s;
    return z;
}

T2 aa, ax, bx, tt, rr;
/* -a*x */
ax.hi = -a * x;
ax.lo = fma (-a, x, -ax.hi);
/* b*x */
bx.hi = b * x;
bx.lo = fma (b, x, -bx.hi);
/* a - a*x, Since x <= 1, |a| >= |a*x| */
aa.hi = a;
aa.lo = 0;
tt = add_x_ge_y (aa, ax);
/* a - a*x + b*x */
rr = ifelse (fabs (tt.hi) >= fabs (bx.hi), add_x_ge_y (tt, bx), add_x_ge_y (bx, tt));
return rr.hi + rr.lo;
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.