diff --git a/src/backend/utils/adt/float.c b/src/backend/utils/adt/float.c
new file mode 100644
index df35557..ebb2d5c
--- a/src/backend/utils/adt/float.c
+++ b/src/backend/utils/adt/float.c
@@ -2401,13 +2401,29 @@ setseed(PG_FUNCTION_ARGS)
  *		float8_stddev_samp	- produce final result for float STDDEV_SAMP()
  *		float8_stddev_pop	- produce final result for float STDDEV_POP()
  *
+ * The naive schoolbook implementation of these aggregates works by
+ * accumulating sum(X) and sum(X^2).  However, this approach suffers from
+ * large rounding errors in the final computation of quantities like the
+ * population variance (N*sum(X^2) - sum(X)^2) / N^2, since each of the
+ * intermediate terms is potentially very large, while the difference is often
+ * quite small.
+ *
+ * Instead we use the Youngs-Cramer algorithm [1] which works by accumulating
+ * Sx=sum(X) and Sxx=sum((X-Sx/N)^2), using a numerically stable algorithm to
+ * incrementally update those quantities.  The final computations of each of
+ * the aggregate values is then trivial and gives more accurate results (for
+ * example, the population variance is just Sxx/N).
+ *
  * The transition datatype for all these aggregates is a 3-element array
- * of float8, holding the values N, sum(X), sum(X*X) in that order.
+ * of float8, holding the values N, Sx, Sxx in that order.
  *
  * Note that we represent N as a float to avoid having to build a special
  * datatype.  Given a reasonable floating-point implementation, there should
  * be no accuracy loss unless N exceeds 2 ^ 52 or so (by which time the
  * user will have doubtless lost interest anyway...)
+ *
+ * [1] Some Results Relevant to Choice of Sum and Sum-of-Product Algorithms,
+ * E. A. Youngs and E. M. Cramer, Technometrics Vol 13, No 3, August 1971.
  */
 
 static float8 *
@@ -2441,6 +2457,16 @@ float8_combine(PG_FUNCTION_ARGS)
 	ArrayType  *transarray2 = PG_GETARG_ARRAYTYPE_P(1);
 	float8	   *transvalues1;
 	float8	   *transvalues2;
+	float8		N1,
+				Sx1,
+				Sxx1,
+				N2,
+				Sx2,
+				Sxx2,
+				tmp,
+				N,
+				Sx,
+				Sxx;
 
 	if (!AggCheckCallContext(fcinfo, NULL))
 		elog(ERROR, "aggregate function called in non-aggregate context");
@@ -2448,9 +2474,51 @@ float8_combine(PG_FUNCTION_ARGS)
 	transvalues1 = check_float8_array(transarray1, "float8_combine", 3);
 	transvalues2 = check_float8_array(transarray2, "float8_combine", 3);
 
-	transvalues1[0] = transvalues1[0] + transvalues2[0];
-	transvalues1[1] = float8_pl(transvalues1[1], transvalues2[1]);
-	transvalues1[2] = float8_pl(transvalues1[2], transvalues2[2]);
+	N1 = transvalues1[0];
+	Sx1 = transvalues1[1];
+	Sxx1 = transvalues1[2];
+
+	N2 = transvalues2[0];
+	Sx2 = transvalues2[1];
+	Sxx2 = transvalues2[2];
+
+	/*--------------------
+	 * The transition values combine using a generalization of the
+	 * Youngs-Cramer algorithm as follows:
+	 *
+	 *	N = N1 + N2
+	 *	Sx = Sx1 + Sx2
+	 *	Sxx = Sxx1 + Sxx2 + N1 * N2 * (Sx1/N1 - Sx2/N2)^2 / N;
+	 *
+	 * It's worth handling the special cases N1 = 0 and N2 = 0 separately
+	 * since those cases are trivial, and we then don't need to worry about
+	 * division-by-zero errors in the general case.
+	 *--------------------
+	 */
+	if (N1 == 0.0)
+	{
+		N = N2;
+		Sx = Sx2;
+		Sxx = Sxx2;
+	}
+	else if (N2 == 0.0)
+	{
+		N = N1;
+		Sx = Sx1;
+		Sxx = Sxx1;
+	}
+	else
+	{
+		N = N1 + N2;
+		Sx = float8_pl(Sx1, Sx2);
+		tmp = Sx1 / N1 - Sx2 / N2;
+		Sxx = Sxx1 + Sxx2 + N1 * N2 * tmp * tmp / N;
+		check_float8_val(Sxx, isinf(Sxx1) || isinf(Sxx2), true);
+	}
+
+	transvalues1[0] = N;
+	transvalues1[1] = Sx;
+	transvalues1[2] = Sxx;
 
 	PG_RETURN_ARRAYTYPE_P(transarray1);
 }
@@ -2461,8 +2529,28 @@ float8_accum(PG_FUNCTION_ARGS)
 	ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
 	float8		newval = PG_GETARG_FLOAT8(1);
 	float8	   *transvalues;
+	float8		N,
+				Sx,
+				Sxx,
+				tmp;
 
 	transvalues = check_float8_array(transarray, "float8_accum", 3);
+	N = transvalues[0];
+	Sx = transvalues[1];
+	Sxx = transvalues[2];
+
+	/*
+	 * Use the Youngs-Cramer algorithm to incorporate the new value into the
+	 * transition values.
+	 */
+	N += 1.0;
+	Sx = float8_pl(Sx, newval);
+	if (transvalues[0] > 0.0)
+	{
+		tmp = newval * N - Sx;
+		Sxx += tmp * tmp / (N * transvalues[0]);
+		check_float8_val(Sxx, isinf(transvalues[2]) || isinf(newval), true);
+	}
 
 	/*
 	 * If we're invoked as an aggregate, we can cheat and modify our first
@@ -2471,9 +2559,9 @@ float8_accum(PG_FUNCTION_ARGS)
 	 */
 	if (AggCheckCallContext(fcinfo, NULL))
 	{
-		transvalues[0] = transvalues[0] + 1.0;
-		transvalues[1] = float8_pl(transvalues[1], newval);
-		transvalues[2] = float8_pl(transvalues[2], float8_mul(newval, newval));
+		transvalues[0] = N;
+		transvalues[1] = Sx;
+		transvalues[2] = Sxx;
 
 		PG_RETURN_ARRAYTYPE_P(transarray);
 	}
@@ -2482,9 +2570,9 @@ float8_accum(PG_FUNCTION_ARGS)
 		Datum		transdatums[3];
 		ArrayType  *result;
 
-		transvalues[0] = transvalues[0] + 1.0;
-		transvalues[1] = float8_pl(transvalues[1], newval);
-		transvalues[2] = float8_pl(transvalues[2], float8_mul(newval, newval));
+		transdatums[0] = Float8GetDatumFast(N);
+		transdatums[1] = Float8GetDatumFast(Sx);
+		transdatums[2] = Float8GetDatumFast(Sxx);
 
 		result = construct_array(transdatums, 3,
 								 FLOAT8OID,
@@ -2502,8 +2590,28 @@ float4_accum(PG_FUNCTION_ARGS)
 	/* do computations as float8 */
 	float8		newval = PG_GETARG_FLOAT4(1);
 	float8	   *transvalues;
+	float8		N,
+				Sx,
+				Sxx,
+				tmp;
 
 	transvalues = check_float8_array(transarray, "float4_accum", 3);
+	N = transvalues[0];
+	Sx = transvalues[1];
+	Sxx = transvalues[2];
+
+	/*
+	 * Use the Youngs-Cramer algorithm to incorporate the new value into the
+	 * transition values.
+	 */
+	N += 1.0;
+	Sx = float8_pl(Sx, newval);
+	if (transvalues[0] > 0.0)
+	{
+		tmp = newval * N - Sx;
+		Sxx += tmp * tmp / (N * transvalues[0]);
+		check_float8_val(Sxx, isinf(transvalues[2]) || isinf(newval), true);
+	}
 
 	/*
 	 * If we're invoked as an aggregate, we can cheat and modify our first
@@ -2512,9 +2620,9 @@ float4_accum(PG_FUNCTION_ARGS)
 	 */
 	if (AggCheckCallContext(fcinfo, NULL))
 	{
-		transvalues[0] = transvalues[0] + 1.0;
-		transvalues[1] = float8_pl(transvalues[1], newval);
-		transvalues[2] = float8_pl(transvalues[2], float8_mul(newval, newval));
+		transvalues[0] = N;
+		transvalues[1] = Sx;
+		transvalues[2] = Sxx;
 
 		PG_RETURN_ARRAYTYPE_P(transarray);
 	}
@@ -2523,9 +2631,9 @@ float4_accum(PG_FUNCTION_ARGS)
 		Datum		transdatums[3];
 		ArrayType  *result;
 
-		transvalues[0] = transvalues[0] + 1.0;
-		transvalues[1] = float8_pl(transvalues[1], newval);
-		transvalues[2] = float8_pl(transvalues[2], float8_mul(newval, newval));
+		transdatums[0] = Float8GetDatumFast(N);
+		transdatums[1] = Float8GetDatumFast(Sx);
+		transdatums[2] = Float8GetDatumFast(Sxx);
 
 		result = construct_array(transdatums, 3,
 								 FLOAT8OID,
@@ -2541,18 +2649,18 @@ float8_avg(PG_FUNCTION_ARGS)
 	ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
 	float8	   *transvalues;
 	float8		N,
-				sumX;
+				Sx;
 
 	transvalues = check_float8_array(transarray, "float8_avg", 3);
 	N = transvalues[0];
-	sumX = transvalues[1];
-	/* ignore sumX2 */
+	Sx = transvalues[1];
+	/* ignore Sxx */
 
 	/* SQL defines AVG of no values to be NULL */
 	if (N == 0.0)
 		PG_RETURN_NULL();
 
-	PG_RETURN_FLOAT8(sumX / N);
+	PG_RETURN_FLOAT8(Sx / N);
 }
 
 Datum
@@ -2561,27 +2669,22 @@ float8_var_pop(PG_FUNCTION_ARGS)
 	ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
 	float8	   *transvalues;
 	float8		N,
-				sumX,
-				sumX2,
-				numerator;
+				Sxx;
 
 	transvalues = check_float8_array(transarray, "float8_var_pop", 3);
 	N = transvalues[0];
-	sumX = transvalues[1];
-	sumX2 = transvalues[2];
+	/* ignore Sx */
+	Sxx = transvalues[2];
 
 	/* Population variance is undefined when N is 0, so return NULL */
 	if (N == 0.0)
 		PG_RETURN_NULL();
 
-	numerator = N * sumX2 - sumX * sumX;
-	check_float8_val(numerator, isinf(sumX2) || isinf(sumX), true);
-
-	/* Watch out for roundoff error producing a negative numerator */
-	if (numerator <= 0.0)
+	/* Watch out for roundoff error producing a negative variance */
+	if (Sxx <= 0.0)
 		PG_RETURN_FLOAT8(0.0);
 
-	PG_RETURN_FLOAT8(numerator / (N * N));
+	PG_RETURN_FLOAT8(Sxx / N);
 }
 
 Datum
@@ -2590,27 +2693,22 @@ float8_var_samp(PG_FUNCTION_ARGS)
 	ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
 	float8	   *transvalues;
 	float8		N,
-				sumX,
-				sumX2,
-				numerator;
+				Sxx;
 
 	transvalues = check_float8_array(transarray, "float8_var_samp", 3);
 	N = transvalues[0];
-	sumX = transvalues[1];
-	sumX2 = transvalues[2];
+	/* ignore Sx */
+	Sxx = transvalues[2];
 
 	/* Sample variance is undefined when N is 0 or 1, so return NULL */
 	if (N <= 1.0)
 		PG_RETURN_NULL();
 
-	numerator = N * sumX2 - sumX * sumX;
-	check_float8_val(numerator, isinf(sumX2) || isinf(sumX), true);
-
-	/* Watch out for roundoff error producing a negative numerator */
-	if (numerator <= 0.0)
+	/* Watch out for roundoff error producing a negative variance */
+	if (Sxx <= 0.0)
 		PG_RETURN_FLOAT8(0.0);
 
-	PG_RETURN_FLOAT8(numerator / (N * (N - 1.0)));
+	PG_RETURN_FLOAT8(Sxx / (N - 1.0));
 }
 
 Datum
@@ -2619,27 +2717,22 @@ float8_stddev_pop(PG_FUNCTION_ARGS)
 	ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
 	float8	   *transvalues;
 	float8		N,
-				sumX,
-				sumX2,
-				numerator;
+				Sxx;
 
 	transvalues = check_float8_array(transarray, "float8_stddev_pop", 3);
 	N = transvalues[0];
-	sumX = transvalues[1];
-	sumX2 = transvalues[2];
+	/* ignore Sx */
+	Sxx = transvalues[2];
 
 	/* Population stddev is undefined when N is 0, so return NULL */
 	if (N == 0.0)
 		PG_RETURN_NULL();
 
-	numerator = N * sumX2 - sumX * sumX;
-	check_float8_val(numerator, isinf(sumX2) || isinf(sumX), true);
-
-	/* Watch out for roundoff error producing a negative numerator */
-	if (numerator <= 0.0)
+	/* Watch out for roundoff error producing a negative variance */
+	if (Sxx <= 0.0)
 		PG_RETURN_FLOAT8(0.0);
 
-	PG_RETURN_FLOAT8(sqrt(numerator / (N * N)));
+	PG_RETURN_FLOAT8(sqrt(Sxx / N));
 }
 
 Datum
@@ -2648,27 +2741,22 @@ float8_stddev_samp(PG_FUNCTION_ARGS)
 	ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
 	float8	   *transvalues;
 	float8		N,
-				sumX,
-				sumX2,
-				numerator;
+				Sxx;
 
 	transvalues = check_float8_array(transarray, "float8_stddev_samp", 3);
 	N = transvalues[0];
-	sumX = transvalues[1];
-	sumX2 = transvalues[2];
+	/* ignore Sx */
+	Sxx = transvalues[2];
 
 	/* Sample stddev is undefined when N is 0 or 1, so return NULL */
 	if (N <= 1.0)
 		PG_RETURN_NULL();
 
-	numerator = N * sumX2 - sumX * sumX;
-	check_float8_val(numerator, isinf(sumX2) || isinf(sumX), true);
-
-	/* Watch out for roundoff error producing a negative numerator */
-	if (numerator <= 0.0)
+	/* Watch out for roundoff error producing a negative variance */
+	if (Sxx <= 0.0)
 		PG_RETURN_FLOAT8(0.0);
 
-	PG_RETURN_FLOAT8(sqrt(numerator / (N * (N - 1.0))));
+	PG_RETURN_FLOAT8(sqrt(Sxx / (N - 1.0)));
 }
 
 /*
@@ -2676,9 +2764,14 @@ float8_stddev_samp(PG_FUNCTION_ARGS)
  *		SQL2003 BINARY AGGREGATES
  *		=========================
  *
+ * As with the preceding aggregates, we use the Youngs-Cramer algorithm to
+ * reduce rounding errors in the aggregate final functions.
+ *
  * The transition datatype for all these aggregates is a 6-element array of
- * float8, holding the values N, sum(X), sum(X*X), sum(Y), sum(Y*Y), sum(X*Y)
- * in that order.  Note that Y is the first argument to the aggregates!
+ * float8, holding the values N, Sx=sum(X), Sxx=sum((X-Sx/N)^2), Sy=sum(Y),
+ * Syy=sum((Y-Sy/N)^2), Sxy=sum((X-Sx/N)*(Y-Sy/N)) in that order.
+ *
+ * Note that Y is the first argument to all these aggregates!
  *
  * It might seem attractive to optimize this by having multiple accumulator
  * functions that only calculate the sums actually needed.  But on most
@@ -2695,32 +2788,43 @@ float8_regr_accum(PG_FUNCTION_ARGS)
 	float8		newvalX = PG_GETARG_FLOAT8(2);
 	float8	   *transvalues;
 	float8		N,
-				sumX,
-				sumX2,
-				sumY,
-				sumY2,
-				sumXY;
+				Sx,
+				Sxx,
+				Sy,
+				Syy,
+				Sxy,
+				factor,
+				tmp1,
+				tmp2;
 
 	transvalues = check_float8_array(transarray, "float8_regr_accum", 6);
 	N = transvalues[0];
-	sumX = transvalues[1];
-	sumX2 = transvalues[2];
-	sumY = transvalues[3];
-	sumY2 = transvalues[4];
-	sumXY = transvalues[5];
+	Sx = transvalues[1];
+	Sxx = transvalues[2];
+	Sy = transvalues[3];
+	Syy = transvalues[4];
+	Sxy = transvalues[5];
 
+	/*
+	 * Use the Youngs-Cramer algorithm to incorporate the new values into the
+	 * transition values.
+	 */
 	N += 1.0;
-	sumX += newvalX;
-	check_float8_val(sumX, isinf(transvalues[1]) || isinf(newvalX), true);
-	sumX2 += newvalX * newvalX;
-	check_float8_val(sumX2, isinf(transvalues[2]) || isinf(newvalX), true);
-	sumY += newvalY;
-	check_float8_val(sumY, isinf(transvalues[3]) || isinf(newvalY), true);
-	sumY2 += newvalY * newvalY;
-	check_float8_val(sumY2, isinf(transvalues[4]) || isinf(newvalY), true);
-	sumXY += newvalX * newvalY;
-	check_float8_val(sumXY, isinf(transvalues[5]) || isinf(newvalX) ||
-					 isinf(newvalY), true);
+	Sx = float8_pl(Sx, newvalX);
+	Sy = float8_pl(Sy, newvalY);
+	if (transvalues[0] > 0.0)
+	{
+		factor = 1.0 / (N * transvalues[0]);
+		tmp1 = newvalX * N - Sx;
+		Sxx += tmp1 * tmp1 * factor;
+		check_float8_val(Sxx, isinf(transvalues[2]) || isinf(newvalX), true);
+		tmp2 = newvalY * N - Sy;
+		Syy += tmp2 * tmp2 * factor;
+		check_float8_val(Syy, isinf(transvalues[4]) || isinf(newvalY), true);
+		Sxy += tmp1 * tmp2 * factor;
+		check_float8_val(Sxy, isinf(transvalues[5]) || isinf(newvalX) ||
+						 isinf(newvalY), true);
+	}
 
 	/*
 	 * If we're invoked as an aggregate, we can cheat and modify our first
@@ -2730,11 +2834,11 @@ float8_regr_accum(PG_FUNCTION_ARGS)
 	if (AggCheckCallContext(fcinfo, NULL))
 	{
 		transvalues[0] = N;
-		transvalues[1] = sumX;
-		transvalues[2] = sumX2;
-		transvalues[3] = sumY;
-		transvalues[4] = sumY2;
-		transvalues[5] = sumXY;
+		transvalues[1] = Sx;
+		transvalues[2] = Sxx;
+		transvalues[3] = Sy;
+		transvalues[4] = Syy;
+		transvalues[5] = Sxy;
 
 		PG_RETURN_ARRAYTYPE_P(transarray);
 	}
@@ -2744,11 +2848,11 @@ float8_regr_accum(PG_FUNCTION_ARGS)
 		ArrayType  *result;
 
 		transdatums[0] = Float8GetDatumFast(N);
-		transdatums[1] = Float8GetDatumFast(sumX);
-		transdatums[2] = Float8GetDatumFast(sumX2);
-		transdatums[3] = Float8GetDatumFast(sumY);
-		transdatums[4] = Float8GetDatumFast(sumY2);
-		transdatums[5] = Float8GetDatumFast(sumXY);
+		transdatums[1] = Float8GetDatumFast(Sx);
+		transdatums[2] = Float8GetDatumFast(Sxx);
+		transdatums[3] = Float8GetDatumFast(Sy);
+		transdatums[4] = Float8GetDatumFast(Syy);
+		transdatums[5] = Float8GetDatumFast(Sxy);
 
 		result = construct_array(transdatums, 6,
 								 FLOAT8OID,
@@ -2773,6 +2877,26 @@ float8_regr_combine(PG_FUNCTION_ARGS)
 	ArrayType  *transarray2 = PG_GETARG_ARRAYTYPE_P(1);
 	float8	   *transvalues1;
 	float8	   *transvalues2;
+	float8		N1,
+				Sx1,
+				Sxx1,
+				Sy1,
+				Syy1,
+				Sxy1,
+				N2,
+				Sx2,
+				Sxx2,
+				Sy2,
+				Syy2,
+				Sxy2,
+				tmp1,
+				tmp2,
+				N,
+				Sx,
+				Sxx,
+				Sy,
+				Syy,
+				Sxy;
 
 	if (!AggCheckCallContext(fcinfo, NULL))
 		elog(ERROR, "aggregate function called in non-aggregate context");
@@ -2780,12 +2904,75 @@ float8_regr_combine(PG_FUNCTION_ARGS)
 	transvalues1 = check_float8_array(transarray1, "float8_regr_combine", 6);
 	transvalues2 = check_float8_array(transarray2, "float8_regr_combine", 6);
 
-	transvalues1[0] = transvalues1[0] + transvalues2[0];
-	transvalues1[1] = float8_pl(transvalues1[1], transvalues2[1]);
-	transvalues1[2] = float8_pl(transvalues1[2], transvalues2[2]);
-	transvalues1[3] = float8_pl(transvalues1[3], transvalues2[3]);
-	transvalues1[4] = float8_pl(transvalues1[4], transvalues2[4]);
-	transvalues1[5] = float8_pl(transvalues1[5], transvalues2[5]);
+	N1 = transvalues1[0];
+	Sx1 = transvalues1[1];
+	Sxx1 = transvalues1[2];
+	Sy1 = transvalues1[3];
+	Syy1 = transvalues1[4];
+	Sxy1 = transvalues1[5];
+
+	N2 = transvalues2[0];
+	Sx2 = transvalues2[1];
+	Sxx2 = transvalues2[2];
+	Sy2 = transvalues2[3];
+	Syy2 = transvalues2[4];
+	Sxy2 = transvalues2[5];
+
+	/*--------------------
+	 * The transition values combine using a generalization of the
+	 * Youngs-Cramer algorithm as follows:
+	 *
+	 *	N = N1 + N2
+	 *	Sx = Sx1 + Sx2
+	 *	Sxx = Sxx1 + Sxx2 + N1 * N2 * (Sx1/N1 - Sx2/N2)^2 / N
+	 *	Sy = Sy1 + Sy2
+	 *	Syy = Syy1 + Syy2 + N1 * N2 * (Sy1/N1 - Sy2/N2)^2 / N
+	 *	Sxy = Sxy1 + Sxy2 + N1 * N2 * (Sx1/N1 - Sx2/N2) * (Sy1/N1 - Sy2/N2) / N
+	 *
+	 * It's worth handling the special cases N1 = 0 and N2 = 0 separately
+	 * since those cases are trivial, and we then don't need to worry about
+	 * division-by-zero errors in the general case.
+	 *--------------------
+	 */
+	if (N1 == 0.0)
+	{
+		N = N2;
+		Sx = Sx2;
+		Sxx = Sxx2;
+		Sy = Sy2;
+		Syy = Syy2;
+		Sxy = Sxy2;
+	}
+	else if (N2 == 0.0)
+	{
+		N = N1;
+		Sx = Sx1;
+		Sxx = Sxx1;
+		Sy = Syy1;
+		Syy = Syy1;
+		Sxy = Sxy1;
+	}
+	else
+	{
+		N = N1 + N2;
+		Sx = float8_pl(Sx1, Sx2);
+		tmp1 = Sx1 / N1 - Sx2 / N2;
+		Sxx = Sxx1 + Sxx2 + N1 * N2 * tmp1 * tmp1 / N;
+		check_float8_val(Sxx, isinf(Sxx1) || isinf(Sxx2), true);
+		Sy = float8_pl(Sy1, Sy2);
+		tmp2 = Sy1 / N1 - Sy2 / N2;
+		Syy = Syy1 + Syy2 + N1 * N2 * tmp2 * tmp2 / N;
+		check_float8_val(Syy, isinf(Syy1) || isinf(Syy2), true);
+		Sxy = Sxy1 + Sxy2 + N1 * N2 * tmp1 * tmp2 / N;
+		check_float8_val(Sxy, isinf(Sxy1) || isinf(Sxy2), true);
+	}
+
+	transvalues1[0] = N;
+	transvalues1[1] = Sx;
+	transvalues1[2] = Sxx;
+	transvalues1[3] = Sy;
+	transvalues1[4] = Syy;
+	transvalues1[5] = Sxy;
 
 	PG_RETURN_ARRAYTYPE_P(transarray1);
 }
@@ -2797,27 +2984,21 @@ float8_regr_sxx(PG_FUNCTION_ARGS)
 	ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
 	float8	   *transvalues;
 	float8		N,
-				sumX,
-				sumX2,
-				numerator;
+				Sxx;
 
 	transvalues = check_float8_array(transarray, "float8_regr_sxx", 6);
 	N = transvalues[0];
-	sumX = transvalues[1];
-	sumX2 = transvalues[2];
+	Sxx = transvalues[2];
 
 	/* if N is 0 we should return NULL */
 	if (N < 1.0)
 		PG_RETURN_NULL();
 
-	numerator = N * sumX2 - sumX * sumX;
-	check_float8_val(numerator, isinf(sumX2) || isinf(sumX), true);
-
-	/* Watch out for roundoff error producing a negative numerator */
-	if (numerator <= 0.0)
+	/* Watch out for roundoff error producing a negative result */
+	if (Sxx <= 0.0)
 		PG_RETURN_FLOAT8(0.0);
 
-	PG_RETURN_FLOAT8(numerator / N);
+	PG_RETURN_FLOAT8(Sxx);
 }
 
 Datum
@@ -2826,27 +3007,21 @@ float8_regr_syy(PG_FUNCTION_ARGS)
 	ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
 	float8	   *transvalues;
 	float8		N,
-				sumY,
-				sumY2,
-				numerator;
+				Syy;
 
 	transvalues = check_float8_array(transarray, "float8_regr_syy", 6);
 	N = transvalues[0];
-	sumY = transvalues[3];
-	sumY2 = transvalues[4];
+	Syy = transvalues[4];
 
 	/* if N is 0 we should return NULL */
 	if (N < 1.0)
 		PG_RETURN_NULL();
 
-	numerator = N * sumY2 - sumY * sumY;
-	check_float8_val(numerator, isinf(sumY2) || isinf(sumY), true);
-
-	/* Watch out for roundoff error producing a negative numerator */
-	if (numerator <= 0.0)
+	/* Watch out for roundoff error producing a negative result */
+	if (Syy <= 0.0)
 		PG_RETURN_FLOAT8(0.0);
 
-	PG_RETURN_FLOAT8(numerator / N);
+	PG_RETURN_FLOAT8(Syy);
 }
 
 Datum
@@ -2855,28 +3030,19 @@ float8_regr_sxy(PG_FUNCTION_ARGS)
 	ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
 	float8	   *transvalues;
 	float8		N,
-				sumX,
-				sumY,
-				sumXY,
-				numerator;
+				Sxy;
 
 	transvalues = check_float8_array(transarray, "float8_regr_sxy", 6);
 	N = transvalues[0];
-	sumX = transvalues[1];
-	sumY = transvalues[3];
-	sumXY = transvalues[5];
+	Sxy = transvalues[5];
 
 	/* if N is 0 we should return NULL */
 	if (N < 1.0)
 		PG_RETURN_NULL();
 
-	numerator = N * sumXY - sumX * sumY;
-	check_float8_val(numerator, isinf(sumXY) || isinf(sumX) ||
-					 isinf(sumY), true);
-
 	/* A negative result is valid here */
 
-	PG_RETURN_FLOAT8(numerator / N);
+	PG_RETURN_FLOAT8(Sxy);
 }
 
 Datum
@@ -2885,17 +3051,17 @@ float8_regr_avgx(PG_FUNCTION_ARGS)
 	ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
 	float8	   *transvalues;
 	float8		N,
-				sumX;
+				Sx;
 
 	transvalues = check_float8_array(transarray, "float8_regr_avgx", 6);
 	N = transvalues[0];
-	sumX = transvalues[1];
+	Sx = transvalues[1];
 
 	/* if N is 0 we should return NULL */
 	if (N < 1.0)
 		PG_RETURN_NULL();
 
-	PG_RETURN_FLOAT8(sumX / N);
+	PG_RETURN_FLOAT8(Sx / N);
 }
 
 Datum
@@ -2904,17 +3070,17 @@ float8_regr_avgy(PG_FUNCTION_ARGS)
 	ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
 	float8	   *transvalues;
 	float8		N,
-				sumY;
+				Sy;
 
 	transvalues = check_float8_array(transarray, "float8_regr_avgy", 6);
 	N = transvalues[0];
-	sumY = transvalues[3];
+	Sy = transvalues[3];
 
 	/* if N is 0 we should return NULL */
 	if (N < 1.0)
 		PG_RETURN_NULL();
 
-	PG_RETURN_FLOAT8(sumY / N);
+	PG_RETURN_FLOAT8(Sy / N);
 }
 
 Datum
@@ -2923,26 +3089,17 @@ float8_covar_pop(PG_FUNCTION_ARGS)
 	ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
 	float8	   *transvalues;
 	float8		N,
-				sumX,
-				sumY,
-				sumXY,
-				numerator;
+				Sxy;
 
 	transvalues = check_float8_array(transarray, "float8_covar_pop", 6);
 	N = transvalues[0];
-	sumX = transvalues[1];
-	sumY = transvalues[3];
-	sumXY = transvalues[5];
+	Sxy = transvalues[5];
 
 	/* if N is 0 we should return NULL */
 	if (N < 1.0)
 		PG_RETURN_NULL();
 
-	numerator = N * sumXY - sumX * sumY;
-	check_float8_val(numerator, isinf(sumXY) || isinf(sumX) ||
-					 isinf(sumY), true);
-
-	PG_RETURN_FLOAT8(numerator / (N * N));
+	PG_RETURN_FLOAT8(Sxy / N);
 }
 
 Datum
@@ -2951,26 +3108,17 @@ float8_covar_samp(PG_FUNCTION_ARGS)
 	ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
 	float8	   *transvalues;
 	float8		N,
-				sumX,
-				sumY,
-				sumXY,
-				numerator;
+				Sxy;
 
 	transvalues = check_float8_array(transarray, "float8_covar_samp", 6);
 	N = transvalues[0];
-	sumX = transvalues[1];
-	sumY = transvalues[3];
-	sumXY = transvalues[5];
+	Sxy = transvalues[5];
 
 	/* if N is <= 1 we should return NULL */
 	if (N < 2.0)
 		PG_RETURN_NULL();
 
-	numerator = N * sumXY - sumX * sumY;
-	check_float8_val(numerator, isinf(sumXY) || isinf(sumX) ||
-					 isinf(sumY), true);
-
-	PG_RETURN_FLOAT8(numerator / (N * (N - 1.0)));
+	PG_RETURN_FLOAT8(Sxy / (N - 1.0));
 }
 
 Datum
@@ -2979,38 +3127,24 @@ float8_corr(PG_FUNCTION_ARGS)
 	ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
 	float8	   *transvalues;
 	float8		N,
-				sumX,
-				sumX2,
-				sumY,
-				sumY2,
-				sumXY,
-				numeratorX,
-				numeratorY,
-				numeratorXY;
+				Sxx,
+				Syy,
+				Sxy;
 
 	transvalues = check_float8_array(transarray, "float8_corr", 6);
 	N = transvalues[0];
-	sumX = transvalues[1];
-	sumX2 = transvalues[2];
-	sumY = transvalues[3];
-	sumY2 = transvalues[4];
-	sumXY = transvalues[5];
+	Sxx = transvalues[2];
+	Syy = transvalues[4];
+	Sxy = transvalues[5];
 
 	/* if N is 0 we should return NULL */
 	if (N < 1.0)
 		PG_RETURN_NULL();
 
-	numeratorX = N * sumX2 - sumX * sumX;
-	check_float8_val(numeratorX, isinf(sumX2) || isinf(sumX), true);
-	numeratorY = N * sumY2 - sumY * sumY;
-	check_float8_val(numeratorY, isinf(sumY2) || isinf(sumY), true);
-	numeratorXY = N * sumXY - sumX * sumY;
-	check_float8_val(numeratorXY, isinf(sumXY) || isinf(sumX) ||
-					 isinf(sumY), true);
-	if (numeratorX <= 0 || numeratorY <= 0)
+	if (Sxx <= 0 || Syy <= 0)
 		PG_RETURN_NULL();
 
-	PG_RETURN_FLOAT8(numeratorXY / sqrt(numeratorX * numeratorY));
+	PG_RETURN_FLOAT8(Sxy / sqrt(Sxx * Syy));
 }
 
 Datum
@@ -3019,42 +3153,27 @@ float8_regr_r2(PG_FUNCTION_ARGS)
 	ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
 	float8	   *transvalues;
 	float8		N,
-				sumX,
-				sumX2,
-				sumY,
-				sumY2,
-				sumXY,
-				numeratorX,
-				numeratorY,
-				numeratorXY;
+				Sxx,
+				Syy,
+				Sxy;
 
 	transvalues = check_float8_array(transarray, "float8_regr_r2", 6);
 	N = transvalues[0];
-	sumX = transvalues[1];
-	sumX2 = transvalues[2];
-	sumY = transvalues[3];
-	sumY2 = transvalues[4];
-	sumXY = transvalues[5];
+	Sxx = transvalues[2];
+	Syy = transvalues[4];
+	Sxy = transvalues[5];
 
 	/* if N is 0 we should return NULL */
 	if (N < 1.0)
 		PG_RETURN_NULL();
 
-	numeratorX = N * sumX2 - sumX * sumX;
-	check_float8_val(numeratorX, isinf(sumX2) || isinf(sumX), true);
-	numeratorY = N * sumY2 - sumY * sumY;
-	check_float8_val(numeratorY, isinf(sumY2) || isinf(sumY), true);
-	numeratorXY = N * sumXY - sumX * sumY;
-	check_float8_val(numeratorXY, isinf(sumXY) || isinf(sumX) ||
-					 isinf(sumY), true);
-	if (numeratorX <= 0)
+	if (Sxx <= 0)
 		PG_RETURN_NULL();
 	/* per spec, horizontal line produces 1.0 */
-	if (numeratorY <= 0)
+	if (Syy <= 0)
 		PG_RETURN_FLOAT8(1.0);
 
-	PG_RETURN_FLOAT8((numeratorXY * numeratorXY) /
-					 (numeratorX * numeratorY));
+	PG_RETURN_FLOAT8((Sxy * Sxy) / (Sxx * Syy));
 }
 
 Datum
@@ -3063,33 +3182,22 @@ float8_regr_slope(PG_FUNCTION_ARGS)
 	ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
 	float8	   *transvalues;
 	float8		N,
-				sumX,
-				sumX2,
-				sumY,
-				sumXY,
-				numeratorX,
-				numeratorXY;
+				Sxx,
+				Sxy;
 
 	transvalues = check_float8_array(transarray, "float8_regr_slope", 6);
 	N = transvalues[0];
-	sumX = transvalues[1];
-	sumX2 = transvalues[2];
-	sumY = transvalues[3];
-	sumXY = transvalues[5];
+	Sxx = transvalues[2];
+	Sxy = transvalues[5];
 
 	/* if N is 0 we should return NULL */
 	if (N < 1.0)
 		PG_RETURN_NULL();
 
-	numeratorX = N * sumX2 - sumX * sumX;
-	check_float8_val(numeratorX, isinf(sumX2) || isinf(sumX), true);
-	numeratorXY = N * sumXY - sumX * sumY;
-	check_float8_val(numeratorXY, isinf(sumXY) || isinf(sumX) ||
-					 isinf(sumY), true);
-	if (numeratorX <= 0)
+	if (Sxx <= 0)
 		PG_RETURN_NULL();
 
-	PG_RETURN_FLOAT8(numeratorXY / numeratorX);
+	PG_RETURN_FLOAT8(Sxy / Sxx);
 }
 
 Datum
@@ -3098,33 +3206,26 @@ float8_regr_intercept(PG_FUNCTION_ARGS)
 	ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
 	float8	   *transvalues;
 	float8		N,
-				sumX,
-				sumX2,
-				sumY,
-				sumXY,
-				numeratorX,
-				numeratorXXY;
+				Sx,
+				Sxx,
+				Sy,
+				Sxy;
 
 	transvalues = check_float8_array(transarray, "float8_regr_intercept", 6);
 	N = transvalues[0];
-	sumX = transvalues[1];
-	sumX2 = transvalues[2];
-	sumY = transvalues[3];
-	sumXY = transvalues[5];
+	Sx = transvalues[1];
+	Sxx = transvalues[2];
+	Sy = transvalues[3];
+	Sxy = transvalues[5];
 
 	/* if N is 0 we should return NULL */
 	if (N < 1.0)
 		PG_RETURN_NULL();
 
-	numeratorX = N * sumX2 - sumX * sumX;
-	check_float8_val(numeratorX, isinf(sumX2) || isinf(sumX), true);
-	numeratorXXY = sumY * sumX2 - sumX * sumXY;
-	check_float8_val(numeratorXXY, isinf(sumY) || isinf(sumX2) ||
-					 isinf(sumX) || isinf(sumXY), true);
-	if (numeratorX <= 0)
+	if (Sxx <= 0)
 		PG_RETURN_NULL();
 
-	PG_RETURN_FLOAT8(numeratorXXY / numeratorX);
+	PG_RETURN_FLOAT8((Sy - Sx * Sxy / Sxx) / N);
 }