Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Finish P0811R3 midpoint and lerp #1048

Merged
merged 12 commits into from
Jul 30, 2020
112 changes: 43 additions & 69 deletions stl/inc/cmath
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,10 @@
#include <cstdlib>
#include <xtr1common>

#if _HAS_CXX20
#include <xutility>
#endif // _HAS_CXX20

#pragma pack(push, _CRT_PACKING)
#pragma warning(push, _STL_WARNING_LEVEL)
#pragma warning(disable : _STL_DISABLED_WARNINGS)
Expand Down Expand Up @@ -1238,13 +1242,12 @@ _NODISCARD auto hypot(const _Ty1 _Dx, const _Ty2 _Dy, const _Ty3 _Dz) {

#if _HAS_CXX20
// FUNCTION lerp
// TRANSITION, P0553: lerp is not yet constexpr
template <class _Ty>
_NODISCARD /* constexpr */ _Ty _Common_lerp(const _Ty _ArgA, const _Ty _ArgB, const _Ty _ArgT) noexcept {
_NODISCARD constexpr _Ty _Common_lerp(const _Ty _ArgA, const _Ty _ArgB, const _Ty _ArgT) noexcept {
// on a line intersecting {(0.0, _ArgA), (1.0, _ArgB)}, return the Y value for X == _ArgT

const int _Finite_mask = (int{isfinite(_ArgA)} << 2) | (int{isfinite(_ArgB)} << 1) | int{isfinite(_ArgT)};
if (_Finite_mask == 0b111) {
const bool _T_is_finite = _STD _Is_finite(_ArgT);
if (_T_is_finite && _STD _Is_finite(_ArgA) && _STD _Is_finite(_ArgB)) {
// 99% case, put it first; this block comes from P0811R3
if ((_ArgA <= 0 && _ArgB >= 0) || (_ArgA >= 0 && _ArgB <= 0)) {
// exact, monotonic, bounded, determinate, and (for _ArgA == _ArgB == 0) consistent:
Expand Down Expand Up @@ -1272,96 +1275,67 @@ _NODISCARD /* constexpr */ _Ty _Common_lerp(const _Ty _ArgA, const _Ty _ArgB, co
return _Candidate;
}

if (isnan(_ArgA)) {
return _ArgA;
}

if (isnan(_ArgB)) {
return _ArgB;
}

if (isnan(_ArgT)) {
return _ArgT;
}

switch (_Finite_mask) {
case 0b000:
// All values are infinities
if (_ArgT >= 1) {
return _ArgB;
}

return _ArgA;
case 0b010:
case 0b100:
case 0b110:
// _ArgT is an infinity; return infinity in the "direction" of _ArgA and _ArgB
return _ArgT * (_ArgB - _ArgA);
case 0b001:
// Here _ArgA and _ArgB are infinities
if (_ArgA == _ArgB) {
// same sign, so T doesn't matter
return _ArgA;
}

// Opposite signs, choose the "infinity direction" according to T if it makes sense.
if (_ArgT <= 0) {
if (_STD is_constant_evaluated()) {
if (_STD _Is_nan(_ArgA)) {
statementreply marked this conversation as resolved.
Show resolved Hide resolved
return _ArgA;
}

if (_ArgT >= 1) {
if (_STD _Is_nan(_ArgB)) {
return _ArgB;
}

// Interpolating between infinities of opposite signs doesn't make sense, NaN
if constexpr (sizeof(_Ty) == sizeof(float)) {
return __builtin_nanf("0");
} else {
return __builtin_nan("0");
}
case 0b011:
// _ArgA is an infinity but _ArgB is not
if (_ArgT == 1) {
return _ArgB;
if (_STD _Is_nan(_ArgT)) {
return _ArgT;
}
} else {
if (_STD _Is_nan(_ArgA)) {
// raise FE_INVALID if at least one of _ArgA, _ArgB and _ArgT is signaling NaN
_Ty _Tmp = _ArgA;

if (_ArgT < 1) {
// towards the infinity, return it
return _ArgA;
}
if (_STD _Is_nan(_ArgB) || _STD _Is_nan(_ArgT)) {
_Tmp = _ArgB + _ArgT;
statementreply marked this conversation as resolved.
Show resolved Hide resolved
}

// away from the infinity
return -_ArgA;
case 0b101:
// _ArgA is finite and _ArgB is an infinity
if (_ArgT == 0) {
return _ArgA;
return _ArgA + _Tmp;
}

if (_ArgT > 0) {
// toward the infinity
return _ArgB;
if (_STD _Is_nan(_ArgB) || _STD _Is_nan(_ArgT)) {
return _ArgB + _ArgT;
}
}

return -_ArgB;
case 0b111: // impossible; handled in fast path
default:
_CSTD abort();
if (_T_is_finite) {
// _ArgT is finite, _ArgA and/or _ArgB is infinity
if (_ArgT < 0) {
// if _ArgT < 0: return infinity in the "direction" of _ArgA if that exists, NaN otherwise
return _ArgA - _ArgB;
} else if (_ArgT <= 1) {
// if _ArgT == 0: return _ArgA (infinity) if _ArgB is finite, NaN otherwise
// if 0 < _ArgT < 1: return infinity "between" _ArgA and _ArgB if that exists, NaN otherwise
// if _ArgT == 1: return _ArgB (infinity) if _ArgA is finite, NaN otherwise
return _ArgT * _ArgB + (1 - _ArgT) * _ArgA;
} else {
// if _ArgT > 1: return infinity in the "direction" of _ArgB if that exists, NaN otherwise
return _ArgB - _ArgA;
}
} else {
// _ArgT is an infinity; return infinity in the "direction" of _ArgA and _ArgB if that exists, NaN otherwise
return _ArgT * (_ArgB - _ArgA);
}
}

// As of 2019-06-17 it is unclear whether the "sufficient additional overloads" clause is intended to target lerp;
// LWG-3223 is pending.

_NODISCARD /* constexpr */ inline float lerp(const float _ArgA, const float _ArgB, const float _ArgT) noexcept {
_NODISCARD constexpr inline float lerp(const float _ArgA, const float _ArgB, const float _ArgT) noexcept {
return _Common_lerp(_ArgA, _ArgB, _ArgT);
}

_NODISCARD /* constexpr */ inline double lerp(const double _ArgA, const double _ArgB, const double _ArgT) noexcept {
_NODISCARD constexpr inline double lerp(const double _ArgA, const double _ArgB, const double _ArgT) noexcept {
return _Common_lerp(_ArgA, _ArgB, _ArgT);
}

_NODISCARD /* constexpr */ inline long double lerp(
_NODISCARD constexpr inline long double lerp(
const long double _ArgA, const long double _ArgB, const long double _ArgT) noexcept {
return _Common_lerp(_ArgA, _ArgB, _ArgT);
}
Expand Down
40 changes: 20 additions & 20 deletions stl/inc/numeric
Original file line number Diff line number Diff line change
Expand Up @@ -540,14 +540,9 @@ template <class _Arithmetic>
_NODISCARD constexpr auto _Abs_u(const _Arithmetic _Val) noexcept {
// computes absolute value of _Val (converting to an unsigned integer type if necessary to avoid overflow
// representing the negation of the minimum value)
if constexpr (is_floating_point_v<_Arithmetic>) {
// TRANSITION, P0553: this mishandles NaNs
if (_Val < 0) {
return -_Val;
}
static_assert(is_integral_v<_Arithmetic>);

return _Val;
} else if constexpr (is_signed_v<_Arithmetic>) {
if constexpr (is_signed_v<_Arithmetic>) {
using _Unsigned = make_unsigned_t<_Arithmetic>;
if (_Val < 0) {
// note static_cast to _Unsigned such that _Arithmetic == short returns unsigned short rather than int
Expand Down Expand Up @@ -619,9 +614,24 @@ _NODISCARD constexpr common_type_t<_Mt, _Nt> lcm(const _Mt _Mx, const _Nt _Nx) n
template <class _Ty, enable_if_t<is_arithmetic_v<_Ty> && !is_same_v<remove_cv_t<_Ty>, bool>, int> = 0>
_NODISCARD constexpr _Ty midpoint(const _Ty _Val1, const _Ty _Val2) noexcept {
if constexpr (is_floating_point_v<_Ty>) {
if (_STD is_constant_evaluated()) {
if (_STD _Is_nan(_Val1)) {
return _Val1;
}

if (_STD _Is_nan(_Val2)) {
return _Val2;
}
} else {
if (_STD _Is_nan(_Val1) || _STD _Is_nan(_Val2)) {
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Before I purposely put this test below the high limit tests before the NaN tests because we expect NaNs to be uncommon; you're fixing -fassociative-math and friends but it seems like that would still be fixed if you put this in the old location?

Copy link
Contributor Author

@statementreply statementreply Jul 23, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Performing the high limit test with <= raises undeserved FE_INVALID when the input value is NaN.

Alternative approaches:

  • Test with quiet comparison islessequal
  • Test with bit_cast tricks
  • Just ignore the issue of FE_INVALID

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I see, to raise floating point exceptions. Can you add a comment to places where you do seemingly needless operations to raise the right exceptions and add tests for that? Otherwise that's likely to break. (I understand that requires /fp:except which might be annoying to set up in our tests :( )

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hmmm actually it looks like we already have an /fp:except flavor tested:

PM_CL="/EHsc /MDd /D_ITERATOR_DEBUG_LEVEL=2 /std:c++latest /fp:except /Zc:preprocessor"

so you should just be able to do that.

// raise FE_INVALID if at least one of _Val1 and _Val2 is signaling NaN
return _Val1 + _Val2;
}
}

constexpr _Ty _High_limit = (numeric_limits<_Ty>::max)() / 2;
const auto _Val1_a = _Abs_u(_Val1);
const auto _Val2_a = _Abs_u(_Val2);
const auto _Val1_a = _Float_abs(_Val1);
const auto _Val2_a = _Float_abs(_Val2);
if (_Val1_a <= _High_limit && _Val2_a <= _High_limit) {
// _Val1 and _Val2 are small enough that _Val1 + _Val2 won't overflow

Expand All @@ -637,22 +647,12 @@ _NODISCARD constexpr _Ty midpoint(const _Ty _Val1, const _Ty _Val2) noexcept {
return (_Val1 + _Val2) / 2;
}

// TRANSITION, P0553: the next two branches handle NaNs but don't produce correct behavior under /fp:fast or
// -fassociative-math
if (_Val1 != _Val1) {
return _Val1;
}

if (_Val2 != _Val2) {
return _Val2;
}

// Here at least one of {_Val1, _Val2} has large magnitude.
// Therefore, if one of the values is too small to divide by 2 exactly, the small magnitude is much less than
// one ULP of the result, so we can add it directly without the potentially inexact division by 2.

// In the default rounding mode this less than one ULP difference will always be rounded away, so under
// /fp:precise or /fp:fast we could avoid these tests if we had some means of detecting it in the caller.
// /fp:fast we could avoid these tests if we had some means of detecting it in the caller.
constexpr _Ty _Low_limit = (numeric_limits<_Ty>::min)() * 2;
if (_Val1_a < _Low_limit) {
return _Val1 + _Val2 / 2;
Expand Down
Loading