database/sqlite
1
0
Fork 0
mirror of https://github.com/sqlite/sqlite.git synced 2025-10-21 11:13:54 +03:00
Code Activity

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

Browse Source 
Kahan-Babuska-Neumaier algorithm.

FossilOrigin-Name: ebc5edd3b10c1102b07b9fb0d6837266b81e55504ef883b9b8a7ad5e8ab29dd2
This commit is contained in:
drh
2023-07-06 13:19:10 +00:00
parent 3748b7329f
commit 45d7562843
3 changed files with 81 additions and 22 deletions
Download Patch File Download Diff File Expand all files Collapse all files

85
src/func.c
View File

@@ -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);