P0811R3 midpoint(), lerp()

[StephanTLavavej]

P0811R3 midpoint(), lerp()

This was partially implemented in VS 2019 16.3, but more work remains to be done. For example, lerp() isn't constexpr yet, and statementreply identified significant issues in #65 (comment) below.

Feature-test macro as of WG21-N4842:
#define __cpp_lib_interpolate 201902L

[utilForever]

Can I implement constexpr lerp() for contribution?

[StephanTLavavej]

Sure! 😸 However, we won’t be able to accept significant functional changes until we get our test/CI infrastructure online in the coming month or two (hopefully). I believe that the lack of constexpr isfinite/isnan or constexpr bit_cast is what’s blocking this part.

[utilForever]

If then, can I also implement constexpr isfinite()/isnan() or constexpr bit_cast? 😃

[sylveon]

bit_cast is going to require a compiler intrinsic

[utilForever]

Charles Milette (@sylveon) Oh, It was my mistake. Thanks. :)

[statementreply]

Disclaimer: I am not an expert on numerical computation.

TL;DR {

The following assertions fail.

constexpr double inf = numeric_limits<double>::infinity();
assert(isnan(lerp(inf, inf, -1.0)));
assert(isnan(lerp(inf, inf, 2.0)));
assert(isnan(lerp(-inf, -inf, -1.0)));
assert(isnan(lerp(-inf, -inf, 2.0)));

}

There is some inconsistency in the existing lerp() implementation when handling special input values (NaN, ±∞, ±0, t = ±0, t = 1), especially lerp(+∞, +∞, t) and lerp(-∞, -∞, t). It is desirable that "well-behaved interpolation" handles special cases consistently.

P0811R3 doesn't fully specify the result of lerp(x, ±∞, ±0) and lerp(±∞, x, 1), nor the sign of zero returned by lerp(±0, b, ±0) and lerp(a, ±0, 1). A design choice needs to be made on how to handle these cases. Two of the possible options are:

• Option A: t = ±0 and t = 1 are special. lerp(x, b, ±0) and lerp(a, x, 1) return x.

  • lerp(x, ±∞, ±0) and lerp(±∞, x, 1) return x.

  • lerp(+0, b, ±0) and lerp(a, +0, 1) return +0.

  • lerp(-0, b, ±0) and lerp(a, -0, 1) return -0.

  • Consequently, lerp(x, qNaN, ±0) and lerp(qNaN, x, 1) also return x. (cf. pow(qNaN, ±0) = 1 and hypot(±∞, qNaN) = hypot(qNaN, ±∞) = +∞)

• Option B: t = ±0 and t = 1 are not quite special.

  • lerp(x, ±∞, ±0) and lerp(±∞, x, 1) return NaN, as they involve 0 × ∞.

  • The sign of lerp(±0, b, ±0) and lerp(a, ±0, 1) is the same as if evaluating lerp(a, b, t) as signbit(a) ? (t-1)*(-a) + t*b : (1-t)*a + t*b.

The following is a standardese-ish description of some expected (I believe) behavior in addition to P0811R3 specification, assuming IEEE 754 floating point:

• For both options
  • lerp(a, b, -∞) returns -∞ when a<b
  • lerp(a, b, -∞) returns +∞ when a>b
  • lerp(a, b, -∞) returns NaN and raises FE_INVALID when a==b
  • lerp(a, b, +∞) returns +∞ when a<b
  • lerp(a, b, +∞) returns -∞ when a>b
  • lerp(a, b, +∞) returns NaN and raises FE_INVALID when a==b
  • lerp(a, -∞, t) returns +∞ when a is finite and t<0
  • lerp(a, -∞, t) returns -∞ when a is finite and t>0
  • lerp(a, +∞, t) returns -∞ when a is finite and t<0
  • lerp(a, +∞, t) returns +∞ when a is finite and t>0
  • lerp(-∞, b, t) returns -∞ when b is finite and t<1
  • lerp(-∞, b, t) returns +∞ when b is finite and t>1
  • lerp(+∞, b, t) returns +∞ when b is finite and t<1
  • lerp(+∞, b, t) returns -∞ when b is finite and t>1
  • lerp(-∞, -∞, t) returns NaN and raises FE_INVALID when t<0
  • lerp(-∞, -∞, t) returns -∞ when t>0 && t<1
  • lerp(-∞, -∞, t) returns NaN and raises FE_INVALID when t>1
  • lerp(-∞, +∞, t) returns -∞ when t<0
  • lerp(-∞, +∞, t) returns NaN and raises FE_INVALID when t>0 && t<1
  • lerp(-∞, +∞, t) returns +∞ when t>1
  • lerp(+∞, -∞, t) returns +∞ when t<0
  • lerp(+∞, -∞, t) returns NaN and raises FE_INVALID when t>0 && t<1
  • lerp(+∞, -∞, t) returns -∞ when t>1
  • lerp(+∞, +∞, t) returns NaN and raises FE_INVALID when t<0
  • lerp(+∞, +∞, t) returns +∞ when t>0 && t<1
  • lerp(+∞, +∞, t) returns NaN and raises FE_INVALID when t>1
• For option A only
  • lerp(a, b, ±0) returns a for any a and b, even when b is a quiet NaN
  • lerp(a, b, +1) returns b for any a and b, even when a is a quiet NaN
• For option B only
  • lerp(a, ±∞, ±0) returns NaN and raises FE_INVALID when a is not a NaN
  • lerp(±∞, b, +1) returns NaN and raises FE_INVALID when b is not a NaN

The full tables of expected (I believe) and implemented special case lerp(a, b, t) behavior are given below. The chosen option should be applied consistently to all input cases. Some cases are listed multiple times for easier reference.

a	b	t	Expected	ImplementedA
-∞	-∞	-∞	qNaN, signal invalid	-∞, no exceptionC
-∞	Finite	-∞	-∞
-∞	+∞	-∞	-∞
Finite	-∞	-∞	+∞
Finite	Finite b<a	-∞	+∞
Finite	b = a	-∞	qNaN, signal invalid
Finite	Finite b>a	-∞	-∞
Finite	+∞	-∞	-∞
+∞	-∞	-∞	+∞
+∞	Finite	-∞	+∞
+∞	+∞	-∞	qNaN, signal invalid	+∞, no exceptionC
-∞	-∞	+∞	qNaN, signal invalid	-∞, no exceptionC
-∞	Finite	+∞	+∞
-∞	+∞	+∞	+∞
Finite	-∞	+∞	-∞
Finite	Finite b<a	+∞	-∞
Finite	b = a	+∞	qNaN, signal invalid
Finite	Finite b>a	+∞	+∞
Finite	+∞	+∞	+∞
+∞	-∞	+∞	-∞
+∞	Finite	+∞	-∞
+∞	+∞	+∞	qNaN, signal invalid	+∞, no exceptionC

a	b	t	Option A	Option B	ImplementedA
±∞	Finite	±0	a	a
Finite	±∞	±0	a	qNaN, signal invalid	Option A
±∞	±∞	±0	a	qNaN, signal invalid	Option A
Finite	±∞	1	b	b
±∞	Finite	1	b	qNaN, signal invalid	Option A
±∞	±∞	1	b	qNaN, signal invalid	Option A

a	b	t	Option A	Option B	ImplementedA
Finite	-∞	-∞	+∞	+∞
Finite	-∞	Finite t<0	+∞	+∞
Finite	-∞	±0	a	qNaN, signal invalid	Option A
Finite	-∞	Finite t>0	-∞	-∞
Finite	-∞	+∞	-∞	-∞
Finite	+∞	-∞	-∞	-∞
Finite	+∞	Finite t<0	-∞	-∞
Finite	+∞	±0	a	qNaN, signal invalid	Option A
Finite	+∞	Finite t>0	+∞	+∞
Finite	+∞	+∞	+∞	+∞
-∞	Finite	-∞	-∞	-∞
-∞	Finite	Finite t<1	-∞	-∞
-∞	Finite	1	b	qNaN, signal invalid	Option A
-∞	Finite	Finite t>1	+∞	+∞
-∞	Finite	+∞	+∞	+∞
+∞	Finite	-∞	+∞	+∞
+∞	Finite	Finite t<1	+∞	+∞
+∞	Finite	1	b	qNaN, signal invalid	Option A
+∞	Finite	Finite t>1	-∞	-∞
+∞	Finite	+∞	-∞	-∞

a	b	t	Option A	Option B	ImplementedA
-∞	-∞	-∞	qNaN, signal invalid	qNaN, signal invalid	-∞, no exceptionC
-∞	-∞	Finite t<0	qNaN, signal invalid	qNaN, signal invalid	-∞, no exceptionC
-∞	-∞	±0	-∞	qNaN, signal invalid	Option A
-∞	-∞	0<t<1	-∞	-∞
-∞	-∞	1	-∞	qNaN, signal invalid	Option A
-∞	-∞	Finite t>1	qNaN, signal invalid	qNaN, signal invalid	-∞, no exceptionC
-∞	-∞	+∞	qNaN, signal invalid	qNaN, signal invalid	-∞, no exceptionC
-∞	+∞	-∞	-∞	-∞
-∞	+∞	Finite t<0	-∞	-∞
-∞	+∞	±0	-∞	qNaN, signal invalid	Option A
-∞	+∞	0<t<1	qNaN, signal invalid	qNaN, signal invalid	qNaN, no exceptionC
-∞	+∞	1	+∞	qNaN, signal invalid	Option A
-∞	+∞	Finite t>1	+∞	+∞
-∞	+∞	+∞	+∞	+∞
+∞	-∞	-∞	+∞	+∞
+∞	-∞	Finite t<0	+∞	+∞
+∞	-∞	±0	+∞	qNaN, signal invalid	Option A
+∞	-∞	0<t<1	qNaN, signal invalid	qNaN, signal invalid	qNaN, no exceptionC
+∞	-∞	1	-∞	qNaN, signal invalid	Option A
+∞	-∞	Finite t>1	-∞	-∞
+∞	-∞	+∞	-∞	-∞
+∞	+∞	-∞	qNaN, signal invalid	qNaN, signal invalid	+∞, no exceptionC
+∞	+∞	Finite t<0	qNaN, signal invalid	qNaN, signal invalid	+∞, no exceptionC
+∞	+∞	±0	+∞	qNaN, signal invalid	Option A
+∞	+∞	0<t<1	+∞	+∞
+∞	+∞	1	+∞	qNaN, signal invalid	Option A
+∞	+∞	Finite t>1	qNaN, signal invalid	qNaN, signal invalid	+∞, no exceptionC
+∞	+∞	+∞	qNaN, signal invalid	qNaN, signal invalid	+∞, no exceptionC

a	b	t	Option A	Option B	ImplementedA
-0	-0	-∞	qNaN, signal invalid	qNaN, signal invalid
-0	-0	Finite t<0	±0B	±0B
-0	-0	-0	-0	±0B	Option B
-0	-0	+0	-0	-0
-0	-0	0<t<1	-0	-0
-0	-0	1	-0	±0B	Option A
-0	-0	Finite t>1	±0B	±0B
-0	-0	+∞	qNaN, signal invalid	qNaN, signal invalid
-0	+0	-∞	qNaN, signal invalid	qNaN, signal invalid
-0	+0	Finite t<0	-0	-0
-0	+0	-0	-0	-0
-0	+0	+0	-0	±0B	Option B
-0	+0	0<t<1	±0B	±0B
-0	+0	1	+0	±0B	Option A
-0	+0	Finite t>1	+0	+0
-0	+0	+∞	qNaN, signal invalid	qNaN, signal invalid
+0	-0	-∞	qNaN, signal invalid	qNaN, signal invalid
+0	-0	Finite t<0	+0	+0
+0	-0	-0	+0	+0
+0	-0	+0	+0	±0B	Option B
+0	-0	0<t<1	±0B	±0B
+0	-0	1	-0	±0B	Option B
+0	-0	Finite t>1	-0	-0
+0	-0	+∞	qNaN, signal invalid	qNaN, signal invalid
+0	+0	-∞	qNaN, signal invalid	qNaN, signal invalid
+0	+0	Finite t<0	±0B	±0B
+0	+0	-0	+0	±0B	Option B
+0	+0	+0	+0	+0
+0	+0	0<t<1	+0	+0
+0	+0	1	+0	±0B	Option B
+0	+0	Finite t>1	±0B	±0B
+0	+0	+∞	qNaN, signal invalid	qNaN, signal invalid

a	b	t	Option A	Option B	ImplementedA
-0	Finite b<0	-0	-0	±0B	Option B
-0	-0	-0	-0	±0B	Option B
-0	+0	-0	-0	-0
-0	Finite b>0	-0	-0	-0
+0	Finite b<0	-0	+0	+0
+0	-0	-0	+0	+0
+0	+0	-0	+0	±0B	Option B
+0	Finite b>0	-0	+0	±0B	Option B
-0	Finite b<0	+0	-0	-0
-0	-0	+0	-0	-0
-0	+0	+0	-0	±0B	Option B
-0	Finite b>0	+0	-0	±0B	Option B
+0	Finite b<0	+0	+0	±0B	Option B
+0	-0	+0	+0	±0B	Option B
+0	+0	+0	+0	+0
+0	Finite b>0	+0	+0	+0
Finite a<0	-0	1	-0	±0B	Option A
-0	-0	1	-0	±0B	Option A
+0	-0	1	-0	±0B	Option B
Finite a>0	-0	1	-0	±0B	Option B
Finite a<0	+0	1	+0	±0B	Option A
-0	+0	1	+0	±0B	Option A
+0	+0	1	+0	±0B	Option B
Finite a>0	+0	1	+0	±0B	Option B

a	b	t	Option A	Option B	ImplementedA
Finite	qNaN	±0	a	qNaN	Option B
±∞	qNaN	±0	a	qNaN	Option B
qNaN	Finite	±0	qNaN	qNaN
qNaN	±∞	±0	qNaN	qNaN
qNaN	qNaN	±0	qNaN	qNaN
Finite	qNaN	1	qNaN	qNaN
±∞	qNaN	1	qNaN	qNaN
qNaN	Finite	1	b	qNaN	Option B
qNaN	±∞	1	b	qNaN	Option B
qNaN	qNaN	1	qNaN	qNaN
Finite	qNaN	Finite t≠0 and t≠1	qNaN	qNaN
±∞	qNaN	Finite t≠0 and t≠1	qNaN	qNaN
qNaN	Finite	Finite t≠0 and t≠1	qNaN	qNaN
qNaN	±∞	Finite t≠0 and t≠1	qNaN	qNaN
qNaN	qNaN	Finite t≠0 and t≠1	qNaN	qNaN
Finite	qNaN	±∞	qNaN	qNaN
±∞	qNaN	±∞	qNaN	qNaN
qNaN	Finite	±∞	qNaN	qNaN
qNaN	±∞	±∞	qNaN	qNaN
qNaN	qNaN	±∞	qNaN	qNaN
Finite	qNaN	qNaN	qNaN	qNaN
±∞	qNaN	qNaN	qNaN	qNaN
qNaN	Finite	qNaN	qNaN	qNaN
qNaN	±∞	qNaN	qNaN	qNaN
qNaN	qNaN	qNaN	qNaN	qNaN

a	b	t	Expected	ImplementedA
sNaN	Any	Any	qNaN, signal invalid	sNaN, no exceptionC
Any	sNaN	Any	qNaN, signal invalid	qNaN or sNaN, no exceptionC
Any	Any	sNaN	qNaN, signal invalid	qNaN or sNaN, no exceptionC

^A Omitted when implemented behavior and expected behavior of both options are all the same.

^B -0 when rounding direction is downward, +0 otherwise.

^C IEEE 754 exception, not C++ exception.

[StephanTLavavej]

Thanks for the careful review; we'll need to resolve that design choice and write corresponding tests before completing this feature.

[BillyONeal]

The following assertions fail.
assert(isnan(lerp(inf, inf, -1.0)));
assert(isnan(lerp(inf, inf, 2.0)));
assert(isnan(lerp(-inf, -inf, -1.0)));
assert(isnan(lerp(-inf, -inf, 2.0)));

If I understand correctly, a large part of the point of lerp is to make those cases not result in NaNs; in fact much of the implementation complexity comes from wanting infinities to pass through here. The linear interpolation between two of the same value is always the same value, so inf / -inf seems to be the correct answer in these places, regardless of the value of T (unless T is a NaN).

[BillyONeal]

(That is to say, the rows that you've highlighted in bold above have the implementation doing what I believe to be the correct behavior; your option A and option B refer to the value of t, but the value of t doesn't matter when a and b have the same value)

[statementreply]

The linear interpolation between two of the same value is always the same value, so inf / -inf seems to be the correct answer in these places, regardless of the value of T (unless T is a NaN).

lerp(x, x, t) is equal to x for all finite x and t, but it doesn't necessarily extend to infinities. When t < 0 or t > 1, lerp(a, b, t) becomes a linear extrapolation of a and b. The limit of such extrapolation of two infinities of the same sign may not exist, so I believe NaN is correct. The same reasoning applies to lerp(+∞, -∞, t) returning NaN when 0 < t < 1.

Example: x - x is 0 for all finite x, but NaN for x = ±∞.

Counter-example: atan2(x, x) is +π/4 for all positive x, including x = +∞. There's a different consideration for this: atan2(y, x) is expected to return the angle in the correct quadrant.

[BillyONeal]

The limit of such extrapolation of two infinities of the same sign may not exist

I don't see how it does not exist. They are the same value, any extrapolated value should be equal to the same value.

The same reasoning applies to lerp(+∞, -∞, t) returning NaN when 0 < t < 1.

That case is different because there's an actual extrapolation going on. If you plot (0, +∞) and (1, -∞), there is no defined value in the range (0, 1), hence the result being NaN. But if you plot (0, +∞), (1, +∞), the resulting plot is a horizontal line at +∞.

Example: x - x is 0 for all finite x, but NaN for x = ±∞.

That's not the behavior I observe:

C:\Users\bion\Desktop>type test.cpp
#include <assert.h>
#include <limits>
#include <math.h>
#include <stdio.h>

int main() {
    const auto inf = std::numeric_limits<double>::infinity();
    assert(!isnan(inf + inf));
    assert(isnan(inf - inf));
    assert(isnan(inf + (-inf)));
    assert(!isnan(inf - (-inf)));
    puts("pass");
}

C:\Users\bion\Desktop>cl /EHsc /W4 /WX .\test.cpp
Microsoft (R) C/C++ Optimizing Compiler Version 19.24.28314 for x86
Copyright (C) Microsoft Corporation.  All rights reserved.

test.cpp
Microsoft (R) Incremental Linker Version 14.24.28314.0
Copyright (C) Microsoft Corporation.  All rights reserved.

/out:test.exe
test.obj

C:\Users\bion\Desktop>.\test.exe
pass

[BillyONeal]

(I guess your "x - x is 0 is NaN for x = ±∞" statement holds, but my point is that infinities of the same sign produce an infinity of the same sign, you only get NaN when you try to mix infinities of different signs, which is I believe what the current implementation does based on

STL/stl/inc/cmath

Line 1158 in 712b797

// same sign, so T doesn't matter

(though we might not set the right floating point exceptions))

[statementreply]

I don't see how it does not exist. They are the same value, any extrapolated value should be equal to the same value.

But if you plot (0, +∞), (1, +∞), the resulting plot is a horizontal line at +∞.

I guess you are referring to the case of e.g. lim_x→+∞lerp(x, x, 2) = +∞. However, lim_{x→+∞, y→+∞}lerp(x, y, 2) may be either fininte, infinite or divergent, and lim_{x→+∞, y→+∞}{plot (0, x), (1, y)} may be diagonal.

The same for (+∞) - (+∞) → NaN: lim_x→+∞(x - x) = 0, but lim_{x→+∞, y→+∞}(x - y) can be different.

[BillyONeal]

However, limx→+∞, y→+∞lerp(x, y, 2) may be either fininte, infinite or divergent

But we aren't talking about limits here. We have concrete values. If floating point math operated in this "limits" domain then ∞ + ∞ would be NaN, not ∞.

The same for (+∞) - (+∞) → NaN

I don't think that case is related at all, because those infinities are of different signs.