Topic: partial_sum


Author: "James Kuyper Jr." <kuyper@wizard.net>
Date: Tue, 22 Jan 2002 09:39:26 CST
Raw View
Section 26.4.3 says that std::partial_sum(first,last, result, binary_op)
must do the following:

> Assigns to every element by an iterator j in the range [result, result+(last-
> first)) a value ... equal to
....
> binary_op(binary_op(..., binary_op(*first, *(first+1)), ...),
> *(first+(i-result)))

The standard specifies almost nothing about binary_op. Almost all that
we know about it is implied the requirement that the above be a legal
expression (with the ...'s suitably expanded, of course). In particular,
we don't know what the return type of binary_op() is. This in turn
implies that different calls to binary_op() in that nested sequence
could invoke different overloads of operator(), each of which could have
a different return type. All that seems to be required is that each of
those return types must be able to be assigned to *result.

The obvious implementation of partial_sum uses a variable to store the
intermediate results in that chain, and the one implementation I checked
does precisely that. It uses a variable of the value type of
InputIterator. It's not guaranteed that the return type of binary_op()
can even be assigned to such a variable. Even if it can be, there's no
gaurantee that doing so would  produce the same results as those
specified by 26.4.3. The only way I can see to get the right result with
such a binary_op() would involve recursion:

template<class _T, class _InputIterator, class _OutputIterator, class
_BinaryOperation>
_void __partial_sum(
    _T const _current_sum&,
    _InputIterator __first,
    _InputIterator const __last,
    _OutputIterator& __result,
    _BinaryOperation& __binary_op
){
    *__result++ = __current_sum;
    if(__last != __first)
        __partial_sum(binary_op(__current_sum, *__first),
            __first+1, __last, __result, __binary_op);
}

template<class _T, class _InputIterator, class _OutputIterator, class
_BinaryOperation>
OutputIterator partial_sum(
    _InputIterator __first,
    _InputIterator __last,
    _OutputIterator __result,
    _BinaryOperation __binary_op
){
    if(__first != __last)
    {
        *__result++ = *__first;
        if(__first+1 != __last)
            __partial_sum(*__first, __first+1, __last, __result,
__binary_op);
    }
    return result;
}


I doubt that this was the intended implementation. What was the intent?
What should an implementation use as the type for storing intermediate
results?

A practical example of the problem would be the following:

    #include <numeric>
    class double_sum {
        double operator()(float x, float y) { return (double)x+y);
        double operator()(double x, float y) { return x+y);
    };
    float array[N];
    // fill in array
    std::partial_sum(array, array+N, array, double_sum());

If partial_sum() were implemented exactly per the formula in 26.4.3,
this would use temporary doubles to accumulate the sums, reducing the
possible loss of precision that could occur due to repeated round-offs
in a long sum. However, the implementation I looked at would not only
fail to gain the benefits of using 'double', it would insert a wasteful
implicit conversion back to float at each step.

Of course, theoretically there could be much worse cases:

    template<int N> class step{
        template<int M> step& operator=(const step<M>&);
        // other details.
    };

    class do_step {
        step<1> operator()(const step<0>&, const step<0>&);
        template<int N> step<N+1>
        operator()(const step<N>&, const step<0>&);
        // more details
    };

    std::vector<step<0> > v;
    // Fill in v
    std::partial_sum<v.begin(), v.end(), v.begin(), do_step());

Even my recursive version couldn't handle this one properly, because the
template recursion in __partial_sum() won't terminate (even though the
function-call recursion does terminate).

---
[ comp.std.c++ is moderated.  To submit articles, try just posting with ]
[ your news-reader.  If that fails, use mailto:std-c++@ncar.ucar.edu    ]
[              --- Please see the FAQ before posting. ---               ]
[ FAQ: http://www.research.att.com/~austern/csc/faq.html                ]