Of course the solution is both easy and obvious. In the common/standard case that the numbers are all very close to the average, a simple recursively dividing sum works perfectly :

template < typename _t > _t recursive_sum(const _t * begin,const _t * end) { int count = end - begin; switch(count) { case 0: return 0.0; case 1: return begin[0]; case 2: return begin[0] + begin[1]; case 3: return begin[0] + begin[1] + begin[2]; case 4: return (begin[0] + begin[1]) + (begin[2] + begin[3]); default: { const _t * mid = begin + (count/2); return recursive_sum(begin,mid) + recursive_sum(mid,end); } } }

ADDENDUM : Okay this is a pretty dumb way to do it. It actually works fine and is plenty fast but the better way is this Kahan Summation thing. BTW the PDF linked at the bottom of that page is an awesome awesome doc on floating point.

(UPDATED 3-04-08) : Here's a version that roughly matches the signature of std::acumulate ; the type inference from accum is a little bit of a subtle thing that can bite you in code if you're not careful. BTW I flipped the sign compared to Wikipedia cuz it just makes it really obvious what's happening IMO.

templatestatic inline T kahan_accumulate(Iter begin, Iter end, T accum) { T err = 0; for(;begin != end;++begin) { T compensated = *begin + err; T tmp = accum + compensated; err = accum - tmp; err += compensated; accum = tmp; } return accum; } template static inline typename std::iterator_traits ::value_type kahan_accumulate(Iter begin, Iter end) { return kahan_accumulate(begin, end, typename std::iterator_traits ::value_type()); }

There are a bunch of ACM papers on this stuff; if you're serious you should go there. For our rough graphics kind of use this Kahan thing is going to be fine. The "Cascading Accumulators" method seems to be the awesome way to go but it's a bit complex to implement. Apparently Kahan also has less error than the recursive way, though in practice they're both tiny.

## No comments:

Post a Comment