Experimental attempt to boost the accuracy of sum() using the

Kahan-Babuska-Neumaier algorithm.

FossilOrigin-Name: ebc5edd3b10c1102b07b9fb0d6837266b81e55504ef883b9b8a7ad5e8ab29dd2

diff --git a/manifest b/manifest
index 5aa49021d55c5acb246c9dac6f22a28e3993720a..3b64c316f5fca9a6c1724691416b82a1d0d208bb 100644 (file)
--- a/manifest
+++ b/manifest
@@ -1,5 +1,5 @@
-C Use\s"volatile"\sisntead\sof\s"#pragma"\sto\sget\sfloating\spoint\scalculations\nworking\scorrectly\swhen\scompiling\swith\sGCC\sfor\sx86\smachines.
-D 2023-07-06T00:55:06.161
+C Experimental\sattempt\sto\sboost\sthe\saccuracy\sof\ssum()\susing\sthe\nKahan-Babuska-Neumaier\salgorithm.
+D 2023-07-06T13:19:10.993
 F .fossil-settings/empty-dirs dbb81e8fc0401ac46a1491ab34a7f2c7c0452f2f06b54ebb845d024ca8283ef1
 F .fossil-settings/ignore-glob 35175cdfcf539b2318cb04a9901442804be81cd677d8b889fcc9149c21f239ea
 F LICENSE.md df5091916dbb40e6e9686186587125e1b2ff51f022cc334e886c19a0e9982724
@@ -590,7 +590,7 @@ F src/delete.c cd5f5cd06ed0b6a882ec1a8c2a0d73b3cecb28479ad19e9931c4706c5e2182be
 F src/expr.c 8d1656b65e26af3e34f78e947ac423f0d20c214ed25a67486e433bf16ca6b543
 F src/fault.c 460f3e55994363812d9d60844b2a6de88826e007
 F src/fkey.c a7fcbf7e66d14dbb73cf49f31489ebf66d0e6006c62b95246924a3bae9f37b36
-F src/func.c 6028c160f693bdd018b651b5468a0a8e790f4e01e200796916b2d10a5d3237aa
+F src/func.c baa50a89130c980f2f7da89918da76cfa72708074c72d2fe5b21308a59454d37
 F src/global.c a16553245e315ee0cda8f9b0bf744efef9dc99f86e9d77f58975ea58824ded92
 F src/hash.c 9ee4269fb1d6632a6fecfb9479c93a1f29271bddbbaf215dd60420bcb80c7220
 F src/hash.h 3340ab6e1d13e725571d7cee6d3e3135f0779a7d8e76a9ce0a85971fa3953c51
@@ -2043,9 +2043,11 @@ F vsixtest/vsixtest.tcl 6a9a6ab600c25a91a7acc6293828957a386a8a93
 F vsixtest/vsixtest.vcxproj.data 2ed517e100c66dc455b492e1a33350c1b20fbcdc
 F vsixtest/vsixtest.vcxproj.filters 37e51ffedcdb064aad6ff33b6148725226cd608e
 F vsixtest/vsixtest_TemporaryKey.pfx e5b1b036facdb453873e7084e1cae9102ccc67a0
-P 7b4c16731e7bf6f03f5adf4fcb2008c0b19be473fb1b90b405c217c08916586a 1d972a690fdc70ab40862bd38427d68b48e8802ddf8e5c301f2d58ce2178b6ec
-R edd6abc545b6b11f65503dd16b4a34c2
-T +closed 1d972a690fdc70ab40862bd38427d68b48e8802ddf8e5c301f2d58ce2178b6ec
+P 9427f42687ed6d97c474bf42d0c3e82d6f4b0075e74206adcb5699d72e32140e
+R d2e2f73c16e9792f0c7025c31bf09b42
+T *branch * kahan-babuska-neumaier-summation
+T *sym-kahan-babuska-neumaier-summation *
+T -sym-trunk *
 U drh
-Z 9d851cc5730e074609430b016c224723
+Z 0c513397b66b5a1b4aafb74e58b9b2cf
 # Remove this line to create a well-formed Fossil manifest.

diff --git a/manifest.uuid b/manifest.uuid
index 82d15a017280ed8f27d257e181452960b63a35ca..0dd4b76c06f32fa3e72876be5dcd9d849964dd32 100644 (file)
--- a/manifest.uuid
+++ b/manifest.uuid
@@ -1 +1 @@
-9427f42687ed6d97c474bf42d0c3e82d6f4b0075e74206adcb5699d72e32140e
\ No newline at end of file
+ebc5edd3b10c1102b07b9fb0d6837266b81e55504ef883b9b8a7ad5e8ab29dd2
\ No newline at end of file

diff --git a/src/func.c b/src/func.c
index c6a49de61342b1e2234224e0a57edef93a7de022..2073e9a299622534a32c0f49c600311773a35c8f 100644 (file)
--- a/src/func.c
+++ b/src/func.c
@@ -1671,12 +1671,57 @@ 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 */
   u8 ovrfl;         /* Integer overflow seen */
 };
 
+/*
+** 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 t = pSum->rSum + r;
+  if( fabs(pSum->rSum) > fabs(r) ){
+    pSum->rErr += (pSum->rSum - t) + r;
+  }else{
+    pSum->rErr += (r - t) + pSum->rSum;
+  }
+  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 - (double)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.
 **
@@ -1698,24 +1743,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 +1781,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/2);
+        kahanBabuskaNeumaierStepInt64(p, LARGEST_INT64/2+1);
+      }       
+    }else{
+      kahanBabuskaNeumaierStep(p, -sqlite3_value_double(argv[0]));
     }
   }
 }
@@ -1750,7 +1807,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 +1820,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 +1833,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);