2/20/2008

02-20-08 - 2

It's annoying that one of the most common things you want to do with floats is also one of the least precise : summing or taking an average of an array. Say you have a bunch of values that are all similar, and you want the average and sdev. To get a good sdev you need a very accurate average. If you have say a million or more numbers, each one becomes tiny compared to the sum, and if you just do a normal sequential sum you make total shit. Of course if your numbers are floats you can do your sum in doubles and mostly not have a problem, but if you really need fine accuracy you need a better solution.

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.


template 
static 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:

old rants