How can I interpolate $z_t = x(1-t)+y t$ with single-precision floats so that it satisfies $x\leq z\leq y$, $z_0=x$, $z_1=y$?

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?

• Closely related question and possibly useful answer on Math Stack Exchange – njuffa Aug 7 '17 at 18:41
• @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. – Kirill Aug 9 '17 at 20:50
• 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. – njuffa Aug 10 '17 at 3:54

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.
• 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). – Kirill Aug 1 '17 at 19:25