std::reduce seems to convert results to integers [duplicate]

Question

I'm comparing these two functions:

double polynomials(const vector<double> & coeffs, double x) {
    double sum = 0.0;
    double factor = 1.0;
    for (double coeff: coeffs) {
        sum += coeff * factor;
        factor *= x;
    }
    return sum;
}

and

double algorithm_polynomials(const vector<double> & coeffs, double x) {
    return reduce(execution::seq, cbegin(coeffs), end(coeffs), 0, [factor=1.0, x](double sum, double coeff) mutable {
        double curr_sum = sum + coeff * factor;
        factor *= x;
        return curr_sum;
    });
}

For integer coefficients and integer x, the results of these two functions are equivalent. But for real values, the results diverge. The ones given by the first function are correct, while the ones given by the second one are always integers (and not necessarily close to the real answer).

My guess is that each intermediate value is being converted to an integer before being summed, but I don't see why this should happen, or how I can fix it.

To run this code use

#include <numeric>
#include <vector>
#include <algorithm>
#include <execution>

using namespace std;

gcc 9 or greater with -std=c++17 flag and link to tbb from 2019 are required.

Edit: Thank you all for pointing out the initial value type that I had overlooked.

I feel a bit ashamed both for not seeing it, and for having that badly utilized reduce function (can't have it parallelized due to the mutable lambda).

To redeem myself, I add here a way to make it parallel...

class counter: public std::iterator<
        std::random_access_iterator_tag,   // iterator_category
        size_t,                      // value_type
        size_t,                      // difference_type
        const size_t*,               // pointer
        size_t                       // reference
>{
    size_t num = 0;
public:
    explicit counter(size_t _num) : num(_num) {}
    counter& operator++() {num += 1; return *this;}
    counter operator++(int) {counter retval = *this; ++(*this); return retval;}
    bool operator==(counter other) const {return num == other.num;}
    bool operator!=(counter other) const {return !(*this == other);}
    counter& operator+=(size_t i) { num += i; return *this; }
    counter& operator-=(size_t i) { num -= i; return *this; }
    counter operator +(counter &other) const { return counter(num + other.num);}
    counter operator -(counter &other) const { return counter(num - other.num); }
    counter operator +(size_t i) const { return counter(num + i); }
    counter operator -(size_t i) const {return counter(num - i); }
    reference operator*() const {return num;}
};


double better_algorithm_polinomials(const vector<double> & coeffs, double x) {
    //this has the advantage of being easily parallelized
    return transform_reduce(execution::par, cbegin(coeffs), end(coeffs), counter(0), 0.0, plus{}, [x](double coeff, size_t index) { return coeff * pow<double>(x, index); });
}

Answer 1

The version of std::reduce() that you are calling:

template<class ExecutionPolicy, class ForwardIt, class T, class BinaryOp>
T reduce(ExecutionPolicy&& policy,
         ForwardIt first, ForwardIt last, T init, BinaryOp binary_op);

You can clearly see that the return value uses the same data type as the init parameter, which in your case is being deduced as an int, which is why the result is an int.

To make the return value be a double instead, simply change the literal 0 to 0.0 in the init parameter:

return reduce(execution::seq, cbegin(coeffs), cend(coeffs), 0.0, ...);