use two-pass method in mean() (closes #219)

Author: eddelbuettel

ChangeLog

@@ -1,3 +1,7 @@
+2014-12-30  Dirk Eddelbuettel  <[email protected]>
+
+        * inst/include/Rcpp/sugar/functions/mean.h: Use two-pass method
+
 2014-12-29  Kevin Ushey  <[email protected]>
 
         * inst/include/Rcpp/macros/macros.h: reformat for legibility

inst/NEWS.Rd

@@ -20,6 +20,7 @@
       \code{List}, \code{RawVector} and \code{ExpressionVector}.
       \item Vectors now have a \code{Vector::const\_iterator} that is 'const correct'
       thanks to fix by Romain following bug report in rcpp-devel by Martyn Plummer 
+      \item The \code{mean()} sugar function now uses a more robust two-pass method.
     }
     \item Changes in Rcpp Attributes:
     \itemize{

inst/include/Rcpp/sugar/functions/mean.h

@@ -2,7 +2,7 @@
 //
 // mean.h: Rcpp R/C++ interface class library -- mean
 //
-// Copyright (C) 2011 Dirk Eddelbuettel and Romain Francois
+// Copyright (C) 2011 - 2014  Dirk Eddelbuettel and Romain Francois
 //
 // This file is part of Rcpp.
 //
@@ -28,23 +28,36 @@ namespace sugar{
 template <int RTYPE, bool NA, typename T>
 class Mean : public Lazy< typename Rcpp::traits::storage_type<RTYPE>::type , Mean<RTYPE,NA,T> > {
 public:
-	typedef typename Rcpp::VectorBase<RTYPE,NA,T> VEC_TYPE ;
-	typedef typename Rcpp::traits::storage_type<RTYPE>::type STORAGE ;
+    typedef typename Rcpp::VectorBase<RTYPE,NA,T> VEC_TYPE ;
+    typedef typename Rcpp::traits::storage_type<RTYPE>::type STORAGE ;
 
-	Mean( const VEC_TYPE& object_ ) : object(object_){}
+    Mean( const VEC_TYPE& object_ ) : object(object_){}
 
-	STORAGE get() const {
-		return sum(object).get() / object.size() ;
-	}
+    STORAGE get() const {
+        //return sum(object).get() / object.size() ;
+        NumericVector input = object;
+                
+        int n = input.size();           // double pass (as in summary.c)
+        long double s = std::accumulate(input.begin(), input.end(), 0.0L);
+        s /= n;
+        if (R_FINITE((double)s)) {  
+            long double t = 0.0;
+            for (int i = 0; i < n; i++) {
+                t += input[i] - s;
+            }
+            s += t/n;
+        }
+        return (double)s ;
+    }
 private:
-	const VEC_TYPE& object ;
+    const VEC_TYPE& object ;
 } ;
 
 } // sugar
 
 template <bool NA, typename T>
 inline sugar::Mean<REALSXP,NA,T> mean( const VectorBase<REALSXP,NA,T>& t){
-	return sugar::Mean<REALSXP,NA,T>( t ) ;
+    return sugar::Mean<REALSXP,NA,T>( t ) ;
 }