I recently traced a bug of mine to the fact that cumsum(s)[length(s)] is not always exactly equal to sum(s). For example, x<-1/(12:14) sum(x) - cumsum(x)[3] => 2.8e-17 Floating-point addition is of course not exact, and in particular is not associative, so there are various possible reasons for this. Perhaps sum uses clever summing tricks to get more accurate results? In some quick experiments, it does seem to get more accurate results than cumsum. It might be worth documenting. -s
Check out sum.exact and cumsum.exact in the caTools package.> library(caTools)Loading required package: bitops> x <- 1/(12:14) > sum(x) - cumsum(x)[3][1] 2.775558e-17> sum.exact(x) - cumsum.exact(x)[3][1] 0 On Tue, Feb 17, 2009 at 5:12 PM, Stavros Macrakis <macrakis at alum.mit.edu> wrote:> I recently traced a bug of mine to the fact that cumsum(s)[length(s)] > is not always exactly equal to sum(s). > > For example, > > x<-1/(12:14) > sum(x) - cumsum(x)[3] => 2.8e-17 > > Floating-point addition is of course not exact, and in particular is > not associative, so there are various possible reasons for this. > Perhaps sum uses clever summing tricks to get more accurate results? > In some quick experiments, it does seem to get more accurate results > than cumsum. > > It might be worth documenting. > > -s > > ______________________________________________ > R-help at r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. >
>>>>> "SM" == Stavros Macrakis <macrakis at alum.mit.edu> >>>>> on Wed, 18 Feb 2009 10:00:40 -0500 writes:SM> Nice! Glad to hear it. It sounds as though it is still possible for SM> cumsum(x)[length(x)] to not be exactly equal to sum, though? Well, possible, probably yes, platform-dependently; However I vaguely remember that I didn't see one such case in the few experiments I did. Martin SM> On Wed, Feb 18, 2009 at 8:03 AM, Martin Maechler SM> <maechler at stat.math.ethz.ch> wrote: SM> ... >> o cumsum(x) and cumprod(x) for double precision x now use a long >> double accumulator where available and so more closely match >> sum() and prod() in potentially being more accurate.