## [Octave-cvsupdate] SF.net SVN: octave:[6282] trunk/octave-forge/extra/NaN/src/sumskipnan_mex. cpp

 [Octave-cvsupdate] SF.net SVN: octave:[6282] trunk/octave-forge/extra/NaN/src/sumskipnan_mex. cpp From: - 2009-09-30 08:18:18 ```Revision: 6282 http://octave.svn.sourceforge.net/octave/?rev=6282&view=rev Author: schloegl Date: 2009-09-30 08:18:02 +0000 (Wed, 30 Sep 2009) Log Message: ----------- increase accuracy by using Kahan's summation formula - this reduces the error by a factor of N; first tests show a penalty of about 40% in terms of speed; currently, this is only implemented for stride=1 (summation along 1st dimension) Modified Paths: -------------- trunk/octave-forge/extra/NaN/src/sumskipnan_mex.cpp Modified: trunk/octave-forge/extra/NaN/src/sumskipnan_mex.cpp =================================================================== --- trunk/octave-forge/extra/NaN/src/sumskipnan_mex.cpp 2009-09-29 21:59:09 UTC (rev 6281) +++ trunk/octave-forge/extra/NaN/src/sumskipnan_mex.cpp 2009-09-30 08:18:02 UTC (rev 6282) @@ -22,6 +22,11 @@ // usage: // [o,count,SSQ] = sumskipnan_mex(x,DIM,flag,W); // +// SUMSKIPNAN uses two techniques to reduce errors: +// 1) long double (80bit) instead of 64-bit double is used internally +// 2) The Kahan Summation formula is used to reduce the error margin from N*eps to 2*eps +// The latter is only implemented in case of stride=1 (column vectors only, summation along 1st dimension). +// // Input: // - x data array // - DIM (optional) dimension to sum @@ -50,10 +55,11 @@ inline int __sumskipnan2w__(double *data, size_t Ni, double *s, double *No, char *flag_anyISNAN, double *W); inline int __sumskipnan3w__(double *data, size_t Ni, double *s, double *s2, double *No, char *flag_anyISNAN, double *W); +inline int __sumskipnan2we__(double *data, size_t Ni, double *s, double *No, char *flag_anyISNAN, double *W); +inline int __sumskipnan3we__(double *data, size_t Ni, double *s, double *s2, double *No, char *flag_anyISNAN, double *W); //#define NO_FLAG - #ifdef tmwtypes_h #if (MX_API_VER<0x07020000) typedef int mwSize; @@ -174,7 +180,7 @@ if (D1==1) { double count; - __sumskipnan2w__(LInput+ix1, D2, LOutputSum+ix0, &count, &flag_isNaN, W); + __sumskipnan2we__(LInput+ix1, D2, LOutputSum+ix0, &count, &flag_isNaN, W); } else { for (j=0; j