Parallel numeric ranges algorithms: Add identity interface use cases

[mhoemmen]

Add identity element interface sketch and use cases

Add use cases showing sketches of an interface for specifying an identity element.

Interface sketch

Users specify identity element with a named struct reduce_identity. Its default specialization looks like this.

template<class T>
struct reduce_identity {
  T value{};
};

The default identity is T{}. Users have two ways to use a nondefault value.

1. Specify the value: reduce_identity{value}
2. Specialize reduce_identity<T>

Option (1) supports the use case where the same reduction result type might be used with different binary operators. For example, users might have some parts of their code where they use double values with the normal floating-point algebra, and other parts where they use double values with the min-plus ("tropical") algebra. For the former case, the default reduce_identity<double> would suffice. For the latter case, users would specify the identity as reduce_identity{numeric_limits<double>::infinity()}.

Option (2) lets users specify the value at compile time if they want.

namespace impl {
  inline constexpr T some_value = /* value goes here */;
}
template<class T>
struct reduce_identity<T> : constant_wrapper<impl::some_value> {};

Should the default identity value depend on the binary operation?

It would be more user-friendly if the default identity value could depend on the binary operation. For example, the default identity value for std::multiplies<float>{} should be 1.0f. This would imitate OpenACC and OpenMP behavior.

One way to achieve this goal would be to have special cases for known combinations of binary operator type and reduction result type. For example, the identity for std::multiplies<float> would be 1.0f. The problem with this approach is that users often use the Standard Library binary operator types with their default template argument void. std::multiples<void> just means "use binary operator*." This would give incorrect results for user-defined reduction result types. For example, if users define a min-plus algebra number type MinPlus with underlying type float, its overloaded operator* would actually do addition on float values, and thus should have identity value 0.0f.

Another way to achieve this goal would be to template reduce_identity on the binary operation type BinaryOp, instead of (or in addition to) the reduction result type T. However, this would have at least two issues.

1. Users would not be able to specialize reduce_identity for lambdas, since every lambda declaration declares a unique type.
2. Binary operations with the same return type can themselves have different types.

Here is an example of (2).

float binary_op_1(const float& x, const float& y) {
    return x + y;
}
float binary_op_2(float x, float y) {
    return x + y;
}
static_assert(! std::is_same_v<
  decltype(binary_op_1),
  decltype(binary_op_2)>);
static_assert(std::is_invocable_r_v<float,
  decltype(binary_op_1), float, float>);
static_assert(std::is_invocable_r_v<float,
  decltype(binary_op_2), float, float>);

A result is that users might try to specialize on a function float(float, float), but then define the actual function as float(const float&, const float&). This would give them the wrong specialization of reduce_identity at the point of use without a compile-time error.

Users might try to solve this by writing a partial specialization of reduce_identity constrained on "any binary function invocable with two arguments U and V and returning T." However, correctly constraining U and V is difficult (range value type, reduction result type, or some intermediate type?). More importantly, the partial specialization would take effect for all binary operations on T, not just for the user's binary operation.

We conclude that it's not worth trying to make reduce_identity depend at compile time on the binary operation. The following decision tree explains how users can specify the identity value.

1. If T{} is the identity for the binary operation, rely on default reduce_identity<T>{}.value;

2. otherwise, if T{} is NOT the identity for the binary operation:

   a. If the binary operation is always T operator+(x,y) or std::plus<T>:

    i. specialize `reduce_identity<T>`, or
   
    ii. provide a value `reduce_identity<T>{value}` at the point of use;
   

   b. otherwise, provide a value reduce_identity<T>{value} at the point of use.

If users can specify the identity value, should they also specify an initial value?

The initial value parameter T init of C++17 std::reduce has two purposes.

1. It specifies the actual initial value, which lets users compute partial reductions.
2. It determines std::reduce's return type.

If users could specify the identity value, that would also determine std::reduce's result type. In fact, it would do so redundantly, because for C++17 std::reduce, the initial value type is the return type.

The ranges::fold_left algorithm does not need an identity value, for reasons we explain elsewhere in this proposal. However, it requires a (nonoptional) initial value. The initial value of ranges::fold_left determines the return type of ranges::fold_left, but is not necessarily the same as the return type. The actual return type depends on the binary operation, the range value type, and the initial value type.

It's an important feature of both ranges::fold_left and of our ranges::reduce that we propose here, to be able to deduce a default initial value type from the range type.

What should happen if users could specify both an initial value and an identity value, and those types differ? Should that be ill-formed?

Use cases

Ordinary floating-point plus, range of double.

auto id = reduce_identity<double>{};
auto result1 = ranges::reduce(exec_policy, range, plus{}, id);
// Default identity value is reduce_identity<double>{}
auto result2 = ranges::reduce(exec_policy, range, plus{});

Min reduction (e.g., for min-plus (tropical) algebra) on double.

auto binary_op = [] (double y) { return std::fmin(x, y); }
auto id = reduce_identity{numeric_limits<double>::infinity()}
// No initial value, just identity value
ranges::reduce(exec_policy, range, binary_op, id);

[akukanov]

I like the idea of providing both the customization point and the function argument.

But I think we can still try to explore customizations for the operation in addition to the type. The main motivation for that is to better support the standard binary functions out of the box. I agree with the outlined issues, but possibly the function argument option may serve as a "plan B" for cases where customization is unfeasible, like with lambdas.

Just thinking out loud: if a transparent function object, like std::multiplies<void>, is used as the reduction operator, we still know what types reduce will apply it to: the range value type, the initial value type (if given), and the local result type (which we can perhaps assume/require being the same as the return type of reduce). So, what if reduce_identity is specialized for std::multiplies<void> and std::ranges::range_value_t<R> - can the implementation use that?

And another thought: would reflection help (especially the presumably more powerful C++29 reflection)?

[mhoemmen]

@akukanov wrote:

So, what if reduce_identity is specialized for std::multiplies<void> and std::ranges::range_value_t<R> - can the implementation use that?

It's a good thought! Here are some design goals.

1. Users need to be able to specify a different binary operation and identity value for the same reduction result type (example: the usual algebra vs. the tropical algebra, both on double).
2. It would be user-friendly for *reduce and *scan defaults to behave like fold expression defaults. This means that algorithms need to know both range_value_t<R> and the binary operation in order to deduce a default identity value.
3. Users can use the identity value to specify a reduction result type that differs from range_value_t<R>.
4. This all needs to work for binary operations that are lambdas.
5. Users shouldn't have to specialize a class to specify a nondefault identity value.

Here are some use cases.

// Should this even work?
// Deduce binary op: std::plus{}
// Deduce identity: range_value_t<R>{}
auto result1 = std::ranges::reduce(exec_policy, range);

// Deduce identity: range_value_t<R>{}
auto result2 = std::ranges::reduce(exec_policy, range, std::plus{});

// If range_value_t<R> is arithmetic or std::complex,
// deduce identity: range_value_t<R>(1)
auto result3 = std::ranges::reduce(exec_policy, range, std::multiplies{});

// Should this even work?
// The "identity" is really -Inf if range_value_t<R> has that,
// but this is one of those cases where users might not like that behavior.
auto result4 = std::ranges::reduce(exec_policy, range, std::ranges::min);

// User explicitly opts into -Inf as the identity.
constexpr auto minus_Inf = -std::numeric_limits<float>::infinity();
auto result5 = std::ranges::reduce(exec_policy, range, std::ranges::min,
  reduce_identity{minus_Inf});

// range_value_t<R> is float, but identity value is double
// (even though it's otherwise the default value, zero).
// std::plus<void> should use operator()(double, double) -> double
auto result6 = std::ranges::reduce(exec_policy, range, std::plus{}, reduce_identity{0.0});

We could think about another design that combines the binary operation with the identity value into a single struct: reduce_operation{custom_binary_op, custom_value}. This would encapsulate the two things that belong together into a single argument of reduce. It would also make it easier for the implementation to deduce defaults.

auto result7 = std::ranges::reduce(exec_policy, range,
  reduce_operation{custom_binary_op, custom_value});

constexpr auto minus_Inf = -std::numeric_limits<float>::infinity();
auto result5 = std::ranges::reduce(exec_policy, range,
  reduce_operation{std::ranges::min, minus_Inf};

// range_value_t<R> is float
auto result6 = std::ranges::reduce(exec_policy, range,
  reduce_operation{std::plus{}, /* double */ 0.0};

// Deduce identity: range_value_t<R>{}
auto result2 = std::ranges::reduce(exec_policy, range, std::plus{});

Here is an implementation sketch.

// Keep this an aggregate if possible.
// We can constrain BinaryOp on being invocable with T, T,
// but we don't know range_value_t<R> at this point,
// so we don't know yet if the algorithm is well-formed with BinaryOp.
template<class BinaryOp, class T = invoke_result<BinaryOp, T, T>>
  requires semiregular<T> && invocable<BinaryOp, T, T>
struct reduce_operation {
  BinaryOp op;
  T value{};
};

// std::plus with arithmetic types -> T{}.
//
// We can't specialize for all T, because users might have
// defined std::plus<T> such that T{} is not the identity.
template<class U, class T>
  requires(
    (is_arithmetic_v<T>) &&
    (is_void_v<U> || (is_arithmetic_v<U> && is_convertible_v<U, T>))
  )
struct reduce_operation<std::plus<U>, T> {
  [[no_unique_address]] std::plus<U> op;
  static constexpr T value = T{};
};

// (std::plus or std::multiplies) with std::complex takes more effort,
// since e.g., std::complex<double> + float is not well-formed.

// std::multiplies with arithmetic types -> T(1).
//
// We can't specialize for all T, because users might have
// defined std::multiplies<T> such that T(1) is not the identity.
template<class T, class U>
  requires(
    (is_arithmetic_v<T>) &&
    (is_void_v<U> || (is_arithmetic_v<U> && is_convertible_v<U, T>))
  )
struct reduce_operation<T, std::multiplies<U>> {
  [[no_unique_address]] std::multiplies<U> op;
  static constexpr T value = T(1);
};

// ... more specializations for Standard binary function objects ...

Some design issues:

1. Users can't override default identity values for known binary operators and reduction result types. (Is this a feature or a bug?)
2. This design cannot constrain BinaryOp to be well-formed for the reduce algorithm, because it doesn't know range_value_t<R>.
3. By storing BinaryOp in a struct, it requires that BinaryOp be copy-constructible or move-constructible. That's fine for parallel algorithms, but not strictly necessary for non-parallel algorithms. (On the other hand, if you want non-copyable BinaryOp, maybe you should use fold_* instead?)