Use the Kahan-Babushka-Neumaier algorithm to improve the accuracy of sum().

FossilOrigin-Name: c63e26e705f5e967e14ef6aea8ce226548293ad8d25066069f29fa89673913d2

1 files changed, 72 insertions, 14 deletions

diff --git a/src/func.c b/src/func.c
index c6a49de61..beafb5431 100644
--- a/src/func.c
+++ b/src/func.c
@@ -1671,6 +1671,7 @@ static void loadExt(sqlite3_context *context, int argc, sqlite3_value **argv){
 typedef struct SumCtx SumCtx;
 struct SumCtx {
   double rSum;      /* Running sum as as a double */
+  double rErr;      /* Error term for Kahan-Babushka-Neumaier summation */
   i64 iSum;         /* Running sum as a signed integer */
   i64 cnt;          /* Number of elements summed */
   u8 approx;        /* True if any non-integer value was input to the sum */
@@ -1678,6 +1679,51 @@ struct SumCtx {
 };
 
 /*
+** Do one step of the Kahan-Babushka-Neumaier summation.
+**
+** https://en.wikipedia.org/wiki/Kahan_summation_algorithm
+**
+** Variables are marked "volatile" to defeat c89 x86 floating point
+** optimizations can mess up this algorithm.
+*/
+static void kahanBabuskaNeumaierStep(
+  volatile SumCtx *pSum,
+  volatile double r
+){
+  volatile double s = pSum->rSum;
+  volatile double t = s + r;
+  if( fabs(s) > fabs(r) ){
+    pSum->rErr += (s - t) + r;
+  }else{
+    pSum->rErr += (r - t) + s;
+  }
+  pSum->rSum = t;
+}
+
+/*
+** Add a (possibly large) integer to the running sum.
+*/
+static void kahanBabuskaNeumaierStepInt64(volatile SumCtx *pSum, i64 iVal){
+  volatile double rVal = (double)iVal;
+  kahanBabuskaNeumaierStep(pSum, rVal);
+  if( iVal<=-4503599627370496 || iVal>=+4503599627370496 ){
+    double rDiff = (double)(iVal - (i64)rVal);
+    kahanBabuskaNeumaierStep(pSum, rDiff);
+  }
+}
+
+/*
+** Initialize the Kahan-Babaska-Neumaier sum from a 64-bit integer
+*/
+static void kahanBabuskaNeumaierInit(
+  volatile SumCtx *pSum,
+  i64 iVal
+){
+  pSum->rSum = (double)iVal;
+  pSum->rErr = (double)(iVal - (i64)pSum->rSum);
+}
+
+/*
 ** Routines used to compute the sum, average, and total.
 **
 ** The SUM() function follows the (broken) SQL standard which means
@@ -1698,24 +1744,28 @@ static void sumStep(sqlite3_context *context, int argc, sqlite3_value **argv){
     p->cnt++;
     if( p->approx==0 ){
       if( type!=SQLITE_INTEGER ){
-        p->rSum = (double)p->iSum;
+        kahanBabuskaNeumaierInit(p, p->iSum);
         p->approx = 1;
-        p->rSum += sqlite3_value_double(argv[0]);
+        kahanBabuskaNeumaierStep(p, sqlite3_value_double(argv[0]));
       }else{
         i64 x = p->iSum;
         if( sqlite3AddInt64(&x, sqlite3_value_int64(argv[0]))==0 ){
           p->iSum = x;
         }else{
           p->ovrfl = 1;
-          p->rSum = (double)p->iSum;
+          kahanBabuskaNeumaierInit(p, p->iSum);
           p->approx = 1;
-          p->rSum += sqlite3_value_double(argv[0]);
+          kahanBabuskaNeumaierStep(p, sqlite3_value_double(argv[0]));
         }
       }
     }else{
-      if( type!=SQLITE_INTEGER ) p->ovrfl = 0;
       p->approx = 1;
-      p->rSum += sqlite3_value_double(argv[0]);
+      if( type==SQLITE_INTEGER ){
+        kahanBabuskaNeumaierStepInt64(p, sqlite3_value_int64(argv[0]));
+      }else{
+        p->ovrfl = 0;
+        kahanBabuskaNeumaierStep(p, sqlite3_value_double(argv[0]));
+      }
     }
   }
 }
@@ -1732,10 +1782,18 @@ static void sumInverse(sqlite3_context *context, int argc, sqlite3_value**argv){
   if( ALWAYS(p) && type!=SQLITE_NULL ){
     assert( p->cnt>0 );
     p->cnt--;
-    if( p->approx ){
-      p->rSum -= sqlite3_value_double(argv[0]);
-    }else{
+    if( !p->approx ){
       p->iSum -= sqlite3_value_int64(argv[0]);
+    }else if( type==SQLITE_INTEGER ){
+      i64 iVal = sqlite3_value_int64(argv[0]);
+      if( iVal!=SMALLEST_INT64 ){
+        kahanBabuskaNeumaierStepInt64(p, -iVal);
+      }else{
+        kahanBabuskaNeumaierStepInt64(p, LARGEST_INT64);
+        kahanBabuskaNeumaierStepInt64(p, 1);
+      }       
+    }else{
+      kahanBabuskaNeumaierStep(p, -sqlite3_value_double(argv[0]));
     }
   }
 }
@@ -1750,7 +1808,7 @@ static void sumFinalize(sqlite3_context *context){
       if( p->ovrfl ){
         sqlite3_result_error(context,"integer overflow",-1);
       }else{
-        sqlite3_result_double(context, p->rSum);
+        sqlite3_result_double(context, p->rSum+p->rErr);
       }
     }else{
       sqlite3_result_int64(context, p->iSum);
@@ -1763,9 +1821,9 @@ static void avgFinalize(sqlite3_context *context){
   if( p && p->cnt>0 ){
     double r;
     if( p->approx ){
-      r = p->rSum;
+      r = p->rSum+p->rErr;
     }else{
-      r = sqlite3RealToI64(p->iSum);
+      r = (double)(p->iSum);
     }
     sqlite3_result_double(context, r/(double)p->cnt);
   }
@@ -1776,9 +1834,9 @@ static void totalFinalize(sqlite3_context *context){
   p = sqlite3_aggregate_context(context, 0);
   if( p ){
     if( p->approx ){
-      r = p->rSum;
+      r = p->rSum+p->rErr;
     }else{
-      r = sqlite3RealToI64(p->iSum);
+      r = (double)(p->iSum);
     }
   }
   sqlite3_result_double(context, r);