Improve the accuracy of numeric power() for integer exponents.

This makes the choice of result scale of numeric power() for integer
exponents consistent with the choice for non-integer exponents, and
with the result scale of other numeric functions. Specifically, the
result scale will be at least as large as the scale of either input,
and sufficient to ensure that the result has at least 16 significant
digits.

Formerly, the result scale was based only on the scale of the first
input, without taking into account the weight of the result. For
results with negative weight, that could lead to results with very few
or even no non-zero significant digits (e.g., 10.0 ^ (-18) produced
0.0000000000000000).

Fix this by moving responsibility for the choice of result scale into
power_var_int(), which already has code to estimate the result weight.

Per report by Adrian Klaver and suggested fix by Tom Lane.

No back-patch -- arguably this is a bug fix, but one which is easy to
work around, so it doesn't seem worth the risk of changing query
results in stable branches.

Discussion: https://postgr.es/m/12a40226-70ac-3a3b-3d3a-fdaf9e32d312%40aklaver.com

Author: deanrasheed

src/backend/utils/adt/numeric.c

@@ -571,8 +571,8 @@ static void log_var(const NumericVar *base, const NumericVar *num,
 					NumericVar *result);
 static void power_var(const NumericVar *base, const NumericVar *exp,
 					  NumericVar *result);
-static void power_var_int(const NumericVar *base, int exp, NumericVar *result,
-						  int rscale);
+static void power_var_int(const NumericVar *base, int exp, int exp_dscale,
+						  NumericVar *result);
 static void power_ten_int(int exp, NumericVar *result);
 
 static int	cmp_abs(const NumericVar *var1, const NumericVar *var2);
@@ -10335,13 +10335,8 @@ power_var(const NumericVar *base, const NumericVar *exp, NumericVar *result)
 		{
 			if (expval64 >= PG_INT32_MIN && expval64 <= PG_INT32_MAX)
 			{
-				/* Okay, select rscale */
-				rscale = NUMERIC_MIN_SIG_DIGITS;
-				rscale = Max(rscale, base->dscale);
-				rscale = Max(rscale, NUMERIC_MIN_DISPLAY_SCALE);
-				rscale = Min(rscale, NUMERIC_MAX_DISPLAY_SCALE);
-
-				power_var_int(base, (int) expval64, result, rscale);
+				/* Okay, use power_var_int */
+				power_var_int(base, (int) expval64, exp->dscale, result);
 				return;
 			}
 		}
@@ -10475,19 +10470,76 @@ power_var(const NumericVar *base, const NumericVar *exp, NumericVar *result)
  * power_var_int() -
  *
  *	Raise base to the power of exp, where exp is an integer.
+ *
+ *	Note: this routine chooses dscale of the result.
  */
 static void
-power_var_int(const NumericVar *base, int exp, NumericVar *result, int rscale)
+power_var_int(const NumericVar *base, int exp, int exp_dscale,
+			  NumericVar *result)
 {
 	double		f;
 	int			p;
 	int			i;
+	int			rscale;
 	int			sig_digits;
 	unsigned int mask;
 	bool		neg;
 	NumericVar	base_prod;
 	int			local_rscale;
 
+	/*
+	 * Choose the result scale.  For this we need an estimate of the decimal
+	 * weight of the result, which we obtain by approximating using double
+	 * precision arithmetic.
+	 *
+	 * We also perform crude overflow/underflow tests here so that we can exit
+	 * early if the result is sure to overflow/underflow, and to guard against
+	 * integer overflow when choosing the result scale.
+	 */
+	if (base->ndigits != 0)
+	{
+		/*----------
+		 * Choose f (double) and p (int) such that base ~= f * 10^p.
+		 * Then log10(result) = log10(base^exp) ~= exp * (log10(f) + p).
+		 *----------
+		 */
+		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;
+		}
+
+		f = exp * (log10(f) + p);	/* approximate decimal result weight */
+	}
+	else
+		f = 0;					/* result is 0 or 1 (weight 0), or error */
+
+	/* overflow/underflow tests with fuzz factors */
+	if (f > (SHRT_MAX + 1) * DEC_DIGITS)
+		ereport(ERROR,
+				(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
+				 errmsg("value overflows numeric format")));
+	if (f + 1 < -NUMERIC_MAX_DISPLAY_SCALE)
+	{
+		zero_var(result);
+		result->dscale = NUMERIC_MAX_DISPLAY_SCALE;
+		return;
+	}
+
+	/*
+	 * Choose the result scale in the same way as power_var(), so it has at
+	 * least NUMERIC_MIN_SIG_DIGITS significant digits and is not less than
+	 * either input's display scale.
+	 */
+	rscale = NUMERIC_MIN_SIG_DIGITS - (int) f;
+	rscale = Max(rscale, base->dscale);
+	rscale = Max(rscale, exp_dscale);
+	rscale = Max(rscale, NUMERIC_MIN_DISPLAY_SCALE);
+	rscale = Min(rscale, NUMERIC_MAX_DISPLAY_SCALE);
+
 	/* Handle some common special cases, as well as corner cases */
 	switch (exp)
 	{
@@ -10532,43 +10584,15 @@ power_var_int(const NumericVar *base, int exp, NumericVar *result, int rscale)
 	 * The general case repeatedly multiplies base according to the bit
 	 * pattern of exp.
 	 *
-	 * First we need to estimate the weight of the result so that we know how
-	 * many significant digits are needed.
+	 * The local rscale used for each multiplication is varied to keep a fixed
+	 * number of significant digits, sufficient to give the required result
+	 * scale.
 	 */
-	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.
+	 * underflow test above, together with the choice of rscale, ensures that
+	 * this approximation is necessarily > 0.
 	 */
 	sig_digits = 1 + rscale + (int) f;