1 files changed, 377 insertions, 227 deletions

diff --git a/src/backend/utils/adt/numeric.c b/src/backend/utils/adt/numeric.c
index 1667d8093f1..fe9f3f7a506 100644
--- a/src/backend/utils/adt/numeric.c
+++ b/src/backend/utils/adt/numeric.c
@@ -224,6 +224,7 @@ struct NumericData
  *
  * The value represented by a NumericVar is determined by the sign, weight,
  * ndigits, and digits[] array.
+ *
  * Note: the first digit of a NumericVar's value is assumed to be multiplied
  * by NBASE ** weight.  Another way to say it is that there are weight+1
  * digits before the decimal point.  It is possible to have weight < 0.
@@ -255,6 +256,11 @@ struct NumericData
  * Note that rscale is not stored in variables --- it's figured on-the-fly
  * from the dscales of the inputs.
  *
+ * While we consistently use "weight" to refer to the base-NBASE weight of
+ * a numeric value, it is convenient in some scale-related calculations to
+ * make use of the base-10 weight (ie, the approximate log10 of the value).
+ * To avoid confusion, such a decimal-units weight is called a "dweight".
+ *
  * NB: All the variable-level functions are written in a style that makes it
  * possible to give one and the same variable as argument and destination.
  * This is feasible because the digit buffer is separate from the variable.
@@ -361,20 +367,6 @@ static NumericVar const_zero_point_nine =
 {1, -1, NUMERIC_POS, 1, NULL, const_zero_point_nine_data};
 
 #if DEC_DIGITS == 4
-static NumericDigit const_zero_point_01_data[1] = {100};
-static NumericVar const_zero_point_01 =
-{1, -1, NUMERIC_POS, 2, NULL, const_zero_point_01_data};
-#elif DEC_DIGITS == 2
-static NumericDigit const_zero_point_01_data[1] = {1};
-static NumericVar const_zero_point_01 =
-{1, -1, NUMERIC_POS, 2, NULL, const_zero_point_01_data};
-#elif DEC_DIGITS == 1
-static NumericDigit const_zero_point_01_data[1] = {1};
-static NumericVar const_zero_point_01 =
-{1, -2, NUMERIC_POS, 2, NULL, const_zero_point_01_data};
-#endif
-
-#if DEC_DIGITS == 4
 static NumericDigit const_one_point_one_data[2] = {1, 1000};
 #elif DEC_DIGITS == 2
 static NumericDigit const_one_point_one_data[2] = {1, 10};
@@ -477,7 +469,7 @@ static void floor_var(NumericVar *var, NumericVar *result);
 
 static void sqrt_var(NumericVar *arg, NumericVar *result, int rscale);
 static void exp_var(NumericVar *arg, NumericVar *result, int rscale);
-static void exp_var_internal(NumericVar *arg, NumericVar *result, int rscale);
+static int	estimate_ln_dweight(NumericVar *var);
 static void ln_var(NumericVar *arg, NumericVar *result, int rscale);
 static void log_var(NumericVar *base, NumericVar *num, NumericVar *result);
 static void power_var(NumericVar *base, NumericVar *exp, NumericVar *result);
@@ -2697,7 +2689,7 @@ numeric_ln(PG_FUNCTION_ARGS)
 	Numeric		res;
 	NumericVar	arg;
 	NumericVar	result;
-	int			dec_digits;
+	int			ln_dweight;
 	int			rscale;
 
 	/*
@@ -2709,16 +2701,10 @@ numeric_ln(PG_FUNCTION_ARGS)
 	init_var_from_num(num, &arg);
 	init_var(&result);
 
-	/* Approx decimal digits before decimal point */
-	dec_digits = (arg.weight + 1) * DEC_DIGITS;
-
-	if (dec_digits > 1)
-		rscale = NUMERIC_MIN_SIG_DIGITS - (int) log10(dec_digits - 1);
-	else if (dec_digits < 1)
-		rscale = NUMERIC_MIN_SIG_DIGITS - (int) log10(1 - dec_digits);
-	else
-		rscale = NUMERIC_MIN_SIG_DIGITS;
+	/* Estimated dweight of logarithm */
+	ln_dweight = estimate_ln_dweight(&arg);
 
+	rscale = NUMERIC_MIN_SIG_DIGITS - ln_dweight;
 	rscale = Max(rscale, arg.dscale);
 	rscale = Max(rscale, NUMERIC_MIN_DISPLAY_SCALE);
 	rscale = Min(rscale, NUMERIC_MAX_DISPLAY_SCALE);
@@ -5720,17 +5706,35 @@ mul_var(NumericVar *var1, NumericVar *var2, NumericVar *result,
 	int			carry;
 	int			maxdig;
 	int			newdig;
+	int			var1ndigits;
+	int			var2ndigits;
+	NumericDigit *var1digits;
+	NumericDigit *var2digits;
 	NumericDigit *res_digits;
 	int			i,
-				ri,
 				i1,
 				i2;
 
+	/*
+	 * Arrange for var1 to be the shorter of the two numbers.  This improves
+	 * performance because the inner multiplication loop is much simpler than
+	 * the outer loop, so it's better to have a smaller number of iterations
+	 * of the outer loop.  This also reduces the number of times that the
+	 * accumulator array needs to be normalized.
+	 */
+	if (var1->ndigits > var2->ndigits)
+	{
+		NumericVar *tmp = var1;
+
+		var1 = var2;
+		var2 = tmp;
+	}
+
 	/* copy these values into local vars for speed in inner loop */
-	int			var1ndigits = var1->ndigits;
-	int			var2ndigits = var2->ndigits;
-	NumericDigit *var1digits = var1->digits;
-	NumericDigit *var2digits = var2->digits;
+	var1ndigits = var1->ndigits;
+	var2ndigits = var2->ndigits;
+	var1digits = var1->digits;
+	var2digits = var2->digits;
 
 	if (var1ndigits == 0 || var2ndigits == 0)
 	{
@@ -5748,39 +5752,28 @@ mul_var(NumericVar *var1, NumericVar *var2, NumericVar *result,
 	res_weight = var1->weight + var2->weight + 2;
 
 	/*
-	 * Determine number of result digits to compute.  If the exact result
+	 * Determine the number of result digits to compute.  If the exact result
 	 * would have more than rscale fractional digits, truncate the computation
-	 * with MUL_GUARD_DIGITS guard digits.  We do that by pretending that one
-	 * or both inputs have fewer digits than they really do.
+	 * with MUL_GUARD_DIGITS guard digits, i.e., ignore input digits that
+	 * would only contribute to the right of that.  (This will give the exact
+	 * rounded-to-rscale answer unless carries out of the ignored positions
+	 * would have propagated through more than MUL_GUARD_DIGITS digits.)
+	 *
+	 * Note: an exact computation could not produce more than var1ndigits +
+	 * var2ndigits digits, but we allocate one extra output digit in case
+	 * rscale-driven rounding produces a carry out of the highest exact digit.
 	 */
 	res_ndigits = var1ndigits + var2ndigits + 1;
-	maxdigits = res_weight + 1 + (rscale * DEC_DIGITS) + MUL_GUARD_DIGITS;
-	if (res_ndigits > maxdigits)
+	maxdigits = res_weight + 1 + (rscale + DEC_DIGITS - 1) / DEC_DIGITS +
+		MUL_GUARD_DIGITS;
+	res_ndigits = Min(res_ndigits, maxdigits);
+
+	if (res_ndigits < 3)
 	{
-		if (maxdigits < 3)
-		{
-			/* no useful precision at all in the result... */
-			zero_var(result);
-			result->dscale = rscale;
-			return;
-		}
-		/* force maxdigits odd so that input ndigits can be equal */
-		if ((maxdigits & 1) == 0)
-			maxdigits++;
-		if (var1ndigits > var2ndigits)
-		{
-			var1ndigits -= res_ndigits - maxdigits;
-			if (var1ndigits < var2ndigits)
-				var1ndigits = var2ndigits = (var1ndigits + var2ndigits) / 2;
-		}
-		else
-		{
-			var2ndigits -= res_ndigits - maxdigits;
-			if (var2ndigits < var1ndigits)
-				var1ndigits = var2ndigits = (var1ndigits + var2ndigits) / 2;
-		}
-		res_ndigits = maxdigits;
-		Assert(res_ndigits == var1ndigits + var2ndigits + 1);
+		/* All input digits will be ignored; so result is zero */
+		zero_var(result);
+		result->dscale = rscale;
+		return;
 	}
 
 	/*
@@ -5802,8 +5795,16 @@ mul_var(NumericVar *var1, NumericVar *var2, NumericVar *result,
 	dig = (int *) palloc0(res_ndigits * sizeof(int));
 	maxdig = 0;
 
-	ri = res_ndigits - 1;
-	for (i1 = var1ndigits - 1; i1 >= 0; ri--, i1--)
+	/*
+	 * The least significant digits of var1 should be ignored if they don't
+	 * contribute directly to the first res_ndigits digits of the result that
+	 * we are computing.
+	 *
+	 * Digit i1 of var1 and digit i2 of var2 are multiplied and added to digit
+	 * i1+i2+2 of the accumulator array, so we need only consider digits of
+	 * var1 for which i1 <= res_ndigits - 3.
+	 */
+	for (i1 = Min(var1ndigits - 1, res_ndigits - 3); i1 >= 0; i1--)
 	{
 		int			var1digit = var1digits[i1];
 
@@ -5833,9 +5834,14 @@ mul_var(NumericVar *var1, NumericVar *var2, NumericVar *result,
 			maxdig = 1 + var1digit;
 		}
 
-		/* Add appropriate multiple of var2 into the accumulator */
-		i = ri;
-		for (i2 = var2ndigits - 1; i2 >= 0; i2--)
+		/*
+		 * Add the appropriate multiple of var2 into the accumulator.
+		 *
+		 * As above, digits of var2 can be ignored if they don't contribute,
+		 * so we only include digits for which i1+i2+2 <= res_ndigits - 1.
+		 */
+		for (i2 = Min(var2ndigits - 1, res_ndigits - i1 - 3), i = i1 + i2 + 2;
+			 i2 >= 0; i2--)
 			dig[i--] += var1digit * var2digits[i2];
 	}
 
@@ -6635,120 +6641,81 @@ sqrt_var(NumericVar *arg, NumericVar *result, int rscale)
 /*
  * exp_var() -
  *
- *	Raise e to the power of x
+ *	Raise e to the power of x, computed to rscale fractional digits
  */
 static void
 exp_var(NumericVar *arg, NumericVar *result, int rscale)
 {
 	NumericVar	x;
-	int			xintval;
-	bool		xneg = FALSE;
+	NumericVar	elem;
+	NumericVar	ni;
+	double		val;
+	int			dweight;
+	int			ndiv2;
+	int			sig_digits;
 	int			local_rscale;
 
-	/*----------
-	 * We separate the integral and fraction parts of x, then compute
-	 *		e^x = e^xint * e^xfrac
-	 * where e = exp(1) and e^xfrac = exp(xfrac) are computed by
-	 * exp_var_internal; the limited range of inputs allows that routine
-	 * to do a good job with a simple Taylor series.  Raising e^xint is
-	 * done by repeated multiplications in power_var_int.
-	 *----------
-	 */
 	init_var(&x);
+	init_var(&elem);
+	init_var(&ni);
 
 	set_var_from_var(arg, &x);
 
-	if (x.sign == NUMERIC_NEG)
-	{
-		xneg = TRUE;
-		x.sign = NUMERIC_POS;
-	}
-
-	/* Extract the integer part, remove it from x */
-	xintval = 0;
-	while (x.weight >= 0)
-	{
-		xintval *= NBASE;
-		if (x.ndigits > 0)
-		{
-			xintval += x.digits[0];
-			x.digits++;
-			x.ndigits--;
-		}
-		x.weight--;
-		/* Guard against overflow */
-		if (xintval >= NUMERIC_MAX_RESULT_SCALE * 3)
-			ereport(ERROR,
-					(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
-					 errmsg("argument for function \"exp\" too big")));
-	}
+	/*
+	 * Estimate the dweight of the result using floating point arithmetic, so
+	 * that we can choose an appropriate local rscale for the calculation.
+	 */
+	val = numericvar_to_double_no_overflow(&x);
 
-	/* Select an appropriate scale for internal calculation */
-	local_rscale = rscale + MUL_GUARD_DIGITS * 2;
+	/* Guard against overflow */
+	if (Abs(val) >= NUMERIC_MAX_RESULT_SCALE * 3)
+		ereport(ERROR,
+				(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
+				 errmsg("value overflows numeric format")));
 
-	/* Compute e^xfrac */
-	exp_var_internal(&x, result, local_rscale);
+	/* decimal weight = log10(e^x) = x * log10(e) */
+	dweight = (int) (val * 0.434294481903252);
 
-	/* If there's an integer part, multiply by e^xint */
-	if (xintval > 0)
+	/*
+	 * Reduce x to the range -0.01 <= x <= 0.01 (approximately) by dividing by
+	 * 2^n, to improve the convergence rate of the Taylor series.
+	 */
+	if (Abs(val) > 0.01)
 	{
-		NumericVar	e;
+		NumericVar	tmp;
 
-		init_var(&e);
-		exp_var_internal(&const_one, &e, local_rscale);
-		power_var_int(&e, xintval, &e, local_rscale);
-		mul_var(&e, result, result, local_rscale);
-		free_var(&e);
-	}
+		init_var(&tmp);
+		set_var_from_var(&const_two, &tmp);
 
-	/* Compensate for input sign, and round to requested rscale */
-	if (xneg)
-		div_var_fast(&const_one, result, result, rscale, true);
-	else
-		round_var(result, rscale);
-
-	free_var(&x);
-}
-
-
-/*
- * exp_var_internal() -
- *
- *	Raise e to the power of x, where 0 <= x <= 1
- *
- * NB: the result should be good to at least rscale digits, but it has
- * *not* been rounded off; the caller must do that if wanted.
- */
-static void
-exp_var_internal(NumericVar *arg, NumericVar *result, int rscale)
-{
-	NumericVar	x;
-	NumericVar	xpow;
-	NumericVar	ifac;
-	NumericVar	elem;
-	NumericVar	ni;
-	int			ndiv2 = 0;
-	int			local_rscale;
+		ndiv2 = 1;
+		val /= 2;
 
-	init_var(&x);
-	init_var(&xpow);
-	init_var(&ifac);
-	init_var(&elem);
-	init_var(&ni);
+		while (Abs(val) > 0.01)
+		{
+			ndiv2++;
+			val /= 2;
+			add_var(&tmp, &tmp, &tmp);
+		}
 
-	set_var_from_var(arg, &x);
+		local_rscale = x.dscale + ndiv2;
+		div_var_fast(&x, &tmp, &x, local_rscale, true);
 
-	Assert(x.sign == NUMERIC_POS);
+		free_var(&tmp);
+	}
+	else
+		ndiv2 = 0;
 
-	local_rscale = rscale + 8;
+	/*
+	 * Set the scale for the Taylor series expansion.  The final result has
+	 * (dweight + rscale + 1) significant digits.  In addition, we have to
+	 * raise the Taylor series result to the power 2^ndiv2, which introduces
+	 * an error of up to around log10(2^ndiv2) digits, so work with this many
+	 * extra digits of precision (plus a few more for good measure).
+	 */
+	sig_digits = 1 + dweight + rscale + (int) (ndiv2 * 0.301029995663981);
+	sig_digits = Max(sig_digits, 0) + 8;
 
-	/* Reduce input into range 0 <= x <= 0.01 */
-	while (cmp_var(&x, &const_zero_point_01) > 0)
-	{
-		ndiv2++;
-		local_rscale++;
-		mul_var(&x, &const_zero_point_five, &x, x.dscale + 1);
-	}
+	local_rscale = sig_digits - 1;
 
 	/*
 	 * Use the Taylor series
@@ -6759,36 +6726,122 @@ exp_var_internal(NumericVar *arg, NumericVar *result, int rscale)
 	 * We run the series until the terms fall below the local_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 (;;)
-	{
-		add_var(&ni, &const_one, &ni);
-		mul_var(&xpow, &x, &xpow, local_rscale);
-		mul_var(&ifac, &ni, &ifac, 0);
-		div_var_fast(&xpow, &ifac, &elem, local_rscale, true);
-
-		if (elem.ndigits == 0)
-			break;
+	mul_var(&x, &x, &elem, local_rscale);
+	set_var_from_var(&const_two, &ni);
+	div_var_fast(&elem, &ni, &elem, local_rscale, true);
 
+	while (elem.ndigits != 0)
+	{
 		add_var(result, &elem, result);
+
+		mul_var(&elem, &x, &elem, local_rscale);
+		add_var(&ni, &const_one, &ni);
+		div_var_fast(&elem, &ni, &elem, local_rscale, true);
 	}
 
-	/* Compensate for argument range reduction */
+	/*
+	 * Compensate for the argument range reduction.  Since the weight of the
+	 * result doubles with each multiplication, we can reduce the local rscale
+	 * as we proceed.
+	 */
 	while (ndiv2-- > 0)
+	{
+		local_rscale = sig_digits - result->weight * 2 * DEC_DIGITS;
+		local_rscale = Max(local_rscale, NUMERIC_MIN_DISPLAY_SCALE);
 		mul_var(result, result, result, local_rscale);
+	}
+
+	/* Round to requested rscale */
+	round_var(result, rscale);
 
 	free_var(&x);
-	free_var(&xpow);
-	free_var(&ifac);
 	free_var(&elem);
 	free_var(&ni);
 }
 
 
 /*
+ * Estimate the dweight of the most significant decimal digit of the natural
+ * logarithm of a number.
+ *
+ * Essentially, we're approximating log10(abs(ln(var))).  This is used to
+ * determine the appropriate rscale when computing natural logarithms.
+ */
+static int
+estimate_ln_dweight(NumericVar *var)
+{
+	int			ln_dweight;
+
+	if (cmp_var(var, &const_zero_point_nine) >= 0 &&
+		cmp_var(var, &const_one_point_one) <= 0)
+	{
+		/*
+		 * 0.9 <= var <= 1.1
+		 *
+		 * ln(var) has a negative weight (possibly very large).  To get a
+		 * reasonably accurate result, estimate it using ln(1+x) ~= x.
+		 */
+		NumericVar	x;
+
+		init_var(&x);
+		sub_var(var, &const_one, &x);
+
+		if (x.ndigits > 0)
+		{
+			/* Use weight of most significant decimal digit of x */
+			ln_dweight = x.weight * DEC_DIGITS + (int) log10(x.digits[0]);
+		}
+		else
+		{
+			/* x = 0.  Since ln(1) = 0 exactly, we don't need extra digits */
+			ln_dweight = 0;
+		}
+
+		free_var(&x);
+	}
+	else
+	{
+		/*
+		 * Estimate the logarithm using the first couple of digits from the
+		 * input number.  This will give an accurate result whenever the input
+		 * is not too close to 1.
+		 */
+		if (var->ndigits > 0)
+		{
+			int			digits;
+			int			dweight;
+			double		ln_var;
+
+			digits = var->digits[0];
+			dweight = var->weight * DEC_DIGITS;
+
+			if (var->ndigits > 1)
+			{
+				digits = digits * NBASE + var->digits[1];
+				dweight -= DEC_DIGITS;
+			}
+
+			/*----------
+			 * We have var ~= digits * 10^dweight
+			 * so ln(var) ~= ln(digits) + dweight * ln(10)
+			 *----------
+			 */
+			ln_var = log((double) digits) + dweight * 2.302585092994046;
+			ln_dweight = (int) log10(Abs(ln_var));
+		}
+		else
+		{
+			/* Caller should fail on ln(0), but for the moment return zero */
+			ln_dweight = 0;
+		}
+	}
+
+	return ln_dweight;
+}
+
+
+/*
  * ln_var() -
  *
  *	Compute the natural log of x
@@ -6814,8 +6867,6 @@ ln_var(NumericVar *arg, NumericVar *result, int rscale)
 				(errcode(ERRCODE_INVALID_ARGUMENT_FOR_LOG),
 				 errmsg("cannot take logarithm of a negative number")));
 
-	local_rscale = rscale + 8;
-
 	init_var(&x);
 	init_var(&xx);
 	init_var(&ni);
@@ -6825,16 +6876,25 @@ ln_var(NumericVar *arg, NumericVar *result, int rscale)
 	set_var_from_var(arg, &x);
 	set_var_from_var(&const_two, &fact);
 
-	/* Reduce input into range 0.9 < x < 1.1 */
+	/*
+	 * Reduce input into range 0.9 < x < 1.1 with repeated sqrt() operations.
+	 *
+	 * The final logarithm will have up to around rscale+6 significant digits.
+	 * Each sqrt() will roughly halve the weight of x, so adjust the local
+	 * rscale as we work so that we keep this many significant digits at each
+	 * step (plus a few more for good measure).
+	 */
 	while (cmp_var(&x, &const_zero_point_nine) <= 0)
 	{
-		local_rscale++;
+		local_rscale = rscale - x.weight * DEC_DIGITS / 2 + 8;
+		local_rscale = Max(local_rscale, NUMERIC_MIN_DISPLAY_SCALE);
 		sqrt_var(&x, &x, local_rscale);
 		mul_var(&fact, &const_two, &fact, 0);
 	}
 	while (cmp_var(&x, &const_one_point_one) >= 0)
 	{
-		local_rscale++;
+		local_rscale = rscale - x.weight * DEC_DIGITS / 2 + 8;
+		local_rscale = Max(local_rscale, NUMERIC_MIN_DISPLAY_SCALE);
 		sqrt_var(&x, &x, local_rscale);
 		mul_var(&fact, &const_two, &fact, 0);
 	}
@@ -6850,6 +6910,8 @@ ln_var(NumericVar *arg, NumericVar *result, int rscale)
 	 * The convergence of this is not as fast as one would like, but is
 	 * tolerable given that z is small.
 	 */
+	local_rscale = rscale + 8;
+
 	sub_var(&x, &const_one, result);
 	add_var(&x, &const_one, &elem);
 	div_var_fast(result, &elem, result, local_rscale, true);
@@ -6896,42 +6958,47 @@ log_var(NumericVar *base, NumericVar *num, NumericVar *result)
 {
 	NumericVar	ln_base;
 	NumericVar	ln_num;
-	int			dec_digits;
+	int			ln_base_dweight;
+	int			ln_num_dweight;
+	int			result_dweight;
 	int			rscale;
-	int			local_rscale;
+	int			ln_base_rscale;
+	int			ln_num_rscale;
 
 	init_var(&ln_base);
 	init_var(&ln_num);
 
-	/* Set scale for ln() calculations --- compare numeric_ln() */
-
-	/* Approx decimal digits before decimal point */
-	dec_digits = (num->weight + 1) * DEC_DIGITS;
-
-	if (dec_digits > 1)
-		rscale = NUMERIC_MIN_SIG_DIGITS - (int) log10(dec_digits - 1);
-	else if (dec_digits < 1)
-		rscale = NUMERIC_MIN_SIG_DIGITS - (int) log10(1 - dec_digits);
-	else
-		rscale = NUMERIC_MIN_SIG_DIGITS;
+	/* Estimated dweights of ln(base), ln(num) and the final result */
+	ln_base_dweight = estimate_ln_dweight(base);
+	ln_num_dweight = estimate_ln_dweight(num);
+	result_dweight = ln_num_dweight - ln_base_dweight;
 
+	/*
+	 * Select the scale of the result so that it will have at least
+	 * NUMERIC_MIN_SIG_DIGITS significant digits and is not less than either
+	 * input's display scale.
+	 */
+	rscale = NUMERIC_MIN_SIG_DIGITS - result_dweight;
 	rscale = Max(rscale, base->dscale);
 	rscale = Max(rscale, num->dscale);
 	rscale = Max(rscale, NUMERIC_MIN_DISPLAY_SCALE);
 	rscale = Min(rscale, NUMERIC_MAX_DISPLAY_SCALE);
 
-	local_rscale = rscale + 8;
-
-	/* Form natural logarithms */
-	ln_var(base, &ln_base, local_rscale);
-	ln_var(num, &ln_num, local_rscale);
+	/*
+	 * Set the scales for ln(base) and ln(num) so that they each have more
+	 * significant digits than the final result.
+	 */
+	ln_base_rscale = rscale + result_dweight - ln_base_dweight + 8;
+	ln_base_rscale = Max(ln_base_rscale, NUMERIC_MIN_DISPLAY_SCALE);
 
-	ln_base.dscale = rscale;
-	ln_num.dscale = rscale;
+	ln_num_rscale = rscale + result_dweight - ln_num_dweight + 8;
+	ln_num_rscale = Max(ln_num_rscale, NUMERIC_MIN_DISPLAY_SCALE);
 
-	/* Select scale for division result */
-	rscale = select_div_scale(&ln_num, &ln_base);
+	/* Form natural logarithms */
+	ln_var(base, &ln_base, ln_base_rscale);
+	ln_var(num, &ln_num, ln_num_rscale);
 
+	/* Divide and round to the required scale */
 	div_var_fast(&ln_num, &ln_base, result, rscale, true);
 
 	free_var(&ln_num);
@@ -6951,7 +7018,7 @@ power_var(NumericVar *base, NumericVar *exp, NumericVar *result)
 {
 	NumericVar	ln_base;
 	NumericVar	ln_num;
-	int			dec_digits;
+	int			ln_dweight;
 	int			rscale;
 	int			local_rscale;
 	double		val;
@@ -6982,7 +7049,7 @@ power_var(NumericVar *base, NumericVar *exp, NumericVar *result)
 	}
 
 	/*
-	 * This avoids log(0) for cases of 0 raised to a non-integer. 0 ^ 0
+	 * This avoids log(0) for cases of 0 raised to a non-integer.  0 ^ 0 is
 	 * handled by power_var_int().
 	 */
 	if (cmp_var(base, &const_zero) == 0)
@@ -6995,49 +7062,50 @@ power_var(NumericVar *base, NumericVar *exp, NumericVar *result)
 	init_var(&ln_base);
 	init_var(&ln_num);
 
-	/* Set scale for ln() calculation --- need extra accuracy here */
-
-	/* Approx decimal digits before decimal point */
-	dec_digits = (base->weight + 1) * DEC_DIGITS;
-
-	if (dec_digits > 1)
-		rscale = NUMERIC_MIN_SIG_DIGITS * 2 - (int) log10(dec_digits - 1);
-	else if (dec_digits < 1)
-		rscale = NUMERIC_MIN_SIG_DIGITS * 2 - (int) log10(1 - dec_digits);
-	else
-		rscale = NUMERIC_MIN_SIG_DIGITS * 2;
-
-	rscale = Max(rscale, base->dscale * 2);
-	rscale = Max(rscale, exp->dscale * 2);
-	rscale = Max(rscale, NUMERIC_MIN_DISPLAY_SCALE * 2);
-	rscale = Min(rscale, NUMERIC_MAX_DISPLAY_SCALE * 2);
+	/*----------
+	 * Decide on the scale for the ln() calculation.  For this we need an
+	 * estimate of the weight of the result, which we obtain by doing an
+	 * initial low-precision calculation of exp * ln(base).
+	 *
+	 * We want result = e ^ (exp * ln(base))
+	 * so result dweight = log10(result) = exp * ln(base) * log10(e)
+	 *----------
+	 */
+	ln_dweight = estimate_ln_dweight(base);
 
-	local_rscale = rscale + 8;
+	local_rscale = 8 - ln_dweight;
+	local_rscale = Max(local_rscale, NUMERIC_MIN_DISPLAY_SCALE);
+	local_rscale = Min(local_rscale, NUMERIC_MAX_DISPLAY_SCALE);
 
 	ln_var(base, &ln_base, local_rscale);
 
 	mul_var(&ln_base, exp, &ln_num, local_rscale);
 
-	/* Set scale for exp() -- compare numeric_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;
+	val *= 0.434294481903252;	/* approximate decimal result weight */
 
 	/* limit to something that won't cause integer overflow */
 	val = Max(val, -NUMERIC_MAX_RESULT_SCALE);
 	val = Min(val, NUMERIC_MAX_RESULT_SCALE);
 
+	/* choose the result scale */
 	rscale = NUMERIC_MIN_SIG_DIGITS - (int) val;
 	rscale = Max(rscale, base->dscale);
 	rscale = Max(rscale, exp->dscale);
 	rscale = Max(rscale, NUMERIC_MIN_DISPLAY_SCALE);
 	rscale = Min(rscale, NUMERIC_MAX_DISPLAY_SCALE);
 
+	/* set the scale for the real exp * ln(base) calculation */
+	local_rscale = rscale + (int) val - ln_dweight + 8;
+	local_rscale = Max(local_rscale, NUMERIC_MIN_DISPLAY_SCALE);
+
+	/* and do the real calculation */
+
+	ln_var(base, &ln_base, local_rscale);
+
+	mul_var(&ln_base, exp, &ln_num, local_rscale);
+
 	exp_var(&ln_num, result, rscale);
 
 	free_var(&ln_num);
@@ -7052,6 +7120,10 @@ power_var(NumericVar *base, NumericVar *exp, NumericVar *result)
 static void
 power_var_int(NumericVar *base, int exp, NumericVar *result, int rscale)
 {
+	double		f;
+	int			p;
+	int			i;
+	int			sig_digits;
 	unsigned int mask;
 	bool		neg;
 	NumericVar	base_prod;
@@ -7086,15 +7158,75 @@ power_var_int(NumericVar *base, int exp, NumericVar *result, int rscale)
 			break;
 	}
 
+	/* Handle the special case where the base is zero */
+	if (base->ndigits == 0)
+	{
+		if (exp < 0)
+			ereport(ERROR,
+					(errcode(ERRCODE_DIVISION_BY_ZERO),
+					 errmsg("division by zero")));
+		zero_var(result);
+		result->dscale = rscale;
+		return;
+	}
+
 	/*
 	 * The general case repeatedly multiplies base according to the bit
-	 * pattern of exp.  We do the multiplications with some extra precision.
+	 * pattern of exp.
+	 *
+	 * First we need to estimate the weight of the result so that we know how
+	 * many significant digits are needed.
+	 */
+	f = base->digits[0];
+	p = base->weight * DEC_DIGITS;
+
+	for (i = 1; i < base->ndigits && i * DEC_DIGITS < 16; i++)
+	{
+		f = f * NBASE + base->digits[i];
+		p -= DEC_DIGITS;
+	}
+
+	/*----------
+	 * We have base ~= f * 10^p
+	 * so log10(result) = log10(base^exp) ~= exp * (log10(f) + p)
+	 *----------
+	 */
+	f = exp * (log10(f) + p);
+
+	/*
+	 * Apply crude overflow/underflow tests so we can exit early if the result
+	 * certainly will overflow/underflow.
+	 */
+	if (f > 3 * SHRT_MAX * DEC_DIGITS)
+		ereport(ERROR,
+				(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
+				 errmsg("value overflows numeric format")));
+	if (f + 1 < -rscale || f + 1 < -NUMERIC_MAX_DISPLAY_SCALE)
+	{
+		zero_var(result);
+		result->dscale = rscale;
+		return;
+	}
+
+	/*
+	 * Approximate number of significant digits in the result.  Note that the
+	 * underflow test above means that this is necessarily >= 0.
+	 */
+	sig_digits = 1 + rscale + (int) f;
+
+	/*
+	 * The multiplications to produce the result may introduce an error of up
+	 * to around log10(abs(exp)) digits, so work with this many extra digits
+	 * of precision (plus a few more for good measure).
+	 */
+	sig_digits += (int) log(Abs(exp)) + 8;
+
+	/*
+	 * Now we can proceed with the multiplications.
 	 */
 	neg = (exp < 0);
 	mask = Abs(exp);
 
-	local_rscale = rscale + MUL_GUARD_DIGITS * 2;
-
 	init_var(&base_prod);
 	set_var_from_var(base, &base_prod);
 
@@ -7105,9 +7237,27 @@ power_var_int(NumericVar *base, int exp, NumericVar *result, int rscale)
 
 	while ((mask >>= 1) > 0)
 	{
+		/*
+		 * Do the multiplications using rscales large enough to hold the
+		 * results to the required number of significant digits, but don't
+		 * waste time by exceeding the scales of the numbers themselves.
+		 */
+		local_rscale = sig_digits - 2 * base_prod.weight * DEC_DIGITS;
+		local_rscale = Min(local_rscale, 2 * base_prod.dscale);
+		local_rscale = Max(local_rscale, NUMERIC_MIN_DISPLAY_SCALE);
+
 		mul_var(&base_prod, &base_prod, &base_prod, local_rscale);
+
 		if (mask & 1)
+		{
+			local_rscale = sig_digits -
+				(base_prod.weight + result->weight) * DEC_DIGITS;
+			local_rscale = Min(local_rscale,
+							   base_prod.dscale + result->dscale);
+			local_rscale = Max(local_rscale, NUMERIC_MIN_DISPLAY_SCALE);
+
 			mul_var(&base_prod, result, result, local_rscale);
+		}
 
 		/*
 		 * When abs(base) > 1, the number of digits to the left of the decimal