diff --git a/src/backend/utils/adt/numeric.c b/src/backend/utils/adt/numeric.c
index 4ea0fec1c1a13e71ae41bfbffc06caffe2b765d3..7f32e2e62438c2aed5da16ffb6e29379e13c7f6f 100644
--- a/src/backend/utils/adt/numeric.c
+++ b/src/backend/utils/adt/numeric.c
@@ -5,7 +5,7 @@
  *
  *	1998 Jan Wieck
  *
- * $Header: /cvsroot/pgsql/src/backend/utils/adt/numeric.c,v 1.54 2002/09/18 21:35:22 tgl Exp $
+ * $Header: /cvsroot/pgsql/src/backend/utils/adt/numeric.c,v 1.55 2002/10/02 19:21:26 tgl Exp $
  *
  * ----------
  */
@@ -88,7 +88,7 @@ typedef struct NumericVar
  * Local data
  * ----------
  */
-static int	global_rscale = NUMERIC_MIN_RESULT_SCALE;
+static int	global_rscale = 0;
 
 /* ----------
  * Some preinitialized variables we need often
@@ -106,6 +106,18 @@ static NumericDigit const_two_data[1] = {2};
 static NumericVar const_two =
 {1, 0, 0, 0, NUMERIC_POS, NULL, const_two_data};
 
+static NumericDigit const_zero_point_one_data[1] = {1};
+static NumericVar const_zero_point_one =
+{1, -1, 1, 1, NUMERIC_POS, NULL, const_zero_point_one_data};
+
+static NumericDigit const_zero_point_nine_data[1] = {9};
+static NumericVar const_zero_point_nine =
+{1, -1, 1, 1, NUMERIC_POS, NULL, const_zero_point_nine_data};
+
+static NumericDigit const_one_point_one_data[2] = {1, 1};
+static NumericVar const_one_point_one =
+{2, 0, 1, 1, NUMERIC_POS, NULL, const_one_point_one_data};
+
 static NumericVar const_nan =
 {0, 0, 0, 0, NUMERIC_NAN, NULL, NULL};
 
@@ -146,6 +158,9 @@ static Numeric make_result(NumericVar *var);
 
 static void apply_typmod(NumericVar *var, int32 typmod);
 
+static double numeric_to_double_no_overflow(Numeric num);
+static double numericvar_to_double_no_overflow(NumericVar *var);
+
 static int	cmp_numerics(Numeric num1, Numeric num2);
 static int	cmp_var(NumericVar *var1, NumericVar *var2);
 static void add_var(NumericVar *var1, NumericVar *var2, NumericVar *result);
@@ -477,8 +492,8 @@ numeric_round(PG_FUNCTION_ARGS)
 	 * Limit the scale value to avoid possible overflow in calculations
 	 * below.
 	 */
-	scale = Min(NUMERIC_MAX_RESULT_SCALE,
-				Max(-NUMERIC_MAX_RESULT_SCALE, scale));
+	scale = Max(scale, -NUMERIC_MAX_RESULT_SCALE);
+	scale = Min(scale, NUMERIC_MAX_RESULT_SCALE);
 
 	/*
 	 * Unpack the argument and round it at the proper digit position
@@ -563,8 +578,8 @@ numeric_trunc(PG_FUNCTION_ARGS)
 	 * Limit the scale value to avoid possible overflow in calculations
 	 * below.
 	 */
-	scale = Min(NUMERIC_MAX_RESULT_SCALE,
-				Max(-NUMERIC_MAX_RESULT_SCALE, scale));
+	scale = Max(scale, -NUMERIC_MAX_RESULT_SCALE);
+	scale = Min(scale, NUMERIC_MAX_RESULT_SCALE);
 
 	/*
 	 * Unpack the argument and truncate it at the proper digit position
@@ -1190,6 +1205,7 @@ numeric_sqrt(PG_FUNCTION_ARGS)
 	Numeric		res;
 	NumericVar	arg;
 	NumericVar	result;
+	int			sweight;
 	int			res_dscale;
 
 	/*
@@ -1199,20 +1215,28 @@ numeric_sqrt(PG_FUNCTION_ARGS)
 		PG_RETURN_NUMERIC(make_result(&const_nan));
 
 	/*
-	 * Unpack the argument, determine the scales like for divide, let
-	 * sqrt_var() do the calculation and return the result.
+	 * Unpack the argument and determine the scales.  We choose a display
+	 * scale to give at least NUMERIC_MIN_SIG_DIGITS significant digits;
+	 * but in any case not less than the input's dscale.
 	 */
 	init_var(&arg);
 	init_var(&result);
 
 	set_var_from_num(num, &arg);
 
-	res_dscale = Max(arg.dscale, NUMERIC_MIN_DISPLAY_SCALE);
+	/* Assume the input was normalized, so arg.weight is accurate */
+	sweight = (arg.weight / 2) - 1;
+
+	res_dscale = NUMERIC_MIN_SIG_DIGITS - sweight;
+	res_dscale = Max(res_dscale, arg.dscale);
+	res_dscale = Max(res_dscale, NUMERIC_MIN_DISPLAY_SCALE);
 	res_dscale = Min(res_dscale, NUMERIC_MAX_DISPLAY_SCALE);
-	global_rscale = Max(arg.rscale, NUMERIC_MIN_RESULT_SCALE);
-	global_rscale = Max(global_rscale, res_dscale + 4);
-	global_rscale = Min(global_rscale, NUMERIC_MAX_RESULT_SCALE);
 
+	global_rscale = res_dscale + 8;
+
+	/*
+	 * Let sqrt_var() do the calculation and return the result.
+	 */
 	sqrt_var(&arg, &result);
 
 	result.dscale = res_dscale;
@@ -1240,6 +1264,7 @@ numeric_exp(PG_FUNCTION_ARGS)
 	NumericVar	arg;
 	NumericVar	result;
 	int			res_dscale;
+	double		val;
 
 	/*
 	 * Handle NaN
@@ -1248,18 +1273,35 @@ numeric_exp(PG_FUNCTION_ARGS)
 		PG_RETURN_NUMERIC(make_result(&const_nan));
 
 	/*
-	 * Same procedure like for sqrt().
+	 * Unpack the argument and determine the scales.  We choose a display
+	 * scale to give at least NUMERIC_MIN_SIG_DIGITS significant digits;
+	 * but in any case not less than the input's dscale.
 	 */
 	init_var(&arg);
 	init_var(&result);
+
 	set_var_from_num(num, &arg);
 
-	res_dscale = Max(arg.dscale, NUMERIC_MIN_DISPLAY_SCALE);
+	/* convert input to float8, ignoring overflow */
+	val = numeric_to_double_no_overflow(num);
+
+	/* log10(result) = num * log10(e), so this is approximately the weight: */
+	val *= 0.434294481903252;
+
+	/* limit to something that won't cause integer overflow */
+	val = Max(val, -NUMERIC_MAX_RESULT_SCALE);
+	val = Min(val, NUMERIC_MAX_RESULT_SCALE);
+
+	res_dscale = NUMERIC_MIN_SIG_DIGITS - (int) val;
+	res_dscale = Max(res_dscale, arg.dscale);
+	res_dscale = Max(res_dscale, NUMERIC_MIN_DISPLAY_SCALE);
 	res_dscale = Min(res_dscale, NUMERIC_MAX_DISPLAY_SCALE);
-	global_rscale = Max(arg.rscale, NUMERIC_MIN_RESULT_SCALE);
-	global_rscale = Max(global_rscale, res_dscale + 4);
-	global_rscale = Min(global_rscale, NUMERIC_MAX_RESULT_SCALE);
 
+	global_rscale = res_dscale + NUMERIC_EXTRA_DIGITS;
+
+	/*
+	 * Let exp_var() do the calculation and return the result.
+	 */
 	exp_var(&arg, &result);
 
 	result.dscale = res_dscale;
@@ -1294,18 +1336,23 @@ numeric_ln(PG_FUNCTION_ARGS)
 	if (NUMERIC_IS_NAN(num))
 		PG_RETURN_NUMERIC(make_result(&const_nan));
 
-	/*
-	 * Same procedure like for sqrt()
-	 */
 	init_var(&arg);
 	init_var(&result);
+
 	set_var_from_num(num, &arg);
 
-	res_dscale = Max(arg.dscale, NUMERIC_MIN_DISPLAY_SCALE);
+	if (arg.weight > 0)
+		res_dscale = NUMERIC_MIN_SIG_DIGITS - (int) log10(arg.weight);
+	else if (arg.weight < 0)
+		res_dscale = NUMERIC_MIN_SIG_DIGITS - (int) log10(- arg.weight);
+	else
+		res_dscale = NUMERIC_MIN_SIG_DIGITS;
+
+	res_dscale = Max(res_dscale, arg.dscale);
+	res_dscale = Max(res_dscale, NUMERIC_MIN_DISPLAY_SCALE);
 	res_dscale = Min(res_dscale, NUMERIC_MAX_DISPLAY_SCALE);
-	global_rscale = Max(arg.rscale, NUMERIC_MIN_RESULT_SCALE);
-	global_rscale = Max(global_rscale, res_dscale + 4);
-	global_rscale = Min(global_rscale, NUMERIC_MAX_RESULT_SCALE);
+
+	global_rscale = res_dscale + NUMERIC_EXTRA_DIGITS;
 
 	ln_var(&arg, &result);
 
@@ -1335,7 +1382,6 @@ numeric_log(PG_FUNCTION_ARGS)
 	NumericVar	arg1;
 	NumericVar	arg2;
 	NumericVar	result;
-	int			res_dscale;
 
 	/*
 	 * Handle NaN
@@ -1344,27 +1390,21 @@ numeric_log(PG_FUNCTION_ARGS)
 		PG_RETURN_NUMERIC(make_result(&const_nan));
 
 	/*
-	 * Initialize things and calculate scales
+	 * Initialize things
 	 */
 	init_var(&arg1);
 	init_var(&arg2);
 	init_var(&result);
+
 	set_var_from_num(num1, &arg1);
 	set_var_from_num(num2, &arg2);
 
-	res_dscale = Max(arg1.dscale + arg2.dscale, NUMERIC_MIN_DISPLAY_SCALE);
-	res_dscale = Min(res_dscale, NUMERIC_MAX_DISPLAY_SCALE);
-	global_rscale = Max(arg1.rscale + arg2.rscale, NUMERIC_MIN_RESULT_SCALE);
-	global_rscale = Max(global_rscale, res_dscale + 4);
-	global_rscale = Min(global_rscale, NUMERIC_MAX_RESULT_SCALE);
-
 	/*
-	 * Call log_var() to compute and return the result
+	 * Call log_var() to compute and return the result; note it handles
+	 * rscale/dscale itself.
 	 */
 	log_var(&arg1, &arg2, &result);
 
-	result.dscale = res_dscale;
-
 	res = make_result(&result);
 
 	free_var(&result);
@@ -1378,7 +1418,7 @@ numeric_log(PG_FUNCTION_ARGS)
 /* ----------
  * numeric_power() -
  *
- *	Raise m to the power of x
+ *	Raise b to the power of x
  * ----------
  */
 Datum
@@ -1390,7 +1430,6 @@ numeric_power(PG_FUNCTION_ARGS)
 	NumericVar	arg1;
 	NumericVar	arg2;
 	NumericVar	result;
-	int			res_dscale;
 
 	/*
 	 * Handle NaN
@@ -1399,27 +1438,21 @@ numeric_power(PG_FUNCTION_ARGS)
 		PG_RETURN_NUMERIC(make_result(&const_nan));
 
 	/*
-	 * Initialize things and calculate scales
+	 * Initialize things
 	 */
 	init_var(&arg1);
 	init_var(&arg2);
 	init_var(&result);
+
 	set_var_from_num(num1, &arg1);
 	set_var_from_num(num2, &arg2);
 
-	res_dscale = Max(arg1.dscale + arg2.dscale, NUMERIC_MIN_DISPLAY_SCALE);
-	res_dscale = Min(res_dscale, NUMERIC_MAX_DISPLAY_SCALE);
-	global_rscale = Max(arg1.rscale + arg2.rscale, NUMERIC_MIN_RESULT_SCALE);
-	global_rscale = Max(global_rscale, res_dscale + 4);
-	global_rscale = Min(global_rscale, NUMERIC_MAX_RESULT_SCALE);
-
 	/*
-	 * Call log_var() to compute and return the result
+	 * Call power_var() to compute and return the result; note it handles
+	 * rscale/dscale itself.
 	 */
 	power_var(&arg1, &arg2, &result);
 
-	result.dscale = res_dscale;
-
 	res = make_result(&result);
 
 	free_var(&result);
@@ -1640,30 +1673,16 @@ Datum
 numeric_float8_no_overflow(PG_FUNCTION_ARGS)
 {
 	Numeric		num = PG_GETARG_NUMERIC(0);
-	char	   *tmp;
 	double		val;
-	char	   *endptr;
 
 	if (NUMERIC_IS_NAN(num))
 		PG_RETURN_FLOAT8(NAN);
 
-	tmp = DatumGetCString(DirectFunctionCall1(numeric_out,
-											  NumericGetDatum(num)));
-
-	/* unlike float8in, we ignore ERANGE from strtod */
-	val = strtod(tmp, &endptr);
-	if (*endptr != '\0')
-	{
-		/* shouldn't happen ... */
-		elog(ERROR, "Bad float8 input format '%s'", tmp);
-	}
-
-	pfree(tmp);
+	val = numeric_to_double_no_overflow(num);
 
 	PG_RETURN_FLOAT8(val);
 }
 
-
 Datum
 float4_numeric(PG_FUNCTION_ARGS)
 {
@@ -1939,6 +1958,9 @@ numeric_variance(PG_FUNCTION_ARGS)
 	init_var(&vsumX2);
 	set_var_from_num(sumX2, &vsumX2);
 
+	/* set rscale for mul_var calls */
+	global_rscale = vsumX.rscale * 2;
+
 	mul_var(&vsumX, &vsumX, &vsumX);	/* now vsumX contains sumX * sumX */
 	mul_var(&vN, &vsumX2, &vsumX2);		/* now vsumX2 contains N * sumX2 */
 	sub_var(&vsumX2, &vsumX, &vsumX2);	/* N * sumX2 - sumX * sumX */
@@ -2018,6 +2040,9 @@ numeric_stddev(PG_FUNCTION_ARGS)
 	init_var(&vsumX2);
 	set_var_from_num(sumX2, &vsumX2);
 
+	/* set rscale for mul_var calls */
+	global_rscale = vsumX.rscale * 2;
+
 	mul_var(&vsumX, &vsumX, &vsumX);	/* now vsumX contains sumX * sumX */
 	mul_var(&vN, &vsumX2, &vsumX2);		/* now vsumX2 contains N * sumX2 */
 	sub_var(&vsumX2, &vsumX, &vsumX2);	/* N * sumX2 - sumX * sumX */
@@ -2785,6 +2810,54 @@ apply_typmod(NumericVar *var, int32 typmod)
 	var->dscale = scale;
 }
 
+/* Convert numeric to float8; if out of range, return +/- HUGE_VAL */
+/* Caller should have eliminated the possibility of NAN */
+static double
+numeric_to_double_no_overflow(Numeric num)
+{
+	char	   *tmp;
+	double		val;
+	char	   *endptr;
+
+	tmp = DatumGetCString(DirectFunctionCall1(numeric_out,
+											  NumericGetDatum(num)));
+
+	/* unlike float8in, we ignore ERANGE from strtod */
+	val = strtod(tmp, &endptr);
+	if (*endptr != '\0')
+	{
+		/* shouldn't happen ... */
+		elog(ERROR, "Bad float8 input format '%s'", tmp);
+	}
+
+	pfree(tmp);
+
+	return val;
+}
+
+/* As above, but work from a NumericVar */
+static double
+numericvar_to_double_no_overflow(NumericVar *var)
+{
+	char	   *tmp;
+	double		val;
+	char	   *endptr;
+
+	tmp = get_str_from_var(var, var->dscale);
+
+	/* unlike float8in, we ignore ERANGE from strtod */
+	val = strtod(tmp, &endptr);
+	if (*endptr != '\0')
+	{
+		/* shouldn't happen ... */
+		elog(ERROR, "Bad float8 input format '%s'", tmp);
+	}
+
+	pfree(tmp);
+
+	return val;
+}
+
 
 /* ----------
  * cmp_var() -
@@ -3073,7 +3146,7 @@ sub_var(NumericVar *var1, NumericVar *var2, NumericVar *result)
  * mul_var() -
  *
  *	Multiplication on variable level. Product of var1 * var2 is stored
- *	in result.
+ *	in result.  Accuracy of result is determined by global_rscale.
  * ----------
  */
 static void
@@ -3158,7 +3231,8 @@ mul_var(NumericVar *var1, NumericVar *var2, NumericVar *result)
 /* ----------
  * div_var() -
  *
- *	Division on variable level.
+ *	Division on variable level.  Accuracy of result is determined by
+ *	global_rscale.
  * ----------
  */
 static void
@@ -3370,33 +3444,69 @@ div_var(NumericVar *var1, NumericVar *var2, NumericVar *result)
 static int
 select_div_scale(NumericVar *var1, NumericVar *var2)
 {
+	int			weight1,
+				weight2,
+				qweight,
+				i;
+	NumericDigit firstdigit1,
+				firstdigit2;
 	int			res_dscale;
 	int			res_rscale;
 
-	/* ----------
-	 * The result scale of a division isn't specified in any
-	 * SQL standard. For Postgres it is the following (where
-	 * SR, DR are the result- and display-scales of the returned
-	 * value, S1, D1, S2 and D2 are the scales of the two arguments,
-	 * The minimum and maximum scales are compile time options from
-	 * numeric.h):
-	 *
-	 *	DR = Min(Max(D1 + D2, MIN_DISPLAY_SCALE), MAX_DISPLAY_SCALE)
-	 *	SR = Min(Max(Max(S1 + S2, DR + 4), MIN_RESULT_SCALE), MAX_RESULT_SCALE)
+	/*
+	 * The result scale of a division isn't specified in any SQL standard.
+	 * For PostgreSQL we select a display scale that will give at least
+	 * NUMERIC_MIN_SIG_DIGITS significant digits, so that numeric gives a
+	 * result no less accurate than float8; but use a scale not less than
+	 * either input's display scale.
 	 *
-	 * By default, any result is computed with a minimum of 34 digits
-	 * after the decimal point or at least with 4 digits more than
-	 * displayed.
-	 * ----------
+	 * The result scale is NUMERIC_EXTRA_DIGITS more than the display scale,
+	 * to provide some guard digits in the calculation.
+	 */
+
+	/* Get the actual (normalized) weight and first digit of each input */
+
+	weight1 = 0;				/* values to use if var1 is zero */
+	firstdigit1 = 0;
+	for (i = 0; i < var1->ndigits; i++)
+	{
+		firstdigit1 = var1->digits[i];
+		if (firstdigit1 != 0)
+		{
+			weight1 = var1->weight - i;
+			break;
+		}
+	}
+
+	weight2 = 0;				/* values to use if var2 is zero */
+	firstdigit2 = 0;
+	for (i = 0; i < var2->ndigits; i++)
+	{
+		firstdigit2 = var2->digits[i];
+		if (firstdigit2 != 0)
+		{
+			weight2 = var2->weight - i;
+			break;
+		}
+	}
+
+	/*
+	 * Estimate weight of quotient.  If the two first digits are equal,
+	 * we can't be sure, but assume that var1 is less than var2.
 	 */
-	res_dscale = var1->dscale + var2->dscale;
+	qweight = weight1 - weight2;
+	if (firstdigit1 <= firstdigit2)
+		qweight--;
+
+	/* Select display scale */
+	res_dscale = NUMERIC_MIN_SIG_DIGITS - qweight;
+	res_dscale = Max(res_dscale, var1->dscale);
+	res_dscale = Max(res_dscale, var2->dscale);
 	res_dscale = Max(res_dscale, NUMERIC_MIN_DISPLAY_SCALE);
 	res_dscale = Min(res_dscale, NUMERIC_MAX_DISPLAY_SCALE);
 
-	res_rscale = var1->rscale + var2->rscale;
-	res_rscale = Max(res_rscale, res_dscale + 4);
-	res_rscale = Max(res_rscale, NUMERIC_MIN_RESULT_SCALE);
-	res_rscale = Min(res_rscale, NUMERIC_MAX_RESULT_SCALE);
+	/* Select result scale */
+	res_rscale = res_dscale + NUMERIC_EXTRA_DIGITS;
 	global_rscale = res_rscale;
 
 	return res_dscale;
@@ -3503,7 +3613,7 @@ floor_var(NumericVar *var, NumericVar *result)
 /* ----------
  * sqrt_var() -
  *
- *	Compute the square root of x using Newtons algorithm
+ *	Compute the square root of x using Newton's algorithm
  * ----------
  */
 static void
@@ -3526,6 +3636,7 @@ sqrt_var(NumericVar *arg, NumericVar *result)
 		set_var_from_var(&const_zero, result);
 		result->rscale = res_rscale;
 		result->sign = NUMERIC_POS;
+		global_rscale = save_global_rscale;
 		return;
 	}
 
@@ -3536,8 +3647,8 @@ sqrt_var(NumericVar *arg, NumericVar *result)
 	init_var(&tmp_val);
 	init_var(&last_val);
 
+	/* Copy arg in case it is the same var as result */
 	set_var_from_var(arg, &tmp_arg);
-	set_var_from_var(result, &last_val);
 
 	/*
 	 * Initialize the result to the first guess
@@ -3553,6 +3664,8 @@ sqrt_var(NumericVar *arg, NumericVar *result)
 	result->rscale = res_rscale;
 	result->sign = NUMERIC_POS;
 
+	set_var_from_var(result, &last_val);
+
 	for (;;)
 	{
 		div_var(&tmp_arg, result, &tmp_val);
@@ -3590,6 +3703,7 @@ exp_var(NumericVar *arg, NumericVar *result)
 	NumericVar	ni;
 	int			d;
 	int			i;
+	int			xintval;
 	int			ndiv2 = 0;
 	bool		xneg = FALSE;
 	int			save_global_rscale;
@@ -3608,32 +3722,42 @@ exp_var(NumericVar *arg, NumericVar *result)
 		x.sign = NUMERIC_POS;
 	}
 
-	save_global_rscale = global_rscale;
-	global_rscale = 0;
+	/* Select an appropriate scale for internal calculation */
+	xintval = 0;
 	for (i = x.weight, d = 0; i >= 0; i--, d++)
 	{
-		global_rscale *= 10;
+		xintval *= 10;
 		if (d < x.ndigits)
-			global_rscale += x.digits[d];
-		if (global_rscale >= 1000)
+			xintval += x.digits[d];
+		if (xintval >= NUMERIC_MAX_RESULT_SCALE)
 			elog(ERROR, "argument for EXP() too big");
 	}
 
-	global_rscale = global_rscale / 2 + save_global_rscale + 8;
+	save_global_rscale = global_rscale;
+	global_rscale += xintval / 2 + 8;
 
-	while (cmp_var(&x, &const_one) > 0)
+	/* Reduce input into range 0 <= x <= 0.1 */
+	while (cmp_var(&x, &const_zero_point_one) > 0)
 	{
 		ndiv2++;
 		global_rscale++;
 		div_var(&x, &const_two, &x);
 	}
 
+	/*
+	 * Use the Taylor series
+	 *
+	 *		exp(x) = 1 + x + x^2/2! + x^3/3! + ...
+	 *
+	 * Given the limited range of x, this should converge reasonably quickly.
+	 * We run the series until the terms fall below the global_rscale limit.
+	 */
 	add_var(&const_one, &x, result);
 	set_var_from_var(&x, &xpow);
 	set_var_from_var(&const_one, &ifac);
 	set_var_from_var(&const_one, &ni);
 
-	for (i = 2;; i++)
+	for (;;)
 	{
 		add_var(&ni, &const_one, &ni);
 		mul_var(&xpow, &x, &xpow);
@@ -3646,10 +3770,13 @@ exp_var(NumericVar *arg, NumericVar *result)
 		add_var(result, &elem, result);
 	}
 
+	/* Compensate for argument range reduction */
 	while (ndiv2-- > 0)
 		mul_var(result, result, result);
 
+	/* Compensate for input sign, and round to caller's global_rscale */
 	global_rscale = save_global_rscale;
+
 	if (xneg)
 		div_var(&const_one, result, result);
 	else
@@ -3679,7 +3806,6 @@ ln_var(NumericVar *arg, NumericVar *result)
 	NumericVar	ni;
 	NumericVar	elem;
 	NumericVar	fact;
-	int			i;
 	int			save_global_rscale;
 
 	if (cmp_var(arg, &const_zero) <= 0)
@@ -3697,18 +3823,31 @@ ln_var(NumericVar *arg, NumericVar *result)
 	set_var_from_var(&const_two, &fact);
 	set_var_from_var(arg, &x);
 
-	while (cmp_var(&x, &const_two) >= 0)
+	/* Reduce input into range 0.9 < x < 1.1 */
+	while (cmp_var(&x, &const_zero_point_nine) <= 0)
 	{
+		global_rscale++;
 		sqrt_var(&x, &x);
 		mul_var(&fact, &const_two, &fact);
 	}
-	set_var_from_str("0.5", &elem);
-	while (cmp_var(&x, &elem) <= 0)
+	while (cmp_var(&x, &const_one_point_one) >= 0)
 	{
+		global_rscale++;
 		sqrt_var(&x, &x);
 		mul_var(&fact, &const_two, &fact);
 	}
 
+	/*
+	 * We use the Taylor series for 0.5 * ln((1+z)/(1-z)),
+	 *
+	 *		z + z^3/3 + z^5/5 + ...
+	 *
+	 * where z = (x-1)/(x+1) is in the range (approximately) -0.053 .. 0.048
+	 * due to the above range-reduction of x.
+	 *
+	 * The convergence of this is not as fast as one would like, but is
+	 * tolerable given that z is small.
+	 */
 	sub_var(&x, &const_one, result);
 	add_var(&x, &const_one, &elem);
 	div_var(result, &elem, result);
@@ -3717,19 +3856,21 @@ ln_var(NumericVar *arg, NumericVar *result)
 
 	set_var_from_var(&const_one, &ni);
 
-	for (i = 2;; i++)
+	for (;;)
 	{
 		add_var(&ni, &const_two, &ni);
 		mul_var(&xx, &x, &xx);
 		div_var(&xx, &ni, &elem);
 
-		if (cmp_var(&elem, &const_zero) == 0)
+		if (elem.ndigits == 0)
 			break;
 
 		add_var(result, &elem, result);
 	}
 
+	/* Compensate for argument range reduction, round to caller's rscale */
 	global_rscale = save_global_rscale;
+
 	mul_var(result, &fact, result);
 
 	free_var(&x);
@@ -3743,7 +3884,9 @@ ln_var(NumericVar *arg, NumericVar *result)
 /* ----------
  * log_var() -
  *
- *	Compute the logarithm of x in a given base
+ *	Compute the logarithm of num in a given base.
+ *
+ *	Note: this routine chooses rscale and dscale of the result.
  * ----------
  */
 static void
@@ -3751,19 +3894,43 @@ log_var(NumericVar *base, NumericVar *num, NumericVar *result)
 {
 	NumericVar	ln_base;
 	NumericVar	ln_num;
-
-	global_rscale += 8;
+	int			save_global_rscale = global_rscale;
+	int			res_dscale;
 
 	init_var(&ln_base);
 	init_var(&ln_num);
 
+	/* Set scale for ln() calculations */
+	if (num->weight > 0)
+		res_dscale = NUMERIC_MIN_SIG_DIGITS - (int) log10(num->weight);
+	else if (num->weight < 0)
+		res_dscale = NUMERIC_MIN_SIG_DIGITS - (int) log10(- num->weight);
+	else
+		res_dscale = NUMERIC_MIN_SIG_DIGITS;
+
+	res_dscale = Max(res_dscale, base->dscale);
+	res_dscale = Max(res_dscale, num->dscale);
+	res_dscale = Max(res_dscale, NUMERIC_MIN_DISPLAY_SCALE);
+	res_dscale = Min(res_dscale, NUMERIC_MAX_DISPLAY_SCALE);
+
+	global_rscale = res_dscale + 8;
+
+	/* Form natural logarithms */
 	ln_var(base, &ln_base);
 	ln_var(num, &ln_num);
 
-	global_rscale -= 8;
+	ln_base.dscale = res_dscale;
+	ln_num.dscale = res_dscale;
+
+	/* Select scale for division result */
+	res_dscale = select_div_scale(&ln_num, &ln_base);
 
 	div_var(&ln_num, &ln_base, result);
 
+	result->dscale = res_dscale;
+
+	global_rscale = save_global_rscale;
+
 	free_var(&ln_num);
 	free_var(&ln_base);
 }
@@ -3773,6 +3940,8 @@ log_var(NumericVar *base, NumericVar *num, NumericVar *result)
  * power_var() -
  *
  *	Raise base to the power of exp
+ *
+ *	Note: this routine chooses rscale and dscale of the result.
  * ----------
  */
 static void
@@ -3780,24 +3949,64 @@ power_var(NumericVar *base, NumericVar *exp, NumericVar *result)
 {
 	NumericVar	ln_base;
 	NumericVar	ln_num;
-	int			save_global_rscale;
-
-	save_global_rscale = global_rscale;
-	global_rscale += global_rscale / 3 + 8;
+	int			save_global_rscale = global_rscale;
+	int			res_dscale;
+	double		val;
 
 	init_var(&ln_base);
 	init_var(&ln_num);
 
+	/* Set scale for ln() calculation --- need extra accuracy here */
+	if (base->weight > 0)
+		res_dscale = NUMERIC_MIN_SIG_DIGITS*2 - (int) log10(base->weight);
+	else if (base->weight < 0)
+		res_dscale = NUMERIC_MIN_SIG_DIGITS*2 - (int) log10(- base->weight);
+	else
+		res_dscale = NUMERIC_MIN_SIG_DIGITS*2;
+
+	res_dscale = Max(res_dscale, base->dscale * 2);
+	res_dscale = Max(res_dscale, exp->dscale * 2);
+	res_dscale = Max(res_dscale, NUMERIC_MIN_DISPLAY_SCALE);
+	res_dscale = Min(res_dscale, NUMERIC_MAX_DISPLAY_SCALE);
+
+	global_rscale = res_dscale + 8;
+
 	ln_var(base, &ln_base);
+
+	ln_base.dscale = res_dscale;
+
 	mul_var(&ln_base, exp, &ln_num);
 
-	global_rscale = save_global_rscale;
+	ln_num.dscale = res_dscale;
+
+	/* Set scale for exp() */
+
+	/* convert input to float8, ignoring overflow */
+	val = numericvar_to_double_no_overflow(&ln_num);
+
+	/* log10(result) = num * log10(e), so this is approximately the weight: */
+	val *= 0.434294481903252;
+
+	/* limit to something that won't cause integer overflow */
+	val = Max(val, -NUMERIC_MAX_RESULT_SCALE);
+	val = Min(val, NUMERIC_MAX_RESULT_SCALE);
+
+	res_dscale = NUMERIC_MIN_SIG_DIGITS - (int) val;
+	res_dscale = Max(res_dscale, base->dscale);
+	res_dscale = Max(res_dscale, exp->dscale);
+	res_dscale = Max(res_dscale, NUMERIC_MIN_DISPLAY_SCALE);
+	res_dscale = Min(res_dscale, NUMERIC_MAX_DISPLAY_SCALE);
+
+	global_rscale = res_dscale + 8;
 
 	exp_var(&ln_num, result);
 
+	result->dscale = res_dscale;
+
+	global_rscale = save_global_rscale;
+
 	free_var(&ln_num);
 	free_var(&ln_base);
-
 }
 
 
diff --git a/src/include/utils/numeric.h b/src/include/utils/numeric.h
index f6050dfc65995e2fa977de119e13b424fc6d8785..a8148d59cea1697fe87080acb7c99dd08745fcf8 100644
--- a/src/include/utils/numeric.h
+++ b/src/include/utils/numeric.h
@@ -5,7 +5,7 @@
  *
  *	1998 Jan Wieck
  *
- * $Header: /cvsroot/pgsql/src/include/utils/numeric.h,v 1.15 2001/11/05 17:46:36 momjian Exp $
+ * $Id: numeric.h,v 1.16 2002/10/02 19:21:26 tgl Exp $
  *
  * ----------
  */
@@ -13,29 +13,36 @@
 #ifndef _PG_NUMERIC_H_
 #define _PG_NUMERIC_H_
 
-/* ----------
- * The hardcoded limits and defaults of the numeric data type
- * ----------
+/*
+ * Hardcoded precision limit - arbitrary, but must be small enough that
+ * dscale values will fit in 14 bits.
  */
 #define NUMERIC_MAX_PRECISION		1000
-#define NUMERIC_DEFAULT_PRECISION	30
-#define NUMERIC_DEFAULT_SCALE		6
 
-
-/* ----------
+/*
  * Internal limits on the scales chosen for calculation results
- * ----------
  */
 #define NUMERIC_MAX_DISPLAY_SCALE	NUMERIC_MAX_PRECISION
-#define NUMERIC_MIN_DISPLAY_SCALE	(NUMERIC_DEFAULT_SCALE + 4)
+#define NUMERIC_MIN_DISPLAY_SCALE	0
 
 #define NUMERIC_MAX_RESULT_SCALE	(NUMERIC_MAX_PRECISION * 2)
-#define NUMERIC_MIN_RESULT_SCALE	(NUMERIC_DEFAULT_PRECISION + 4)
 
+/*
+ * For inherently inexact calculations such as division and square root,
+ * we try to get at least this many significant digits; the idea is to
+ * deliver a result no worse than float8 would.
+ */
+#define NUMERIC_MIN_SIG_DIGITS		16
 
-/* ----------
+/*
+ * Standard number of extra digits carried internally while doing
+ * inexact calculations.
+ */
+#define NUMERIC_EXTRA_DIGITS		4
+
+
+/*
  * Sign values and macros to deal with packing/unpacking n_sign_dscale
- * ----------
  */
 #define NUMERIC_SIGN_MASK	0xC000
 #define NUMERIC_POS			0x0000
@@ -48,7 +55,7 @@
 								NUMERIC_SIGN(n) != NUMERIC_NEG)
 
 
-/* ----------
+/*
  * The Numeric data type stored in the database
  *
  * NOTE: by convention, values in the packed form have been stripped of
@@ -56,7 +63,6 @@
  * in the last byte, if the number of digits is odd).  In particular,
  * if the value is zero, there will be no digits at all!  The weight is
  * arbitrary in that case, but we normally set it to zero.
- * ----------
  */
 typedef struct NumericData
 {
diff --git a/src/test/regress/expected/aggregates.out b/src/test/regress/expected/aggregates.out
index aaee72c01abdacc6520ba49a573041799875d226..7519f6e5392a9ee8d45a7e4ae94b54b70d0d2e89 100644
--- a/src/test/regress/expected/aggregates.out
+++ b/src/test/regress/expected/aggregates.out
@@ -2,15 +2,15 @@
 -- AGGREGATES
 --
 SELECT avg(four) AS avg_1 FROM onek;
-    avg_1     
---------------
- 1.5000000000
+        avg_1        
+---------------------
+ 1.50000000000000000
 (1 row)
 
 SELECT avg(a) AS avg_32 FROM aggtest WHERE a < 100;
-    avg_32     
----------------
- 32.6666666667
+       avg_32       
+--------------------
+ 32.666666666666667
 (1 row)
 
 -- In 7.1, avg(float4) is computed using float8 arithmetic.
@@ -118,9 +118,9 @@ select ten, count(four), sum(DISTINCT four) from onek group by ten;
 (10 rows)
 
 SELECT newavg(four) AS avg_1 FROM onek;
-    avg_1     
---------------
- 1.5000000000
+        avg_1        
+---------------------
+ 1.50000000000000000
 (1 row)
 
 SELECT newsum(four) AS sum_1500 FROM onek;