Enhance the SUM() aggregate (and related AVG() and TOTAL()) so that the running

sum is accurate to about 100 bits.

FossilOrigin-Name: a915f15a916af698e0cef46c8b3e7ed11bda19349179d2d414073cd39c4cce24

6 files changed, 104 insertions, 32 deletions

diff --git a/manifest b/manifest
index 8cd8dc169..95f091d73 100644
--- a/manifest
+++ b/manifest
@@ -1,5 +1,5 @@
-C Improved\saccuracy\sof\sfloating-point\sconversion\sconstants\sas\ssuggested\sby\n[forum:/info/598d32f6135c41c1|forum\spost\s598d32f6135c41c1].
-D 2023-06-28T11:46:28.098
+C Enhance\sthe\sSUM()\saggregate\s(and\srelated\sAVG()\sand\sTOTAL())\sso\sthat\sthe\srunning\nsum\sis\saccurate\sto\sabout\s100\sbits.
+D 2023-06-28T12:02:48.388
 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 ee064b524ac7077939febc2d57a5b4243c9bf0ab69893b0d336ec15d101b96b6
+F src/func.c 8ada46d362a153b8dfb5385c3aaa9a7d75ebf306b33d2e663aa03920126a1bc3
 F src/global.c bd0892ade7289f6e20bff44c07d06371f2ff9b53cea359e7854b9b72f65adc30
 F src/hash.c 9ee4269fb1d6632a6fecfb9479c93a1f29271bddbbaf215dd60420bcb80c7220
 F src/hash.h 3340ab6e1d13e725571d7cee6d3e3135f0779a7d8e76a9ce0a85971fa3953c51
@@ -642,7 +642,7 @@ F src/shell.c.in 5f07ea7ca4262872fcf07ce623a0be7c003f11834b80501874f039319aaf61b
 F src/sqlite.h.in 3076d78836b6dac53b3ab0875fc8fd15bca8077aad4d33c85336e05af6aef8c7
 F src/sqlite3.rc 5121c9e10c3964d5755191c80dd1180c122fc3a8
 F src/sqlite3ext.h da473ce2b3d0ae407a6300c4a164589b9a6bfdbec9462688a8593ff16f3bb6e4
-F src/sqliteInt.h 44353825a42fec1f2740a2a54de328c46d3dc3014be139341d42b4ad9ca2b7d9
+F src/sqliteInt.h f7e904f7fdeccfd5606ca4457122e647dcbaf307ed8615ac4865f8b5f536f77b
 F src/sqliteLimit.h 33b1c9baba578d34efe7dfdb43193b366111cdf41476b1e82699e14c11ee1fb6
 F src/status.c 160c445d7d28c984a0eae38c144f6419311ed3eace59b44ac6dafc20db4af749
 F src/table.c 0f141b58a16de7e2fbe81c308379e7279f4c6b50eb08efeec5892794a0ba30d1
@@ -705,7 +705,7 @@ F src/trigger.c ad6ab9452715fa9a8075442e15196022275b414b9141b566af8cdb7a1605f2b0
 F src/update.c 0aa36561167a7c40d01163238c297297962f31a15a8d742216b3c37cdf25f731
 F src/upsert.c 5303dc6c518fa7d4b280ec65170f465c7a70b7ac2b22491598f6d0b4875b3145
 F src/utf.c ee39565f0843775cc2c81135751ddd93eceb91a673ea2c57f61c76f288b041a0
-F src/util.c 4264102045fdb36e9af3ff361e390a5f7a76342a2bd7069e55d8ad332026d6b5
+F src/util.c a6b41c67ff2a5379b46b77e587b2e0adb2d2ddcc5669691674ca4d28e2755ae4
 F src/vacuum.c 604fcdaebe76f3497c855afcbf91b8fa5046b32de3045bab89cc008d68e40104
 F src/vdbe.c c993304c609326cf625b4ad30cbb0e15a3f64c941cf2c9713d0c360b4abbaa98
 F src/vdbe.h 41485521f68e9437fdb7ec4a90f9d86ab294e9bb8281e33b235915e29122cfc0
@@ -1107,7 +1107,7 @@ F test/fts4umlaut.test fcaca4471de7e78c9d1f7e8976e3e8704d7d8ad979d57a739d00f3f75
 F test/fts4unicode.test 82a9c16b68ba2f358a856226bb2ee02f81583797bc4744061c54401bf1a0f4c9
 F test/fts4upfrom.test f25835162c989dffd5e2ef91ec24c4848cc9973093e2d492d1c7b32afac1b49d
 F test/full.test 6b3c8fb43c6beab6b95438c1675374b95fab245d
-F test/func.test cdebc080122807254173b942bd0aa9464eb859f7a057a0791136072462fb5bcb
+F test/func.test 14ec60fd0cc4182e285e403b91e916e6c45597786347d9d5e59415cb4b130942
 F test/func2.test 772d66227e4e6684b86053302e2d74a2500e1e0f
 F test/func3.test 600a632c305a88f3946d38f9a51efe145c989b2e13bd2b2a488db47fe76bab6a
 F test/func4.test 2285fb5792d593fef442358763f0fd9de806eda47dbc7a5934df57ffdc484c31
@@ -2041,8 +2041,9 @@ F vsixtest/vsixtest.tcl 6a9a6ab600c25a91a7acc6293828957a386a8a93
 F vsixtest/vsixtest.vcxproj.data 2ed517e100c66dc455b492e1a33350c1b20fbcdc
 F vsixtest/vsixtest.vcxproj.filters 37e51ffedcdb064aad6ff33b6148725226cd608e
 F vsixtest/vsixtest_TemporaryKey.pfx e5b1b036facdb453873e7084e1cae9102ccc67a0
-P 9d7813ddcf67f1256e7331c5f390cb122b3f7e35cf46e0b33b610216a5396a00
-R 169687b27435f028f9845b8b7f3b29b8
+P 4943e8a1819e189747eefc414d02c0485e1620deff9cf92664295b21a8a9a83c c66ef2440e4e6c7aa17b50e5a29e543713ccab69aa0a415bac29b35b6116504a
+R 51bf89677ed17bae2e4a623d6d52a0c3
+T +closed c66ef2440e4e6c7aa17b50e5a29e543713ccab69aa0a415bac29b35b6116504a
 U drh
-Z 4fc8fce670e6ce62ea8a2e6be22b995e
+Z 8ab32d51087d44c8fc250f900dd8d446
 # Remove this line to create a well-formed Fossil manifest.
diff --git a/manifest.uuid b/manifest.uuid
index d18d43c4d..749907240 100644
--- a/manifest.uuid
+++ b/manifest.uuid
@@ -1 +1 @@
-4943e8a1819e189747eefc414d02c0485e1620deff9cf92664295b21a8a9a83c
\ No newline at end of file
+a915f15a916af698e0cef46c8b3e7ed11bda19349179d2d414073cd39c4cce24
\ No newline at end of file
diff --git a/src/func.c b/src/func.c
index f0695ac53..c505c37d6 100644
--- a/src/func.c
+++ b/src/func.c
@@ -1670,11 +1670,9 @@ static void loadExt(sqlite3_context *context, int argc, sqlite3_value **argv){
 */
 typedef struct SumCtx SumCtx;
 struct SumCtx {
-  double rSum;      /* Floating point sum */
-  i64 iSum;         /* Integer sum */  
+  double rSum[2];   /* Running sum as a Dekker double-double */
   i64 cnt;          /* Number of elements summed */
-  u8 overflow;      /* True if integer overflow seen */
-  u8 approx;        /* True if non-integer value was input to the sum */
+  u8 approx;        /* True if any non-integer value was input to the sum */
 };
 
 /*
@@ -1695,17 +1693,17 @@ static void sumStep(sqlite3_context *context, int argc, sqlite3_value **argv){
   p = sqlite3_aggregate_context(context, sizeof(*p));
   type = sqlite3_value_numeric_type(argv[0]);
   if( p && type!=SQLITE_NULL ){
+    double y[2];
     p->cnt++;
     if( type==SQLITE_INTEGER ){
       i64 v = sqlite3_value_int64(argv[0]);
-      p->rSum += v;
-      if( (p->approx|p->overflow)==0 && sqlite3AddInt64(&p->iSum, v) ){
-        p->approx = p->overflow = 1;
-      }
+      sqlite3DDFromInt(v, y);
     }else{
-      p->rSum += sqlite3_value_double(argv[0]);
+      y[0] = sqlite3_value_double(argv[0]);
+      y[1] = 0.0;
       p->approx = 1;
     }
+    sqlite3DDAdd(p->rSum[0], p->rSum[1], y[0], y[1], p->rSum);
   }
 }
 #ifndef SQLITE_OMIT_WINDOWFUNC
@@ -1719,16 +1717,17 @@ static void sumInverse(sqlite3_context *context, int argc, sqlite3_value**argv){
   /* p is always non-NULL because sumStep() will have been called first
   ** to initialize it */
   if( ALWAYS(p) && type!=SQLITE_NULL ){
+    double y[2];
     assert( p->cnt>0 );
     p->cnt--;
-    assert( type==SQLITE_INTEGER || p->approx );
-    if( type==SQLITE_INTEGER && p->approx==0 ){
+    if( type==SQLITE_INTEGER ){
       i64 v = sqlite3_value_int64(argv[0]);
-      p->rSum -= v;
-      p->iSum -= v;
+      sqlite3DDFromInt(v, y);
     }else{
-      p->rSum -= sqlite3_value_double(argv[0]);
+      y[0] = sqlite3_value_double(argv[0]);
+      y[1] = 0.0;
     }
+    sqlite3DDSub(p->rSum[0], p->rSum[1], y[0], y[1], p->rSum);
   }
 }
 #else
@@ -1738,12 +1737,18 @@ static void sumFinalize(sqlite3_context *context){
   SumCtx *p;
   p = sqlite3_aggregate_context(context, 0);
   if( p && p->cnt>0 ){
-    if( p->overflow ){
-      sqlite3_result_error(context,"integer overflow",-1);
-    }else if( p->approx ){
-      sqlite3_result_double(context, p->rSum);
+    if( p->approx ){
+      sqlite3_result_double(context, p->rSum[0]+p->rSum[1]);
     }else{
-      sqlite3_result_int64(context, p->iSum);
+      i64 v = (i64)p->rSum[0] + (i64)p->rSum[1];
+      double y[2], z[2];
+      sqlite3DDFromInt(v, y);
+      sqlite3DDSub(y[0], y[1], p->rSum[0], p->rSum[1], z);
+      if( z[0] + z[1] != 0.0 ){
+        sqlite3_result_error(context,"integer overflow",-1);
+      }else{
+        sqlite3_result_int64(context, v);
+      }
     }
   }
 }
@@ -1751,14 +1756,14 @@ static void avgFinalize(sqlite3_context *context){
   SumCtx *p;
   p = sqlite3_aggregate_context(context, 0);
   if( p && p->cnt>0 ){
-    sqlite3_result_double(context, p->rSum/(double)p->cnt);
+    sqlite3_result_double(context, (p->rSum[0]+p->rSum[1])/(double)p->cnt);
   }
 }
 static void totalFinalize(sqlite3_context *context){
   SumCtx *p;
   p = sqlite3_aggregate_context(context, 0);
   /* (double)0 In case of SQLITE_OMIT_FLOATING_POINT... */
-  sqlite3_result_double(context, p ? p->rSum : (double)0);
+  sqlite3_result_double(context, p ? p->rSum[0]+p->rSum[1] : (double)0);
 }
 
 /*
diff --git a/src/sqliteInt.h b/src/sqliteInt.h
index 91b78710c..332048f43 100644
--- a/src/sqliteInt.h
+++ b/src/sqliteInt.h
@@ -1234,6 +1234,7 @@ typedef struct Cte Cte;
 typedef struct CteUse CteUse;
 typedef struct Db Db;
 typedef struct DbFixer DbFixer;
+typedef struct DblDbl DblDbl;
 typedef struct Schema Schema;
 typedef struct Expr Expr;
 typedef struct ExprList ExprList;
@@ -5035,6 +5036,15 @@ int sqlite3FixSrcList(DbFixer*, SrcList*);
 int sqlite3FixSelect(DbFixer*, Select*);
 int sqlite3FixExpr(DbFixer*, Expr*);
 int sqlite3FixTriggerStep(DbFixer*, TriggerStep*);
+
+
+/* Representation of an extended precision floating point value.
+** The actual value is the sum r and rr.  See the dbldbl.c file.
+*/
+void sqlite3DDFromInt(i64,double*);
+void sqlite3DDAdd(double,double,double,double,double*);
+void sqlite3DDSub(double,double,double,double,double*);
+
 int sqlite3RealSameAsInt(double,sqlite3_int64);
 i64 sqlite3RealToI64(double);
 int sqlite3Int64ToText(i64,char*);
diff --git a/src/util.c b/src/util.c
index abd36eda8..256ec7c5c 100644
--- a/src/util.c
+++ b/src/util.c
@@ -1750,3 +1750,43 @@ int sqlite3VListNameToNum(VList *pIn, const char *zName, int nName){
  || defined(SQLITE_ENABLE_STMT_SCANSTATUS)
 # include "hwtime.h"
 #endif
+
+/***************************************************************************
+** Double-Double arithmetic.
+**
+** Reference:
+**   T. J. Dekker, "A Floating-Point Technique for Extending the
+**   Available Precision".  1971-07-26.
+*/
+
+/* Compute z = (i64)x */
+void sqlite3DDFromInt(i64 x, double *z){
+  z[0] = (double)x;
+  z[1] = (double)(x - (i64)z[0]);
+}
+
+/* Compute z = x + y */
+void sqlite3DDAdd(double x, double xx, double y, double yy, double *z){
+  double r, s;
+  r = x + y;
+  if( fabs(x)>fabs(y) ){
+    s = x - r + y + yy + xx;
+  }else{
+    s = y - r + x + xx + yy;
+  }
+  z[0] = r+s;
+  z[1] = r - z[0] + s;
+}
+
+/* Compute z = x - y */
+void sqlite3DDSub(double x, double xx, double y, double yy, double *z){
+  double r, s;
+  r = x - y;
+  if( fabs(x)>fabs(y) ){
+    s = x - r - y - yy + xx;
+  }else{
+    s = -y - r + x + xx - yy;
+  }
+  z[0] = r+s;
+  z[1] = r - z[0] + s;
+}
diff --git a/test/func.test b/test/func.test
index 6ba953647..d5ba5d02a 100644
--- a/test/func.test
+++ b/test/func.test
@@ -1543,6 +1543,22 @@ do_execsql_test func-36.110 {
   SELECT 123 ->> 456
 } {123->>456}
 
-
+# 2023-06-26
+# Enhanced precision of SUM().
+#
+reset_db
+do_execsql_test func-37.100 {
+  WITH c(x) AS (VALUES(9223372036854775807),(9223372036854775807),
+                      (123),(-9223372036854775807),(-9223372036854775807))
+  SELECT sum(x) FROM c;
+} {123}
+do_catchsql_test func-37.110 {
+  WITH c(x) AS (VALUES(9223372036854775807),(1))
+  SELECT sum(x) FROM c;
+} {1 {integer overflow}}
+do_catchsql_test func-37.120 {
+  WITH c(x) AS (VALUES(9223372036854775807),(100),(-101))
+  SELECT sum(x) FROM c;
+} {0 9223372036854775806}
 
 finish_test