Audience: LWG; WG14
S. Davis Herring <>
Los Alamos National Laboratory
February 22, 2019

# History

Changes in r3:

• Made integer `midpoint` a template.
• Fixed example implementation of it and discussion.
• Updated to modern library description conventions.
• Added missing `constexpr` in wording.

Changes in r2:

• Renamed `linear` to `lerp` per LEWG.
• Added `constexpr` and `noexcept` where appropriate.
• Expressed sign computations without multiplication or subtraction for robustness.
• Prohibited spurious NaN results.

Changes in r1:

• Renamed `mid` to `midpoint` and `lint` to `linear` after reflector discussion.
• Extended monotonicity to certain cases with infinite t.

# Introduction

The simple problem of computing a value between two other values is surprisingly subtle in general. This paper proposes library functions to compute the midpoint between two integer, floating-point, or pointer values, as well as a more general routine that interpolates (or extrapolates) between two floating-point values.

These utilities are a pure library extension. With the exception of the pointer versions, they would be valuable and straightforward to implement in C (perhaps via type-generic macros); it would be appropriate to consult with WG14 about including them in the C standard library.

# Midpoint

## Integers

Computing the (integer) midpoint of two integer values via

``(a+b)/2``

can cause overflow for signed or unsigned integers as well as for floating-point values. Java’s binary search implementation had this integer overflow bug for nearly a decade, and Mozilla had the same issue in its JavaScript implementation.

The standard alternative

``a+(b-a)/2``

works for signed integers with the same sign (even if `b<a`), but can overflow if they have different signs. The modular arithmetic of unsigned integers does not produce the value expected for `b<a` because the division inherent to midpoint is not native there; it instead produces the value halfway between a and the smallest modular equivalent to b that is no smaller.

Using the C++20 conversion from unsigned to signed that preserves bit patterns, the library can provide the simple implementation

``````constexpr Integer midpoint(Integer a, Integer b) noexcept {
using U = make_unsigned_t<Integer>;
return a>b ? a-(U(a)-b)/2 : a+(U(b)-a)/2;
}``````

that works for signed or unsigned `Integer`. Note that when `b==a+1` or `b==a-1` (without overflow), the result is a because of truncating division. This can be exploited to round half-integers up or down by supplying `a>=b` or `a<=b` respectively. (The simple `(a+b)/2` always truncates half-integers towards 0, yielding `min(a,b)` when they differ by 1.)

## Pointers

When an array is partitioned via pointer arithmetic (for binary search or parallelization, for example), the array size must not exceed `PTRDIFF_MAX` to avoid undefined behavior ([expr.add]/5). The library can also provide a function template

``````template<class T>
constexpr T* midpoint(T *a, T *b);``````

which is straightforward on common architectures but, it seems, cannot be implemented portably and efficiently in C++. As with integers, when the midpoint lies between two pointer values the one closer to `a` is chosen; for the usual case of `a<b`, this is compatible with the usual half-open ranges by selecting `a` when `[a,b)` is `[a,a+1)`.

## Floating-point types

Each of the midpoint formulas above can cause overflow for floating-point types; the latter can also produce results that are not correctly rounded (by rounding in the subtraction and the addition). A third choice

``a/2+b/2``

prevents overflow but is not correctly rounded for subnormal inputs (due to rounding each separately). The library can easily provide overloads of `midpoint` for floating-point types that is correctly rounded for all inputs by switching between these strategies:

``````Float midpoint(Float a, Float b)
{return isnormal(a) && isnormal(b) ? a/2+b/2 : (a+b)/2;}``````

(P0533R2 proposes making `isnormal` `constexpr`, in which case this implementation could be.)

## Naming

The name `mean` is superior for the floating-point case, but for the other types (and the common application to binary search) the name `midpoint` used above is more appropriate. It would be possible to use both names (despite the possible confusion with the use of `mean` in [rand.dist]), but as it aims to replace the simple expression `a+(b-a)/2` regardless of type, a single name seems best.

# Linear interpolation

Both obvious approaches used in published implementations of floating-point linear interpolation have issues:

1. `a+t*(b-a)` does not in general reproduce b when `t==1`, and can overflow if a and b have the largest exponent and opposite signs.
2. `t*b+(1-t)*a` is not monotonic in general (unless the product ab≤0).

Lacking an obvious, efficient means of obtaining a correctly rounded overall result, the goal is instead to guarantee the basic properties of

1. exactness: `lerp(a,b,0)==a && lerp(a,b,1)==b`
2. monotonicity: `cmp(lerp(a,b,t2),lerp(a,b,t1)) * cmp(t2,t1) * cmp(b,a) >= 0`, where `cmp` is an arithmetic three-way comparison function
3. determinacy: result of NaN only for `lerp(a,a,INFINITY)`
4. boundedness: `t<0 || t>1 || isfinite(lerp(a,b,t))`
5. consistency: `lerp(a,a,t)==a`

given that each argument `isfinite` (`!isnan` is sufficient for for monotonicity). These properties are useful in proofs of correctness of algorithms based on `lerp`. When interpolating, consistency follows from the other properties, but it and monotonicity apply even when extrapolating.

To demonstrate the feasibility of satisfying these properties, a simple implementation that provides all of them is given here:

``````constexpr Float lerp(Float a, Float b, Float t) {
// Exact, monotonic, bounded, determinate, and (for a=b=0) consistent:
if(a<=0 && b>=0 || a>=0 && b<=0) return t*b + (1-t)*a;

if(t==1) return b;                        // exact
// Exact at t=0, monotonic except near t=1,
// bounded, determinate, and consistent:
const Float x = a + t*(b-a);
return t>1 == b>a ? max(b,x) : min(b,x);  // monotonic near t=1
}``````

Putting `b` first in the `min`/`max` prefers it to another equal value (i.e., `-0.` vs. `+0.`), which avoids returning `-0.` for `t==1` but `+0.` for other nearby values of t.

Whether it uses this implementation or not, the library can provide a function satisfying these mathematical properties.

## Naming

The name `lerp` is very widely used in the graphics community (including animation and gaming), although in some uses t is restricted to [0,1].

# Wording

Relative to N4800.

## Feature-test macro

Add at the appropriate place in Table 36:

``__cpp_lib_interpolate           ...  <cmath> <numeric>``

## Midpoint

Add to the end of the synopsis in [numeric.ops.overview]:

``````  // 24.9.14, least common multiple
template<class M, class N>
constexpr common_type_t<M,N> lcm(M m, N n);``````

`  ``// 24.9.15, midpoint`
`  ``template<class T>`
`    ``constexpr T midpoint(T a, T b) noexcept;`
`  ``template<class T>`
`    ``constexpr T* midpoint(T* a, T* b);`

`}`

Add a new subsection after [numeric.ops.lcm]:

### 24.9.15 Midpoint [numeric.ops.midpoint]

``````template<class T>
constexpr T midpoint(T a, T b) noexcept;``````
1. Constraints: `T` is an arithmetic type other than `bool`.
2. Returns: Half the sum of `a` and `b`. If `T` is an integer type and the sum is odd, the result is rounded towards `a`.
3. Remarks: No overflow occurs. If `T` is a floating-point type, at most one inexact operation occurs.
``````template<class T>
constexpr T* midpoint(T* a, T* b);``````
1. Constraints: `T` is a complete object type.
2. Expects: `a` and `b` point to, respectively, elements `x[`i`]` and `x[`j`]` of the same array object `x`. [ Note: An object that is not an array element is considered to belong to a single-element array for this purpose; see [expr.unary.op]. A pointer past the last element of an array x of n elements is considered to be equivalent to a pointer to a hypothetical element x[n] for this purpose; see [basic.compound]. — end note ]
3. Returns: A pointer to `x[`i+(j-i)/2`]`, where the result of the division is truncated towards zero.

## Linear interpolation

Add to the synopsis in [cmath.syn]:

``long double fmal(long double x, long double y, long double z);``

`// 25.8.4, linear interpolation`
`constexpr float lerp(float a, float b, float t);`
`constexpr double lerp(double a, double b, double t);`
`constexpr long double lerp(long double a, long double b, long double t);`

``// 25.8.4, classification / comparison functions``

Add a new subsection after [c.math.hypot3]:

### 25.8.4 Linear interpolation [c.math.lerp]

``````constexpr float lerp(float a, float b, float t);
constexpr double lerp(double a, double b, double t);
constexpr long double lerp(long double a, long double b, long double t);``````
1. Returns: a+t(b-a).
2. Remarks: Let `r` be the value returned. If `isfinite(a) && isfinite(b)`, then:
1. If `t==0`, then `r==a`.
2. If `t==1`, then `r==b`.
3. If `t>=0 && t<=1`, then `isfinite(r)`.
4. If `isfinite(t) && a==b`, then `r==a`.
5. If `isfinite(t) || !isnan(t) && b-a!=0`, then `!isnan(r)`.
3. Let CMP(`x`,`y`) be 1 if `x>y`, −1 if `x<y`, and 0 otherwise. For any `t1` and `t2`, the product of CMP(`lerp(a,b,t2)`,`lerp(a,b,t1)`), CMP(`t2`,`t1`), and CMP(`b`,`a`) is non-negative.

# Acknowledgments

Thanks to Howard Hinnant and Dan Sunderland for debugging the suggested implementation of `midpoint`, and to Dan and Jeff Garland for extra wording review.