Solving the max-subarray problem from CLRS in C++

I'm working through Algorithms Cormen et al to learn more about algorithms and modern C++. I've implemented the maximum subarray problem and would like some feedback. I'd really like to use the modern additions to C++ effectively, so please let me know where things can be improved.

One of my concerns was using std::numeric_limits::infinity in the code. Would this be used in production code? How would one handle something like infinity otherwise? I often see Dijkstra's algorithm described in pseudocode this way (using infinity), but is there something smarter/more standard to use?

Another concern was using tuples to return multiple values from a function. I find that I then have to use std::get quite a lot to access the individual elements of the tuple. Is this acceptable?

#include <iostream>
#include <tuple>
#include <vector>
#include <limits>
#include <math.h>

std::tuple <int,int,double> max_crossing_subarray(std::vector<double> A, int low, int mid, int high);
std::tuple<int,int,double> max_subarray(std::vector<double> A, int low, int high);

int main(){

std::vector<double> A = {13, -3, -25, 20, -3, -16, -23, 18, 20, -7, 12, -5, -22, 15, -4, 7};
//std::vector<double> A = {1, -1, 9};

return 0;

}

std::tuple<int,int,double> max_subarray(std::vector<double> A, int low, int high){

std::tuple<int,int,double> ans;
std::tuple<int,int,double> left_sub;
std::tuple<int,int,double> right_sub;
std::tuple<int,int,double> crossing_sub;

if(high==low){
ans = {low,high,A[low]};
return ans;
}
else{
int mid = floor((low+high)/2);
left_sub = max_subarray(A, low, mid);
right_sub = max_subarray(A, mid+1, high);
crossing_sub = max_crossing_subarray(A, low, mid, high);
}

// is my left subarray the largest?
if(std::get<2>(left_sub)>std::get<2>(right_sub) & std::get<2>(left_sub)>std::get<2>(crossing_sub)){
ans = left_sub;
}
//is my right subarray largest?
else if (std::get<2>(right_sub)>std::get<2>(left_sub) & std::get<2>(right_sub)>std::get<2>(crossing_sub)){
ans = right_sub;
}
// my crossing subarray is largest
else{
ans = crossing_sub;
}

return ans;

}

std::tuple <int,int,double> max_crossing_subarray(std::vector<double> A, int low, int mid, int high){

double inf = std::numeric_limits<double>::infinity();  // is this wise?
double sum = 0;

int max_left = 0;
double left_sum = -inf;
for(int i=mid; i>=low; i--){
sum += A[i];
if(sum>left_sum){
left_sum = sum;
max_left = i;
}
}

int max_right = 0;
double right_sum = -inf;
sum = 0;
for(int j=mid+1; j<=high; j++){
sum += A[j];
if(sum > right_sum){
right_sum = sum;
max_right = j;
}

}

std::tuple<int,int,double> out = {max_left, max_right, left_sum+right_sum};

return out;

}


I'd really like to use the modern additions to C++ effectively, so please let me know where things can be improved.

C++ is suprisingly strong in expressing algorithms in terms of minimal type requirements quite "easily" with templates. To effectively use the standard library and its modern additions one should know that

C++ algorithms are written in terms of iterators (or ranges)

So I would like to C++-ify your implementation of max_subarray and show the usage of some new and cool features of C++17.

Your current version only works for std::vector and with little adaption to any array-like container. But what if you want to search a std::list or some container like std::unordered_map (with some projection onto keys or values)?

With that abstraction I wouldn't call the algorithm max_subarray but max_subrange. The fundamental algorithm is max_crossing_subrange and that is why I will show you how to C++ify this function

One of my concerns was using std::numeric_limits::infinity in the code. Would this be used in production code? How would one handle something like infinity otherwise?

You are right to question that, because it is not okay or generic in any way. What if you don't sum numbers but you are folding some lists instead in your accumulation of neighboured elements? What you do there is a simple search for a maximum and you could simply evaluate the first expression instead. Look as an example at std::max_element, to see it is done there.

Another concern was using tuples to return multiple values from a function. I find that I then have to use std::get quite a lot to access the individual elements of the tuple. Is this acceptable?

Yes, that is okay and should be optimized away completely.

Now lets start step-by-step with a generic mas_crossing_subrange algorithm.

Warning: This does not lead to the optimal algorithm (which is O(N)) but just shows how to express algorithms in a C++-way for your example.

Your max_crossing_subarray uses two loops which are essentially same but reversed in their direction. In fact you can express it as one function using std::reverse_iterator. Lets refactor a sub function which performs the loop. An Iterator version would look like this:

/// \brief Fixing the first position we search an end position which
///        maximises the sum from first to end.
template <
typename I,
typename T = value_type_t<I>
>
constexpr std::pair<I, T>
max_sum_not_empty(I first, I last, const T& init = T{})
noexcept
{
assert(first != last);
auto max = std::make_pair(first, *first);
auto current_sum = max.second;
for (auto iter{std::next(first)}; iter != last; ++iter) {
current_sum = current_sum + *iter;
if (max < current_sum) {
max = std::make_pair(iter, current_sum);
}
}
return max;
}


But this isn't totally generic either. What if we want to generalise the accumulation to not just sums but any given fold operator f: X x X -> X (we can not use f: Y x X -> Y in max_crossing) where summation and less-comparison are default behaviour. It would look like this:

template <
typename I,
typename T = value_type_t<I>,
typename F = std::plus<>,
typename P = std::less<>
>
constexpr std::pair<I, T> max_accumulation_not_empty(
I first, I last, const T& init = T{}, F op = F{}, P pred = P{})
noexcept (
std::is_nothrow_callable_v<F(const T&, reference_t<I>)> &&
std::is_nothrow_callable_v<P(const T&, const T&)>
)
{
assert(first != last);
auto max = std::make_pair(first, std::invoke(op, init, *first));
auto current_sum = max.second;
for (auto iter{std::next(first)}; iter != last; ++iter) {
current_sum = std::invoke(op, current_sum, *iter);
if (std::invoke(pred, max.second, current_sum)) {
max = std::make_pair(iter, current_sum);
}
}
return max;
}


The reversed algorithm can be written as

template <
typename I,
typename T = value_type_t<I>,
typename F = std::plus<>,
typename P = std::less<>
>
constexpr std::pair<I, T> reverse_max_accumulation_not_empty(
I first, I last, const T& init = T{}, F op = F{}, P pred = P{})
noexcept (
std::is_nothrow_callable_v<F(const T&, reference_t<I>)> &&
std::is_nothrow_callable_v<P(const T&, const T&)>
)
{
auto rfirst = std::make_reverse_iterator(last);
auto rlast = std::make_reverse_iterator(first);
auto [rbest, sum] = max_accumulation_not_empty(rfirst, rlast, init,
std::ref(op), std::ref(pred));
auto best = std::next(rbest).base();
return std::make_pair(best, sum);
}


Your max_crossing algorithm can be then given as

template <
typename I,
typename T = value_type_t<I>,
typename F = std::plus<>,
typename P = std::less<>
>
constexpr std::tuple<I, difference_type_t<I>, T>
max_crossing_subrange(I first, I mid, I last, const T& init = T{},
F op = F{}, P pred = P{})
noexcept (
std::is_nothrow_callable_v<F(const T&, reference_t<I>)> &&
std::is_nothrow_callable_v<P(const T&, const T&)>
)
{
assert(first != mid && mid != last);
auto [best_left, left_sum] = reverse_max_accumulation_not_empty(
first, mid, init, std::ref(op), std::ref(pred));
auto [best_right, right_sum] = max_accumulation_not_empty(
mid, last, init, std::ref(op), std::ref(pred));
auto length = std::distance(best_left, best_right);
// here comes a strong assumption: op(left_sum, right_sum) needs
// an operator op: X x X -> X
return {best_left, length, op(left_sum, right_sum)};
}


If this would lead to an optimal algorithm it would look like this in production code. As an exercise: Try to finish this and proof that the assertions always evaluates to true.

• Thanks so much for the comprehensive reply! I'll work on wrapping my head around all of this and implementing the whole method this way. – Matt Mar 6 '17 at 12:54