Stl: 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 https://github.com/microsoft/STL/issues/65#issuecomment-562868441 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]

@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 | Implemented^A
:---: | :---: | :---: | :------: | :---------------------:
-∞ | -∞ | -∞ | qNaN, signal invalid | -∞, no exception^C
-∞ | Finite | -∞ | -∞
-∞ | +∞ | -∞ | -∞
Finite | -∞ | -∞ | +∞
Finite | Finite b<a | -∞ | +∞
Finite | b = a | -∞ | qNaN, signal invalid
Finite | Finite b>a | -∞ | -∞
Finite | +∞ | -∞ | -∞
+∞ | -∞ | -∞ | +∞
+∞ | Finite | -∞ | +∞
+∞ | +∞ | -∞ | qNaN, signal invalid | +∞, no exception^C
-∞ | -∞ | +∞ | qNaN, signal invalid | -∞, no exception^C
-∞ | Finite | +∞ | +∞
-∞ | +∞ | +∞ | +∞
Finite | -∞ | +∞ | -∞
Finite | Finite b<a | +∞ | -∞
Finite | b = a | +∞ | qNaN, signal invalid
Finite | Finite b>a | +∞ | +∞
Finite | +∞ | +∞ | +∞
+∞ | -∞ | +∞ | -∞
+∞ | Finite | +∞ | -∞
+∞ | +∞ | +∞ | qNaN, signal invalid | +∞, no exception^C

a | b | t | Option A | Option B | Implemented^A
:---: | :---: | :---: | :------: | :------: | :---------------------:
±∞ | 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 | Implemented^A
:---: | :---: | :---: | :------: | :------: | :---------------------:
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 | Implemented^A
:---: | :---: | :---: | :------: | :------: | :---------------------:
-∞ | -∞ | -∞ | qNaN, signal invalid | qNaN, signal invalid | -∞, no exception^C
-∞ | -∞ | Finite t<0 | qNaN, signal invalid | qNaN, signal invalid | -∞, no exception^C
-∞ | -∞ | ±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 exception^C
-∞ | -∞ | +∞ | qNaN, signal invalid | qNaN, signal invalid | -∞, no exception^C
-∞ | +∞ | -∞ | -∞ | -∞
-∞ | +∞ | Finite t<0 | -∞ | -∞
-∞ | +∞ | ±0 | -∞ | qNaN, signal invalid | Option A
-∞ | +∞ | 0<t<1 | qNaN, signal invalid | qNaN, signal invalid | qNaN, no exception^C
-∞ | +∞ | 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 exception^C
+∞ | -∞ | 1 | -∞ | qNaN, signal invalid | Option A
+∞ | -∞ | Finite t>1 | -∞ | -∞
+∞ | -∞ | +∞ | -∞ | -∞
+∞ | +∞ | -∞ | qNaN, signal invalid | qNaN, signal invalid | +∞, no exception^C
+∞ | +∞ | Finite t<0 | qNaN, signal invalid | qNaN, signal invalid | +∞, no exception^C
+∞ | +∞ | ±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 exception^C
+∞ | +∞ | +∞ | qNaN, signal invalid | qNaN, signal invalid | +∞, no exception^C

a | b | t | Option A | Option B | Implemented^A
:---: | :---: | :---: | :------: | :------: | :---------------------:
-0 | -0 | -∞ | qNaN, signal invalid | qNaN, signal invalid
-0 | -0 | Finite t<0 | ±0^B | ±0^B
-0 | -0 | -0 | -0 | ±0^B | Option B
-0 | -0 | +0 | -0 | -0
-0 | -0 | 0<t<1 | -0 | -0
-0 | -0 | 1 | -0 | ±0^B | Option A
-0 | -0 | Finite t>1 | ±0^B | ±0^B
-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 | ±0^B | Option B
-0 | +0 | 0<t<1 | ±0^B | ±0^B
-0 | +0 | 1 | +0 | ±0^B | 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 | ±0^B | Option B
+0 | -0 | 0<t<1 | ±0^B | ±0^B
+0 | -0 | 1 | -0 | ±0^B | 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 | ±0^B | ±0^B
+0 | +0 | -0 | +0 | ±0^B | Option B
+0 | +0 | +0 | +0 | +0
+0 | +0 | 0<t<1 | +0 | +0
+0 | +0 | 1 | +0 | ±0^B | Option B
+0 | +0 | Finite t>1 | ±0^B | ±0^B
+0 | +0 | +∞ | qNaN, signal invalid | qNaN, signal invalid

a | b | t | Option A | Option B | Implemented^A
:---: | :---: | :---: | :------: | :------: | :---------------------:
-0 | Finite b<0 | -0 | -0 | ±0^B | Option B
-0 | -0 | -0 | -0 | ±0^B | 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 | ±0^B | Option B
+0 | Finite b>0 | -0 | +0 | ±0^B | Option B
-0 | Finite b<0 | +0 | -0 | -0
-0 | -0 | +0 | -0 | -0
-0 | +0 | +0 | -0 | ±0^B | Option B
-0 | Finite b>0 | +0 | -0 | ±0^B | Option B
+0 | Finite b<0 | +0 | +0 | ±0^B | Option B
+0 | -0 | +0 | +0 | ±0^B | Option B
+0 | +0 | +0 | +0 | +0
+0 | Finite b>0 | +0 | +0 | +0
Finite a<0 | -0 | 1 | -0 | ±0^B | Option A
-0 | -0 | 1 | -0 | ±0^B | Option A
+0 | -0 | 1 | -0 | ±0^B | Option B
Finite a>0 | -0 | 1 | -0 | ±0^B | Option B
Finite a<0 | +0 | 1 | +0 | ±0^B | Option A
-0 | +0 | 1 | +0 | ±0^B | Option A
+0 | +0 | 1 | +0 | ±0^B | Option B
Finite a>0 | +0 | 1 | +0 | ±0^B | Option B

a | b | t | Option A | Option B | Implemented^A
:---: | :---: | :---: | :------: | :------: | :---------------------:
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 | Implemented^A
:---: | :---: | :---: | :------: | :---------------------:
sNaN | Any | Any | qNaN, signal invalid | sNaN, no exception^C
Any | sNaN | Any | qNaN, signal invalid | qNaN or sNaN, no exception^C
Any | Any | sNaN | qNaN, signal invalid | qNaN or sNaN, no exception^C

^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 https://github.com/microsoft/STL/blob/712b7971bdddd72cdaeecccfb82d5f8fd62e23ac/stl/inc/cmath#L1158 (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.

[BillyONeal]

I pinged the L[E]WG mailing lists for the designers of lerp to comment on this issue.

[hfinkel]

I guess you are referring to the case of e.g. limx→+∞lerp(x, x, 2) = +∞. However, limx→+∞, y→+∞lerp(x, y, 2) may be either fininte, infinite or divergent, ...

@statementreply, if I have two functions, both diverging to positive infinity, under what conditions might the interpolation between the two not diverge to positive infinity?

[statementreply]

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 ∞.

This is, arguably, how infinities work in IEEE floating point. (And lim_{x→+∞, y→+∞}(x + y) = +∞, hence (+∞) + (+∞) = +∞)

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

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

Extrapolation (extend the difference beyond the boundary), unlike interpolation (blend between the boundaries), involves the difference of the two values; so I believe the subtraction case is relevant.

Try thinking about lerp(a, b, t) as b + (t-1) * (b-a) when t > 1, and as a + (-t) * (a - b) when t < 0.

[statementreply]

@statementreply, if I have two functions, both diverging to positive infinity, under what conditions might the interpolation between the two not diverge to positive infinity?

It is about the extrapolation case (t < 0 or t > 1).
a(x) = 3x, b(x) = x, lerp(a, b, 2) = -x diverges to negative infinity.

[opensdh]

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 best rigorous mathematical basis for IEC 60559 is the affinely extended real number system, which is in fact based on limits. In that system, ∞ + ∞ = ∞ because lim_{x→+∞, y→+∞}(x + y) = ∞ (for all paths of approach, which is the definition of such a compound limit). This is why ∞ / ∞ → NaN: while obviously lim_t→+∞(t/t) = 1, for arbitrary _x_(_t_) and _y_(_t_) the limit can be 1, 0, ∞, or anything else. (NaN is of course not part of the extended real line; it's what 60559 uses to tell you that the result was undefined, and has its infectious propagation to draw attention to the failure.)

There are exceptions, however: std::pow(std::numeric_limits<double>::infinity(),0)==1 (with 60559 support) despite the fact that lim_t→+∞((2^t)^1/t) = 2, because of things like lim_t→+∞((_a_ _t_)^b/t) = 1 for all (finite) a and b. So it's reasonable to ask what behavior we want here.

P0811R3 deliberately leaves some less-useful cases unspecified to allow efficient implementation of the more-useful cases. Of course, if implementors establish that some additional guarantee does not carry a significant cost in any high-quality implementation, it's obviously a net benefit to provide it. (We might need to consider conflicting available sets of such guarantees, though.)

All that said, the basis for the function is either of the expressions a+t*(b-a) or (1-t)*a+t*b. In regard to the particular question posed to LEWG about lerp(inf,inf,t), the first expression is NaN for all _t_ and the second is NaN for _t_ outside (0,1). From a user's (not the author's) perspective, I'm fine with saying that that's sufficient evidence to make it be NaN there; I don't have a strong opinion about _t_ in {0,1} (for infinite _a_ and/or _b_!) because the limit behavior is again path-dependent (unless you fix _t_ first).

[hfinkel]

@statementreply, if I have two functions, both diverging to positive infinity, under what conditions might the interpolation between the two not diverge to positive infinity?

It is about the extrapolation case (t < 0 or t > 1).
a(x) = 3x, b(x) = x, lerp(a, b, 2) = -x diverges to negative infinity.

Thanks, I agree. NaN seems appropriate for the extrapolation cases.

[BillyONeal]

Try thinking about lerp(a, b, t) as b + (t-1) * (b-a) when t > 1, and as a + (-t) * (a - b) when t < 0.
...
All that said, the basis for the function is either of the expressions a+t(b-a) or (1-t)a+t*b.

My understanding was that neither of those expressions are intrinsically related to what lerp attempts to calculate; they were only a means to the end of "given an unknown function y = mx+b that crosses (0, b), (1, a), what is the value of y given x == t?". The only solution for a==b==infinity is the function f(x) = infinity, which does not depend on the value of x (or t). But for opposite signed infinities, there are no values that can be assigned to produce a true expression, hence, NaN.

That said, I'm not going to argue with the designers of the feature. I find returning NaN for those cases to be extremely surprising, but I'm outvoted since it isn't my feature :)

[BillyONeal]

P0811R3 deliberately leaves some less-useful cases unspecified to allow efficient implementation of the more-useful cases.

We need an fpclassify on each of the inputs to correctly handle NaN inputs so I don't think this restriction helps much...

[statementreply]

@opensdh thanks for the clarification.

There are exceptions, however: std::pow(std::numeric_limits<double>::infinity(),0)==1 (with 60559 support) despite the fact that limt→+∞((2t)1/t) = 2, because of things like limt→+∞((_a_ _t_)b/t) = 1 for all (finite) a and b. So it's reasonable to ask what behavior we want here.

I don't have a strong opinion about _t_ in {0,1} (for infinite _a_ and/or _b_!) because the limit behavior is again path-dependent (unless you fix _t_ first).

From IEEE 754 (draft):

NOTE—This standard defines several operations to raise x to a given power:
pown(x, n) accepts integral powers n only, for any x
pow(x, y) behaves the same as pown(x, n) when y contains a value which is equal to an integral n
powr(x, y) is defined by considering exp(y × log(x)), and thus its domain excludes negative x

pow is an interesting case here: y ∈ Z has special power; for y ∈ Z (and infinite x), y is fixed first when taking the limit. powr is the version without special power. Similarly, for lerp and t ∈ {0,1}, either the pow approach (i.e. my option A) or the powr approch (i.e. my option B) could be taken.

[BillyONeal]

FWIW between the two options, option B seems the most consistent with other operations (where NaN is infectious) so I'd probably choose that, but as we've already established I apparently don't know what people actually want this feature to do :)

[opensdh]

My understanding was that neither of those expressions are intrinsically related to what lerp attempts to calculate; they were only a means to the end of "given an unknown function y = mx+b that crosses (0, b), (1, a), what is the value of y given x == t?".

You're right about the intent in general, but for most cases where _a_ and/or _b_ is infinite there is no such function (of that simple form) in the affine reals. It thus seems more natural to appeal to how 60559 treats other cases (including the basic arithmetic in the expressions I mentioned) to determine how to handle such exceptional cases (if we bother to handle them in any consistent way).

We need an fpclassify on each of the inputs to correctly handle NaN _inputs_ so I don't think this restriction helps much...

If any of the inputs is NaN, none of the P0811R3 conditions apply (note that at least one of CMP(t2,t1) or CMP(b,a) is 0 in that case), so I'm not sure how NaN inputs can cause any trouble. (But if the function does in fact have effectively the same speed regardless of how many special cases we support, I'm all for nailing down every case like C does for pow in N2310:F.10.4.4.)

[BillyONeal]

I'm not sure how NaN inputs can cause any trouble

We need to preserve the usual rules for NaNs, that being that if any input is a NaN, the result is also a NaN. We do need the real fpclassify family rather than something like x == x test because some of the implementations we target will remove such a test entirely in some modes (e.g. msvc++ under /fp:fast, or clang -fassociative-math). Note that support for such modes is a place where our implementation differs from some other common implementations -- I believe gcc's -ffast-math also breaks the isnan function, but we consider that the escape hatch that lets programs that would otherwise be happy with these more lax modes more clearly specify when they actually need floating point model guarantees. (@manbearian and his team or other floating point experts can probably specify what guarantees we intend to provide more clearly; my level of understanding is effectively 'don't say x==x')

I'm all for nailing down every case like C does for pow in N2310:F.10.4.4

I remain frustrated that the standard refuses to specify anything about how floating point stuff should behave out of fear of depending on IEEE754 / IEC60559. We should specify what '754 machines do, then say something like "non '754 machines do something as close to this as reasonable" like we did for filesystem. Then we'll at least have the language we need to specify how things should work on most implementations.

[opensdh]

We need to preserve the usual rules for NaNs, that being that if any input is a NaN, the result is also a NaN.

Is that behavior actually required anywhere? I don't see any blanket wording about it, but of course blanket wording is good at hiding. I do see that my proof-of-principle implementation doesn't have that property, but only because of its trick with std::min and std::max preferring their first argument.

Note that support for such modes is a place where our implementation differs from some other common implementations -- I believe gcc's -ffast-math also breaks the isnan _function_, but we consider that the escape hatch that lets programs that would otherwise be happy with these more lax modes more clearly specify when they actually need floating point model guarantees.

It seems that you then want an intrinsic __builtin_isnan_no_really that can be whatever fast check without being elided.

I remain frustrated that the standard refuses to specify anything about how floating point stuff should behave out of fear of depending on IEEE754 / IEC60559. We should specify what '754 machines do, then say something like "non '754 machines do something as close to this as reasonable" like we did for filesystem. Then we'll at least have the language we need to specify how things should work on most implementations.

I certainly agree: that approach seems like an improvement on C's Annex F (which is already much better than the nothing that we have). Just as soon as I finish rewriting the rest of the core language, I'll get right on that.

[jorgbrown]

Wanted this here for the record and easy linkage: lerp spec

Re: inf - inf being NaN, that's true... however inf + NaN is inf.

What about lerp(1, -1e18, 1e-18) ?
What about the similar lerp(1, -pow(2.0, 60), pow(2.0, -60)) ?
(I've been playing with these on godbolt because they have catastrophic cancellation)

[vinniefalco]

From @pdimov

RE: lerp(infinity, infinity, t)

Clear NaN for t < 0 or t > 1, clear inf for t between 0 and 1, non-obvious NaN at 0 and 1.

A plot of (0, infinity) -> (1, infinity) is not a horizontal line, the line can have any inclination. "Infinity" is a shorthand for "when t approaches x, the functions grows without a bound", but one "infinity" can grow as 2/|t-x|, another can grow as 5/|t-x| or 1/(t-x)^2.

Or in terms of floating point "math": lerp(a, b, t) is not NaN when at least one of the three methods of computing it isn't NaN.

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

When a=b=inf, a-b and b-a are NaN, so the first two are NaN. The third one for t between 0 and 1 is inft + inf(1-t) = inf + inf = inf. For 0 or 1, we get inf * 0, which is NaN. For t < 0 or t > 1, we get inf - inf, which is NaN.

non-obvious NaN at 0 and 1.

V. pedantically speaking, since 0 has a sign, lerp(inf, inf, +0) is inf, and only lerp(inf, inf, -0) is NaN. But, since there's no way to determine the sign of 1-t when t=1, it's probably better to disregard this for symmetry reasons.

[BillyONeal]

It seems that you then want an intrinsic __builtin_isnan_no_really that can be whatever fast check without being elided.

That's std::isnan for us.

[opensdh]

Re: inf - inf being NaN, that's true... however inf + NaN is inf.

It is? Not on the implementations I have…

What about lerp(1, -1e18, 1e-18) ?
What about the similar lerp(1, -pow(2.0, 60), pow(2.0, -60)) ?
(I've been playing with these on godbolt because they have catastrophic cancellation)

What about them? I'd be glad to have a correctly rounded lerp (which would then _a fortiori_ be exact for 0 and 1), but I don't know of an algorithm with that property.

Your alternative implementation did make me realize that (for IEC 60559) the min/max step you omitted is needed only for non-default rounding modes; with round-to-nearest, the correction for t==1 (exactly) never introduces a non-monotonicity (because of monotonicity of the fundamental operations and the fact that x*std::nextafter(1,0)==std::nextafter(x,0) for all normal _x_).

[jorgbrown]

inf + NaN is inf.
Not on the implementations I have...

I stand corrected and you are right: inf + NaN is NaN. I swear I tried and got inf when I tried it the first time. I even checked 754-2008 and it agrees with you as well. https://godbolt.org/z/EAPqWf Sigh

[jorgbrown]

What about lerp(1, -1e18, 1e-18) ?
What about the similar lerp(1, -pow(2.0, 60), pow(2.0, -60)) ?

What about them?

one quintillionth of the way from zero to one quintillion is one.
And one quintillionth of the way from -1 to (one quintillion - 1) is zero.
But one quintillionth of the way from -1 to one quintillion is NOT zero. It is one quintillionth, because -1 + (X+1)/X is 1/X.

But lerp(1, -1e18, 1e-18) produces 0.0 using libc++, thanks to catastrophic cancellation.

If we use this definition for lerp:

double lerp(double a, double b, double frac) {
  return std::fma(-a, frac, std::fma(b, frac, a));
}

Then we get a more accurate representation for lerp(1, -1e18, 1e-18) : 7.2542424054621928e-17, or roughly 1e-18. (It's not very close to 1e-18 because 1e18 and 1e-18 aren't very precise.)

lerp(1, -pow(2.0, 60), pow(2.0, -60)) is more precise: 8.6736173798840355e-19, that is, pow(2.0, -60). See https://godbolt.org/z/STQp57

[BillyONeal]

We need to preserve the usual rules for NaNs, that being that if any input is a NaN, the result is also a NaN.

Is that behavior actually required anywhere? I don't see any blanket wording about it, but of course blanket wording is good at hiding. I do see that my proof-of-principle implementation doesn't have that property, but only because of its trick with std::min and std::max preferring their first argument.

I don't think the standard says anything about it because the standard doesn't say much about how floats are supposed to work, but the usual rule for anything floating point is NaN comes in NaN goes out. (As you said above, infectious to draw attention)