Add lerp to the standard math library

[FunnyGlaggle]

Add lerp to the standard math library

[gaymeowing]

Could we get this merged? Doesn't seem like theres any objections to it.

[vegorov-rbx]

Suggested wording change has never been applied. I would say there's no rush in getting it in.
Team is also focused on delivering other projects, so no need to rush things.

Additionally, the implementation is hard to follow with respect to the guarantees presented and whether or not we need them.
It doesn't match any of the popular implementations mentioned in the Rust issue.
It is also not as precise as the one from C++ (which I don't think is valuable to have).
We currently don't have time to do an analysis on numerical claims of that implementation, maybe the author can give a link to a paper that does that? I think that while special-casing at 3 points makes it exact, but might actually remove monotonicity.

Feels like that the developers will be better off implementing their usual custom implementation of a simple x + (y - x) * alpha.
Current implementation is complex enough with 3 condition statements to not be implemented in native codegen, so a custom function will win on performance.
My suggestion would be that developers who need precise numerical properties can have their own implementation with them, but that the one in Luau can have some properties missing, but be equal to the general expectation of what lerp does.

Here's the objection you wanted :)

[AxisAngles]

lerp(a0, a1, t):

• t can be much higher precision when near 0
  • so when t is near 0, we don't ever want to have a function of the form A + (1 - t)*B
• we want exactness, lerp(a0, a1, 0) == a0 and lerp(a0, a1, 1) == a1
  • when t is closer to 0, we want higher accuracy around a0
  • when t is closer to 1, we want higher accuracy around a1
• we want functions of the form A + t*B, as this guarantees monotonic behavior

Some things to consider for building intuition:

-- (100000000 samples)

-- given x, y = sortabs(random()*randomsign(), random()*randomsign()):
--   x - y + y ~= x: ~49.2%
--   y - x + x ~= y: ~20.6%

-- given x, y = sortabs(random(), random()):
--   x - y + y ~= x: ~30.1%
--   y - x + x ~= y: ~3.98%

-- given x, y = sortabs(random()*1.5, random()*1.5):
--   x - y + y ~= x: ~25.4%
--   y - x + x ~= y: ~1.63%

-- given floorlog2(x) == floorlog2(y) >= floorlog2(y - x)
--   x - y + y ~= x: 0%
--   y - x + x ~= y: 0%

Taking all this into account, we build our first example:

• Obvious solution is to split when t < 1/2
  • at the interchange point, t = 1/2, we want to guarantee monotonic behavior
  • so we bound the result with respect to (a0 + a1)/2

--[[
guarantees exactness because
	a0 + 0*(a1 - a0) = a0
	a1 + (1 - 1)*(a1 - a0) = a1 + 0*(a1 - a0) = a1
guarantees consistency because
	a + (a - a)*t = a + 0*t = a
guarantees monotonicity because
	we bound the results such that
	lerp(a0, a1, 1/2 - 2^-54) <= lerp(a0, a1, 1/2)
]]
local function lerp1(a0, a1, t)
	local m = (a0 + a1)/2
	if t < 1/2 then
		local a = a0 + (a1 - a0)*t
		if a0 < a1 then
			return math.min(m, a)
		else
			return math.max(m, a)
		end
	else
		local a = a1 + (a1 - a0)*(t - 1)
		if a0 < a1 then
			return math.max(m, a)
		else
			return math.min(m, a)
		end
	end
end

Then we build a simpler case, to make sure what we are doing is necessary:

• The relevant question is
  • can a0 + (1/2 - 2^-54)*(a1 - a0) > a1 - 1/2*(a1 - a0) ever be true?
  • From a lot of random sampling, it appears to be impossible, or at least very rare

-- this guarantees exactness and consistency
-- this appears to guarantee monotonicity intrinsically
local function lerp2(a0, a1, t)
	if t < 1/2 then
		return a0 + (a1 - a0)*t
	else
		return a1 + (a1 - a0)*(t - 1)
	end
end

[dphblox]

so from that, seems like at most all we'd need is to conditionally swap a <-> b and invert t to reach most/all the guarantees we wanted? aka no reason to go more complex than lerp2?

[dphblox]

blending in @/vegorov-rbx's comments, would it even be worth special casing anything? would the only issue with using a0 + (a1 - a0) * t be that t has less precision around 1?

[AxisAngles]

t has less precision around 1, but there's no way to extract more precision from t in the first place, what they give is what we get to work with.

Edit:
using only a0 + (a1 - a0)*t is problematic because often times, when t = 1, it does not return a1

[AxisAngles]

From a user perspective, if I want more accuracy around a1 than a0, I would need to find some accurate way of generating s = 1 - t directly, and then use math.lerp(a1, a0, s)

[dphblox]

Oh right this is about the precision of (a1 - a0), right? We want to take advantage of the higher density of floating point numbers near 0 so we would rather "multiply down". e.g.: x * 0.01 will be more precise than x - (x * 0.99).

All makes sense - lerp2 seems like the nicest "correct" solution. I suppose it just depends on whether we would want to use a conditional or whether that's deemed too heavy - will let the Luau folk speak to that.

[AxisAngles]

In my experience as a user, the most important guarantees lerp can have is for
lerp(a0, a1, 0) == a0 and lerp(a0, a1, 1) == a1

If conditionals are not allowed, then the optimal solution seems to be
(1 - t)*a0 + t*a1
but this is neither monotonic nor consistent.

[FunnyGlaggle]

Oh this is actually really clever, so the only "inaccuracy" we get is near 0.5

[FunnyGlaggle]

Thank you for this new insight, I will try and push this forward :)

[zeux]

Some thoughts:

1. Do we need math.lerp? When this was raised at some point a couple years ago I was uncertain: math.lerp is useful wrt ergonomics, but a few uncertain properties around its implementation make it difficult to pick "the best one". Notably, Rust still hasn't added one. C++ has added a function that has guarantees that you'd want but the implementation is the one you'd never want if performance is even a little valuable.

However, the decision has been made to merge math.map; as such, I think math.lerp should be added: lerp is useful more often than map, math was added with no regard for any numerical properties, and while map can be converted into lerp by using 0,1 limits, that use is not ergonomic and not as performant (short of codegen cases where in theory map can be optimized optimally via partial constant folding).

So I support merging some RFC here in principle.

2. Proposed implementation is redundant: there is no need to check x==0 or x==y, this doesn't give any more guarantees. So I think this RFC should be amended to simplify the proposed lowering.

3. There's five properties I think we could consider:

a. exactness: lerp(0) = a, lerp(1) = b
b. monotonicity: lerp(t1) >= lerp(t2) for t1 >= t2
c. determinacy: lerp(a, b, t) is NaN only if t is infinite or a/b/t are NaN
d. boundedness: lerp(t) in [a, b]
e. consistency: lerp(x, x, t) == x for all t

C++ implementation satisfies all 5. However, the result is not useful wrt performance. In my opinion, determinacy is not particularly interesting as violating it requires fp overflow which is too strong of a constraint and not practically useful; exactness and consistency are important; boundedness is sometimes important; monotonicity is nice to have, as I do not know of practical applications where monotonicity is required.

This makes me wonder: why not choose if t == 1 then b else a + (b - a) * t as a suggested implementation? It is exact and consistent; it is bounded for reasonable values of a/b (is it always bounded?) and monotonic (for t in [0, 1]). If b-a overflows you get a bad lerp but that seems too strong of a constraint to satisfy when performance is important.

C++ implementation in https://www.open-std.org/jtc1/sc22/wg21/docs/papers/2019/p0811r3.html is basically that for a/b of the same sign, with an extra case for t=1+eps; I don't fully understand the value of monotonicity outside of the 0..1 interval, and as I mentioned above I don't know why monotonicity is important inside the interval, but certainly that property is not worth the extra branches.

I've empirically validated that the implementation above works well on 32-bit floats by bruteforcing one of a/b/t and fixing a few other values (for example, I validated boundedness by fixing a, fixing t at prevfloat(1), and bruteforcing all b >= a). This SO question https://math.stackexchange.com/questions/907327/accurate-floating-point-linear-interpolation has some discussion in the comments and comes to the same conclusion afaict.

I like the t<0.5 form in principle (note that it's mentioned in the comments of the issue linked above), but it has a steep implementation cost: you need to do two conditional selects on t<0.5 condition (switching a/b and t/t-1), plus compute t-1, and load two constants (0.5 & 1) into registers. A branch around 0.5 will behave well on microbenchmarks but cause branch mispredictions in some cases. Compared to this, t==1 is just one conditional select & load of just one constant, so it's cheaper.

4. Because the exact mechanics of lerp are tricky, I think it's probably fine to specify the above as a suggested implementation, and note that we want to guarantee exactness (and maybe consistency?), but not have explicit language around monotonicity and other properties.

Interestingly, the formula a + (b - a) * t without corrections is not guaranteed to be bounded, however it's often bounded for interesting intervals that you choose to interpolate where boundedness is important (notably, intervals 0..X, -X..X, and in general +X..+Y seem to be bounded in practice?). But it's not exact, and it feels like it's fine to leave boundedness unspecified for now.

[zeux]

Oh and cost & properties of a few implementations; lerp 0 is a + (b - a) * t, lerp 1 is a * (1 - t) + b * t, lerp 2 is the implementation I suggest above, lerp 3 is the implementation with t<0.5 branch, lerp 4 is lerp 1 that enforces consistency via explicit a == b branch (from Rust discussion thread). Cost is in instruction count on A64/X64 respectively; X64 lowering is as generated by clang on Apple platform, not guaranteed to be optimal...

// lerp 0 (1mul): 3/3 insn; exact 0, monotonic 1, bounded 0, consistent 1
// lerp 1 (2mul): 5/5 insn; exact 1, monotonic 0, bounded 0, consistent 0
// lerp 2 (1mul t==1): 6/9 insn; exact 1, monotonic 1, bounded 1, consistent 1
// lerp 3 (1mul t<.5): 9/16 insn; exact 1, monotonic 1, bounded 1, consistent 1
// lerp 4 (2mul a==b): 7/9 insn; exact 1, monotonic 0, bounded 0, consistent 1

[zeux]

I've submitted #86 which intends to replace this RFC while providing the same functionality via math.lerp with the adjusted implementation following comments above and more fuzz testing.

[vegorov-rbx]

Closing this as #86 has a better implementation if we were to include it.