Skip to content

Commit 0700965

Browse files
tglsfdctglsfdc
committed
Repair roundoff-error problem for stddev/variance results near zero,
per complaint from Kemin Zhou. Fix lack of precision in numeric stddev/variance.
1 parent 63cc56d commit 0700965

File tree

2 files changed

+98
-41
lines changed

2 files changed

+98
-41
lines changed

‎src/backend/utils/adt/float.c

Lines changed: 19 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -8,7 +8,7 @@
88
*
99
*
1010
* IDENTIFICATION
11-
* $Header: /cvsroot/pgsql/src/backend/utils/adt/float.c,v 1.77 2001/11/05 17:46:29 momjian Exp $
11+
* $Header: /cvsroot/pgsql/src/backend/utils/adt/float.c,v 1.78 2001/12/11 02:02:12 tgl Exp $
1212
*
1313
*-------------------------------------------------------------------------
1414
*/
@@ -1580,7 +1580,8 @@ float8_variance(PG_FUNCTION_ARGS)
15801580
float8 *transvalues;
15811581
float8 N,
15821582
sumX,
1583-
sumX2;
1583+
sumX2,
1584+
numerator;
15841585

15851586
transvalues = check_float8_array(transarray, "float8_variance");
15861587
N = transvalues[0];
@@ -1594,7 +1595,13 @@ float8_variance(PG_FUNCTION_ARGS)
15941595
if (N <= 1.0)
15951596
PG_RETURN_FLOAT8(0.0);
15961597

1597-
PG_RETURN_FLOAT8((N * sumX2 - sumX * sumX) / (N * (N - 1.0)));
1598+
numerator = N * sumX2 - sumX * sumX;
1599+
1600+
/* Watch out for roundoff error producing a negative numerator */
1601+
if (numerator <= 0.0)
1602+
PG_RETURN_FLOAT8(0.0);
1603+
1604+
PG_RETURN_FLOAT8(numerator / (N * (N - 1.0)));
15981605
}
15991606

16001607
Datum
@@ -1604,7 +1611,8 @@ float8_stddev(PG_FUNCTION_ARGS)
16041611
float8 *transvalues;
16051612
float8 N,
16061613
sumX,
1607-
sumX2;
1614+
sumX2,
1615+
numerator;
16081616

16091617
transvalues = check_float8_array(transarray, "float8_stddev");
16101618
N = transvalues[0];
@@ -1618,7 +1626,13 @@ float8_stddev(PG_FUNCTION_ARGS)
16181626
if (N <= 1.0)
16191627
PG_RETURN_FLOAT8(0.0);
16201628

1621-
PG_RETURN_FLOAT8(sqrt((N * sumX2 - sumX * sumX) / (N * (N - 1.0))));
1629+
numerator = N * sumX2 - sumX * sumX;
1630+
1631+
/* Watch out for roundoff error producing a negative numerator */
1632+
if (numerator <= 0.0)
1633+
PG_RETURN_FLOAT8(0.0);
1634+
1635+
PG_RETURN_FLOAT8(sqrt(numerator / (N * (N - 1.0))));
16221636
}
16231637

16241638

‎src/backend/utils/adt/numeric.c

Lines changed: 79 additions & 36 deletions
Original file line numberDiff line numberDiff line change
@@ -5,7 +5,7 @@
55
*
66
* 1998 Jan Wieck
77
*
8-
* $Header: /cvsroot/pgsql/src/backend/utils/adt/numeric.c,v 1.48 2001/11/05 17:46:29 momjian Exp $
8+
* $Header: /cvsroot/pgsql/src/backend/utils/adt/numeric.c,v 1.49 2001/12/11 02:02:12 tgl Exp $
99
*
1010
* ----------
1111
*/
@@ -159,6 +159,7 @@ static void add_var(NumericVar *var1, NumericVar *var2, NumericVar *result);
159159
static void sub_var(NumericVar *var1, NumericVar *var2, NumericVar *result);
160160
static void mul_var(NumericVar *var1, NumericVar *var2, NumericVar *result);
161161
static void div_var(NumericVar *var1, NumericVar *var2, NumericVar *result);
162+
static int select_div_scale(NumericVar *var1, NumericVar *var2);
162163
static void mod_var(NumericVar *var1, NumericVar *var2, NumericVar *result);
163164
static void ceil_var(NumericVar *var, NumericVar *result);
164165
static void floor_var(NumericVar *var, NumericVar *result);
@@ -999,28 +1000,7 @@ numeric_div(PG_FUNCTION_ARGS)
9991000
set_var_from_num(num1, &arg1);
10001001
set_var_from_num(num2, &arg2);
10011002

1002-
/* ----------
1003-
* The result scale of a division isn't specified in any
1004-
* SQL standard. For Postgres it is the following (where
1005-
* SR, DR are the result- and display-scales of the returned
1006-
* value, S1, D1, S2 and D2 are the scales of the two arguments,
1007-
* The minimum and maximum scales are compile time options from
1008-
* numeric.h):
1009-
*
1010-
* DR = MIN(MAX(D1 + D2, MIN_DISPLAY_SCALE), MAX_DISPLAY_SCALE)
1011-
* SR = MIN(MAX(MAX(S1 + S2, MIN_RESULT_SCALE), DR + 4), MAX_RESULT_SCALE)
1012-
*
1013-
* By default, any result is computed with a minimum of 34 digits
1014-
* after the decimal point or at least with 4 digits more than
1015-
* displayed.
1016-
* ----------
1017-
*/
1018-
res_dscale = MAX(arg1.dscale + arg2.dscale, NUMERIC_MIN_DISPLAY_SCALE);
1019-
res_dscale = MIN(res_dscale, NUMERIC_MAX_DISPLAY_SCALE);
1020-
global_rscale = MAX(arg1.rscale + arg2.rscale,
1021-
NUMERIC_MIN_RESULT_SCALE);
1022-
global_rscale = MAX(global_rscale, res_dscale + 4);
1023-
global_rscale = MIN(global_rscale, NUMERIC_MAX_RESULT_SCALE);
1003+
res_dscale = select_div_scale(&arg1, &arg2);
10241004

10251005
/*
10261006
* Do the divide, set the display scale and return the result
@@ -1884,6 +1864,7 @@ numeric_variance(PG_FUNCTION_ARGS)
18841864
vsumX,
18851865
vsumX2,
18861866
vNminus1;
1867+
int div_dscale;
18871868

18881869
/* We assume the input is array of numeric */
18891870
deconstruct_array(transarray,
@@ -1924,10 +1905,21 @@ numeric_variance(PG_FUNCTION_ARGS)
19241905
mul_var(&vsumX, &vsumX, &vsumX); /* now vsumX contains sumX * sumX */
19251906
mul_var(&vN, &vsumX2, &vsumX2); /* now vsumX2 contains N * sumX2 */
19261907
sub_var(&vsumX2, &vsumX, &vsumX2); /* N * sumX2 - sumX * sumX */
1927-
mul_var(&vN, &vNminus1, &vNminus1); /* N * (N - 1) */
1928-
div_var(&vsumX2, &vNminus1, &vsumX); /* variance */
19291908

1930-
res = make_result(&vsumX);
1909+
if (cmp_var(&vsumX2, &const_zero) <= 0)
1910+
{
1911+
/* Watch out for roundoff error producing a negative numerator */
1912+
res = make_result(&const_zero);
1913+
}
1914+
else
1915+
{
1916+
mul_var(&vN, &vNminus1, &vNminus1); /* N * (N - 1) */
1917+
div_dscale = select_div_scale(&vsumX2, &vNminus1);
1918+
div_var(&vsumX2, &vNminus1, &vsumX); /* variance */
1919+
vsumX.dscale = div_dscale;
1920+
1921+
res = make_result(&vsumX);
1922+
}
19311923

19321924
free_var(&vN);
19331925
free_var(&vNminus1);
@@ -1951,6 +1943,7 @@ numeric_stddev(PG_FUNCTION_ARGS)
19511943
vsumX,
19521944
vsumX2,
19531945
vNminus1;
1946+
int div_dscale;
19541947

19551948
/* We assume the input is array of numeric */
19561949
deconstruct_array(transarray,
@@ -1991,11 +1984,22 @@ numeric_stddev(PG_FUNCTION_ARGS)
19911984
mul_var(&vsumX, &vsumX, &vsumX); /* now vsumX contains sumX * sumX */
19921985
mul_var(&vN, &vsumX2, &vsumX2); /* now vsumX2 contains N * sumX2 */
19931986
sub_var(&vsumX2, &vsumX, &vsumX2); /* N * sumX2 - sumX * sumX */
1994-
mul_var(&vN, &vNminus1, &vNminus1); /* N * (N - 1) */
1995-
div_var(&vsumX2, &vNminus1, &vsumX); /* variance */
1996-
sqrt_var(&vsumX, &vsumX); /* stddev */
19971987

1998-
res = make_result(&vsumX);
1988+
if (cmp_var(&vsumX2, &const_zero) <= 0)
1989+
{
1990+
/* Watch out for roundoff error producing a negative numerator */
1991+
res = make_result(&const_zero);
1992+
}
1993+
else
1994+
{
1995+
mul_var(&vN, &vNminus1, &vNminus1); /* N * (N - 1) */
1996+
div_dscale = select_div_scale(&vsumX2, &vNminus1);
1997+
div_var(&vsumX2, &vNminus1, &vsumX); /* variance */
1998+
vsumX.dscale = div_dscale;
1999+
sqrt_var(&vsumX, &vsumX); /* stddev */
2000+
2001+
res = make_result(&vsumX);
2002+
}
19992003

20002004
free_var(&vN);
20012005
free_var(&vNminus1);
@@ -3318,6 +3322,50 @@ div_var(NumericVar *var1, NumericVar *var2, NumericVar *result)
33183322
}
33193323

33203324

3325+
/*
3326+
* Default scale selection for division
3327+
*
3328+
* Returns the appropriate display scale for the division result,
3329+
* and sets global_rscale to the result scale to use during div_var.
3330+
*
3331+
* Note that this must be called before div_var.
3332+
*/
3333+
static int
3334+
select_div_scale(NumericVar *var1, NumericVar *var2)
3335+
{
3336+
int res_dscale;
3337+
int res_rscale;
3338+
3339+
/* ----------
3340+
* The result scale of a division isn't specified in any
3341+
* SQL standard. For Postgres it is the following (where
3342+
* SR, DR are the result- and display-scales of the returned
3343+
* value, S1, D1, S2 and D2 are the scales of the two arguments,
3344+
* The minimum and maximum scales are compile time options from
3345+
* numeric.h):
3346+
*
3347+
* DR = MIN(MAX(D1 + D2, MIN_DISPLAY_SCALE), MAX_DISPLAY_SCALE)
3348+
* SR = MIN(MAX(MAX(S1 + S2, DR + 4), MIN_RESULT_SCALE), MAX_RESULT_SCALE)
3349+
*
3350+
* By default, any result is computed with a minimum of 34 digits
3351+
* after the decimal point or at least with 4 digits more than
3352+
* displayed.
3353+
* ----------
3354+
*/
3355+
res_dscale = var1->dscale + var2->dscale;
3356+
res_dscale = MAX(res_dscale, NUMERIC_MIN_DISPLAY_SCALE);
3357+
res_dscale = MIN(res_dscale, NUMERIC_MAX_DISPLAY_SCALE);
3358+
3359+
res_rscale = var1->rscale + var2->rscale;
3360+
res_rscale = MAX(res_rscale, res_dscale + 4);
3361+
res_rscale = MAX(res_rscale, NUMERIC_MIN_RESULT_SCALE);
3362+
res_rscale = MIN(res_rscale, NUMERIC_MAX_RESULT_SCALE);
3363+
global_rscale = res_rscale;
3364+
3365+
return res_dscale;
3366+
}
3367+
3368+
33213369
/* ----------
33223370
* mod_var() -
33233371
*
@@ -3343,12 +3391,7 @@ mod_var(NumericVar *var1, NumericVar *var2, NumericVar *result)
33433391
*/
33443392
save_global_rscale = global_rscale;
33453393

3346-
div_dscale = MAX(var1->dscale + var2->dscale, NUMERIC_MIN_DISPLAY_SCALE);
3347-
div_dscale = MIN(div_dscale, NUMERIC_MAX_DISPLAY_SCALE);
3348-
global_rscale = MAX(var1->rscale + var2->rscale,
3349-
NUMERIC_MIN_RESULT_SCALE);
3350-
global_rscale = MAX(global_rscale, div_dscale + 4);
3351-
global_rscale = MIN(global_rscale, NUMERIC_MAX_RESULT_SCALE);
3394+
div_dscale = select_div_scale(var1, var2);
33523395

33533396
div_var(var1, var2, &tmp);
33543397

0 commit comments

Comments
 (0)