4
$\begingroup$

Given two (here and below: single-precision, IEEE 32-bit floats) normalized floating-point numbers $x, y$ (perhaps of reasonable range: my counterexamples don't have unusual magnitudes), and another floating-point number $0\leq t\leq 1$, I would like to evaluate the linear interpolant $$ z_t = x(1-t) + yt $$ in such a way that: $$ \begin{gathered} x \leq z_t \leq y,\\ z_0 = x, \qquad z_1 = y, \end{gathered} $$ as well as $$ z_t \leq z_s \quad\text{whenever}\quad 0\leq t \leq s \leq 1. $$

The direct formula, $z_t = x(1-t)+yt$ fails. For example, with $x=2.0000076$, $y=2.9007149$, we have in single-precision arithmetic $$ \mathrm{fl}(z_t) < x, \qquad \text{for}\quad2.9802326\times 10^{-8}\leq t<4.109652\times 10^{-8}. $$

The obvious other formula $z_t = x + (y-x)t$ satisfies $z_t \geq x$, but sometimes $z_1 > y$ for $$ x=1.0000609, \quad y=3.9999998, \quad z_1 = 4.0. $$

I considered $$ z_t = \mathrm{fma}(x, 1-t, y t), $$ and it has counterexamples x=y=1.9210567f0, t=0.10707888f0, z_t = 1.9210566f0 < x and x=-1.4928777f-36, y=-1.4825014f-36, t=2.9802322f-8, z_t=-1.4928778f-36 < x.

If I allow using if-statements, then simply evaluating one of $$ x+(y-x)t\quad(t\leq\tfrac12), \qquad y - (y-x)(1-t)\quad(t\geq\tfrac12) $$ works as defined above, but this is quite inefficient because of the branching on $t$. There is also sometimes an unexpected discontinuity at $t=\frac12$: $$ \begin{gathered} x = 0.9921611, \quad y = 4.0\\ x+\tfrac12(y-x) = 2.4960806, \qquad y-\tfrac12 (y-x) = 2.4960804. \end{gathered} $$

Does there exist an efficient simple formula that gives the linear interpolant with the standard mathematical properties?

$\endgroup$
  • $\begingroup$ Closely related question and possibly useful answer on Math Stack Exchange $\endgroup$ – njuffa Aug 7 '17 at 18:41
  • $\begingroup$ @njuffa I don't see any conclusive resolution there, but thanks for pointing to it. As a side note, they seemed to be using random numbers+exhaustive search for testing, but that can be a little unreliable—I think a good SMT solver like Z3 is a better way to go these days. $\endgroup$ – Kirill Aug 9 '17 at 20:50
  • $\begingroup$ Few people are equipped to craft a formal mechanically-assisted proof. I am not aware of a mathematically proven solution in the literature (publications by S. Boldo would be one obvious place to look), so you might need to tackle a proof yourself if one is required. Pedro's suggested method showed no monotonicity failures in my testing; that's good enough for me. Note that a two-path solution need not be inefficient, depending on processor architecture compilers can predicate code or use select / conditional-move instructions to remove the branch. $\endgroup$ – njuffa Aug 10 '17 at 3:54
3
$\begingroup$

The following link might give some insight https://fgiesen.wordpress.com/2012/08/15/linear-interpolation-past-present-and-future/

The suggested variant is: lerp_3(t, a, b) = fma(t, b, fnms(t, a, a))

(fnms is fused negative multiply subract, so it can be reexpressed as fma(-t, a, a), which should be equivalent from an accuracy, but not a performance, standpoint)

This is more accurate than ever including the $1-t$ explicitly, where I would expect most of the practical precision loss, so fma(t, y, fma(-t, x, x)) in your variables.

I don't have a proof for the $x \leq y \leq z$ inequality, though.

$\endgroup$
  • $\begingroup$ Thank you for the answer. I understand why this formula is appealing and looks better than what I have in the question, but I think in single precision this formula is not monotonic, counterexample: (x,y,t1,t2) = (1.0f0, 3.2732666f0, 2.9802322f-8, 3.352763f-8), z(t1) = 1.0000001f0 > z(t2) = 1.0f0. I changed the question to make it explicit that this is undesirable (for obvious mathematical reasons). $\endgroup$ – Kirill Aug 1 '17 at 19:25

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.