open-dsi
/
postgres
mirror of https://github.com/postgres/postgres
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Improve type numeric's calculations for ln(), log(), exp(), pow().

Browse Source 
Set the "rscales" for intermediate-result calculations to ensure that
suitable numbers of significant digits are maintained throughout.  The
previous coding hadn't thought this through in any detail, and as a result
could deliver results with many inaccurate digits, or in the worst cases
even fail with divide-by-zero errors as a result of losing all nonzero
digits of intermediate results.

In exp_var(), get rid entirely of the logic that separated the calculation
into integer and fractional parts: that was neither accurate nor
particularly fast.  The existing range-reduction method of dividing by 2^n
can be applied across the full input range instead of only 0..1, as long as
we are careful to set an appropriate rscale for each step.

Also fix the logic in mul_var() for shortening the calculation when the
caller asks for fewer output digits than an exact calculation would
require.  This bug doesn't affect simple multiplications since that code
path asks for an exact result, but it does contribute to accuracy issues
in the transcendental math functions.

In passing, improve performance of mul_var() a bit by forcing the shorter
input to be on the left, thus reducing the number of iterations of the
outer loop and probably also reducing the number of carry-propagation
steps needed.

This is arguably a bug fix, but in view of the lack of field complaints,
it does not seem worth the risk of back-patching.

Dean Rasheed
pull/14/head^2
Tom Lane 10 years ago
parent e57646e962
commit 7d9a4737c2
5 changed files with 2853 additions and 227 deletions
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Show Stats Download Patch File Download Diff File
  1. 604
      src/backend/utils/adt/numeric.c
  2. 302
      src/test/regress/expected/numeric.out
  3. 1376
      src/test/regress/expected/numeric_big.out
  4. 103
      src/test/regress/sql/numeric.sql
  5. 695
      src/test/regress/sql/numeric_big.sql

604
src/backend/utils/adt/numeric.c
Unescape View File

@ -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.
@ -360,20 +366,6 @@ static NumericDigit const_zero_point_nine_data[1] = {9};
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
@ -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,35 +6726,121 @@ 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() -
*
@ -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

302
src/test/regress/expected/numeric.out
Unescape View File

@ -1460,6 +1460,163 @@ select 10.0 ^ 2147483647 as overflows;
ERROR: value overflows numeric format
select 117743296169.0 ^ 1000000000 as overflows;
ERROR: value overflows numeric format
-- cases that used to return inaccurate results
select 3.789 ^ 21;
?column?
--------------------------------
1409343026052.8716016316022141
(1 row)
select 3.789 ^ 35;
?column?
----------------------------------------
177158169650516670809.3820586142670135
(1 row)
select 1.2 ^ 345;
?column?
-----------------------------------------------
2077446682327378559843444695.5827049735727869
(1 row)
select 0.12 ^ (-20);
?column?
--------------------------------------
2608405330458882702.5529619561355838
(1 row)
-- cases that used to error out
select 0.12 ^ (-25);
?column?
-------------------------------------------
104825960103961013959336.4983657883169110
(1 row)
select 0.5678 ^ (-85);
?column?
----------------------------------------
782333637740774446257.7719390061997396
(1 row)
--
-- Tests for raising to non-integer powers
--
-- special cases
select 0.0 ^ 0.0;
?column?
--------------------
1.0000000000000000
(1 row)
select (-12.34) ^ 0.0;
?column?
--------------------
1.0000000000000000
(1 row)
select 12.34 ^ 0.0;
?column?
--------------------
1.0000000000000000
(1 row)
select 0.0 ^ 12.34;
?column?
--------------------
0.0000000000000000
(1 row)
-- invalid inputs
select 0.0 ^ (-12.34);
ERROR: zero raised to a negative power is undefined
select (-12.34) ^ 1.2;
ERROR: a negative number raised to a non-integer power yields a complex result
-- cases that used to generate inaccurate results
select 32.1 ^ 9.8;
?column?
--------------------
580429286790711.10
(1 row)
select 32.1 ^ (-9.8);
?column?
----------------------------------
0.000000000000001722862754788209
(1 row)
select 12.3 ^ 45.6;
?column?
------------------------------------------------------
50081010321492803393171165777624533697036806969694.9
(1 row)
select 12.3 ^ (-45.6);
?column?
---------------------------------------------------------------------
0.00000000000000000000000000000000000000000000000001996764828785491
(1 row)
-- big test
select 1.234 ^ 5678;
?column?
----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
307239295662090741644584872593956173493568238595074141254349565406661439636598896798876823220904084953233015553994854875890890858118656468658643918169805277399402542281777901029346337707622181574346585989613344285010764501017625366742865066948856161360224801370482171458030533346309750557140549621313515752078638620714732831815297168231790779296290266207315344008883935010274044001522606235576584215999260117523114297033944018699691024106823438431754073086813382242140602291215149759520833200152654884259619588924545324.5973362312547382
(1 row)
--
-- Tests for EXP()
--
-- special cases
select exp(0.0);
exp
--------------------
1.0000000000000000
(1 row)
select exp(1.0);
exp
--------------------
2.7182818284590452
(1 row)
select exp(1.0::numeric(71,70));
exp
--------------------------------------------------------------------------
2.7182818284590452353602874713526624977572470936999595749669676277240766
(1 row)
-- cases that used to generate inaccurate results
select exp(32.999);
exp
---------------------
214429043492155.053
(1 row)
select exp(-32.999);
exp
----------------------------------
0.000000000000004663547361468248
(1 row)
select exp(123.456);
exp
------------------------------------------------------------
413294435277809344957685441227343146614594393746575438.725
(1 row)
select exp(-123.456);
exp
-------------------------------------------------------------------------
0.000000000000000000000000000000000000000000000000000002419582541264601
(1 row)
-- big test
select exp(1234.5678);
exp
----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
146549072930959479983482138503979804217622199675223653966270157446954995433819741094410764947112047906012815540251009949604426069672532417736057033099274204598385314594846509975629046864798765888104789074984927709616261452461385220475510438783429612447831614003668421849727379202555580791042606170523016207262965336641214601082882495255771621327088265411334088968112458492660609809762865582162764292604697957813514621259353683899630997077707406305730694385703091201347848855199354307506425820147289848677003277208302716466011827836279231.9667
(1 row)
--
-- Tests for generate_series
--
@ -1550,3 +1707,148 @@ select * from generate_series(1::numeric, 3::numeric) i, generate_series(1,5,i)
3 | 4
(10 rows)
--
-- Tests for LN()
--
-- Invalid inputs
select ln(-12.34);
ERROR: cannot take logarithm of a negative number
select ln(0.0);
ERROR: cannot take logarithm of zero
-- Some random tests
select ln(1.2345678e-28);
ln
-----------------------------------------
-64.26166165451762991204894255882820859
(1 row)
select ln(0.0456789);
ln
---------------------
-3.0861187944847439
(1 row)
select ln(0.349873948359354029493948309745709580730482050975);
ln
-----------------------------------------------------
-1.050182336912082775693991697979750253056317885460
(1 row)
select ln(0.99949452);
ln
-------------------------
-0.00050560779808326467
(1 row)
select ln(1.00049687395);
ln
------------------------
0.00049675054901370394
(1 row)
select ln(1234.567890123456789);
ln
--------------------
7.1184763012977896
(1 row)
select ln(5.80397490724e5);
ln
--------------------
13.271468476626518
(1 row)
select ln(9.342536355e34);
ln
--------------------
80.522470935524187
(1 row)
--
-- Tests for LOG() (base 10)
--
-- invalid inputs
select log(-12.34);
ERROR: cannot take logarithm of a negative number
CONTEXT: SQL function "log" statement 1
select log(0.0);
ERROR: cannot take logarithm of zero
CONTEXT: SQL function "log" statement 1
-- some random tests
select log(1.234567e-89);
log
-----------------------------------------------------------------------------------------------------
-88.90848533591373725637496492944925187293052336306443143312825869985819779294142441287021741054275
(1 row)
select log(3.4634998359873254962349856073435545);
log
--------------------------------------
0.5395151714070134409152404011959981
(1 row)
select log(9.999999999999999999);
log
----------------------
1.000000000000000000
(1 row)
select log(10.00000000000000000);
log
---------------------
1.00000000000000000
(1 row)
select log(10.00000000000000001);
log
---------------------
1.00000000000000000
(1 row)
select log(590489.45235237);
log
-------------------
5.771212144411727
(1 row)
--
-- Tests for LOG() (arbitrary base)
--
-- invalid inputs
select log(-12.34, 56.78);
ERROR: cannot take logarithm of a negative number
select log(-12.34, -56.78);
ERROR: cannot take logarithm of a negative number
select log(12.34, -56.78);
ERROR: cannot take logarithm of a negative number
select log(0.0, 12.34);
ERROR: cannot take logarithm of zero
select log(12.34, 0.0);
ERROR: cannot take logarithm of zero
select log(1.0, 12.34);
ERROR: division by zero
-- some random tests
select log(1.23e-89, 6.4689e45);
log
------------------------------------------------------------------------------------------------
-0.5152489207781856983977054971756484879653568168479201885425588841094788842469115325262329756
(1 row)
select log(0.99923, 4.58934e34);
log
---------------------
-103611.55579544132
(1 row)
select log(1.000016, 8.452010e18);
log
--------------------
2723830.2877097365
(1 row)
select log(3.1954752e47, 9.4792021e-73);
log
-------------------------------------------------------------------------------------
-1.51613372350688302142917386143459361608600157692779164475351842333265418126982165
(1 row)

1376
src/test/regress/expected/numeric_big.out
View File

File diff suppressed because it is too large Load Diff

103
src/test/regress/sql/numeric.sql
Unescape View File

@ -860,6 +860,57 @@ select 10.0 ^ -2147483647 as rounds_to_zero;
select 10.0 ^ 2147483647 as overflows;
select 117743296169.0 ^ 1000000000 as overflows;
-- cases that used to return inaccurate results
select 3.789 ^ 21;
select 3.789 ^ 35;
select 1.2 ^ 345;
select 0.12 ^ (-20);
-- cases that used to error out
select 0.12 ^ (-25);
select 0.5678 ^ (-85);
--
-- Tests for raising to non-integer powers
--
-- special cases
select 0.0 ^ 0.0;
select (-12.34) ^ 0.0;
select 12.34 ^ 0.0;
select 0.0 ^ 12.34;
-- invalid inputs
select 0.0 ^ (-12.34);
select (-12.34) ^ 1.2;
-- cases that used to generate inaccurate results
select 32.1 ^ 9.8;
select 32.1 ^ (-9.8);
select 12.3 ^ 45.6;
select 12.3 ^ (-45.6);
-- big test
select 1.234 ^ 5678;
--
-- Tests for EXP()
--
-- special cases
select exp(0.0);
select exp(1.0);
select exp(1.0::numeric(71,70));
-- cases that used to generate inaccurate results
select exp(32.999);
select exp(-32.999);
select exp(123.456);
select exp(-123.456);
-- big test
select exp(1234.5678);
--
-- Tests for generate_series
--
@ -880,3 +931,55 @@ select (i / (10::numeric ^ 131071))::numeric(1,0)
select * from generate_series(1::numeric, 3::numeric) i, generate_series(i,3) j;
select * from generate_series(1::numeric, 3::numeric) i, generate_series(1,i) j;
select * from generate_series(1::numeric, 3::numeric) i, generate_series(1,5,i) j;
--
-- Tests for LN()
--
-- Invalid inputs
select ln(-12.34);
select ln(0.0);
-- Some random tests
select ln(1.2345678e-28);
select ln(0.0456789);
select ln(0.349873948359354029493948309745709580730482050975);
select ln(0.99949452);
select ln(1.00049687395);
select ln(1234.567890123456789);
select ln(5.80397490724e5);
select ln(9.342536355e34);
--
-- Tests for LOG() (base 10)
--
-- invalid inputs
select log(-12.34);
select log(0.0);
-- some random tests
select log(1.234567e-89);
select log(3.4634998359873254962349856073435545);
select log(9.999999999999999999);
select log(10.00000000000000000);
select log(10.00000000000000001);
select log(590489.45235237);
--
-- Tests for LOG() (arbitrary base)
--
-- invalid inputs
select log(-12.34, 56.78);
select log(-12.34, -56.78);
select log(12.34, -56.78);
select log(0.0, 12.34);
select log(12.34, 0.0);
select log(1.0, 12.34);
-- some random tests
select log(1.23e-89, 6.4689e45);
select log(0.99923, 4.58934e34);
select log(1.000016, 8.452010e18);
select log(3.1954752e47, 9.4792021e-73);

695
src/test/regress/sql/numeric_big.sql
Unescape View File

@ -645,3 +645,698 @@ SELECT t1.id1, t1.result, t2.expected
FROM num_result t1, num_exp_power_10_ln t2
WHERE t1.id1 = t2.id
AND t1.result != t2.expected;
--
-- Test code path for raising to integer powers
--
-- base less than 1
--
-- bc(1) results computed with a scale of 500 and truncated using the script
-- below, and then rounded by hand to match the precision of POW():
--
-- for p in {-20..20}
-- do
-- b="0.084738"
-- r=$(bc -ql <<< "scale=500 ; $b^$p" | head -n 1)
-- echo "($b, $p, $r),"
-- done
WITH t(b, p, bc_result) AS (VALUES
(0.084738, -20, 2744326694304960114888.7859130502035257),
(0.084738, -19, 232548755422013710215.4459407000481464),
(0.084738, -18, 19705716436950597776.2364581230406798),
(0.084738, -17, 1669822999434319754.3627249884302211),
(0.084738, -16, 141497461326065387.3451885900696001),
(0.084738, -15, 11990211877848128.7928565907453178),
(0.084738, -14, 1016026574105094.7376490817865767),
(0.084738, -13, 86096059836517.5178789078924309),
(0.084738, -12, 7295607918426.8214300228969888),
(0.084738, -11, 618215223791.6519943372802450),
(0.084738, -10, 52386321633.6570066961524534),
(0.084738, -9, 4439112122.5928274334185666),
(0.084738, -8, 376161483.0442710110530225),
(0.084738, -7, 31875171.7502054369346110),
(0.084738, -6, 2701038.3037689083149651),
(0.084738, -5, 228880.5837847697527935),
(0.084738, -4, 19394.8829087538193122),
(0.084738, -3, 1643.4835879219811409),
(0.084738, -2, 139.2655122733328379),
(0.084738, -1, 11.8010809790176780),
(0.084738, 0, 1),
(0.084738, 1, .084738),
(0.084738, 2, .007180528644),
(0.084738, 3, .0006084636362353),
(0.084738, 4, .0000515599916073),
(0.084738, 5, .0000043690905688),
(0.084738, 6, .0000003702279966),
(0.084738, 7, .0000000313723800),
(0.084738, 8, .0000000026584327),
(0.084738, 9, .0000000002252703),
(0.084738, 10, .0000000000190890),
(0.084738, 11, .0000000000016176),
(0.084738, 12, .0000000000001371),
(0.084738, 13, .0000000000000116),
(0.084738, 14, .0000000000000010),
(0.084738, 15, .0000000000000001),
(0.084738, 16, .0000000000000000),
(0.084738, 17, .0000000000000000),
(0.084738, 18, .0000000000000000),
(0.084738, 19, .0000000000000000),
(0.084738, 20, .0000000000000000))
SELECT b, p, bc_result, b^p AS power, b^p - bc_result AS diff FROM t;
-- base greater than 1
--
-- bc(1) results computed with a scale of 500 and truncated using the script
-- below, and then rounded by hand to match the precision of POW():
--
-- for p in {-20..20}
-- do
-- b="37.821637"
-- r=$(bc -ql <<< "scale=500 ; $b^$p" | head -n 1)
-- echo "($b, $p, $r),"
-- done
WITH t(b, p, bc_result) AS (VALUES
(37.821637, -20, .0000000000000000),
(37.821637, -19, .0000000000000000),
(37.821637, -18, .0000000000000000),
(37.821637, -17, .0000000000000000),
(37.821637, -16, .0000000000000000),
(37.821637, -15, .0000000000000000),
(37.821637, -14, .0000000000000000),
(37.821637, -13, .0000000000000000),
(37.821637, -12, .0000000000000000),
(37.821637, -11, .0000000000000000),
(37.821637, -10, .0000000000000002),
(37.821637, -9, .0000000000000063),
(37.821637, -8, .0000000000002388),
(37.821637, -7, .0000000000090327),
(37.821637, -6, .0000000003416316),
(37.821637, -5, .0000000129210673),
(37.821637, -4, .0000004886959182),
(37.821637, -3, .0000184832796213),
(37.821637, -2, .0006990678924066),
(37.821637, -1, .0264398920649574),
(37.821637, 0, 1),
(37.821637, 1, 37.821637),
(37.821637, 2, 1430.476225359769),
(37.821637, 3, 54102.9525326873775219),
(37.821637, 4, 2046262.2313195326271135),
(37.821637, 5, 77392987.3197773940323425),
(37.821637, 6, 2927129472.7542235178972258),
(37.821637, 7, 110708828370.5116321107718772),
(37.821637, 8, 4187189119324.7924539711577286),
(37.821637, 9, 158366346921451.9852944363360812),
(37.821637, 10, 5989674486279224.5007355092228730),
(37.821637, 11, 226539294168214309.7083246628376531),
(37.821637, 12, 8568086950266418559.9938312759931069),
(37.821637, 13, 324059074417413536066.1494087598581043),
(37.821637, 14, 12256444679171401239980.3109258799733927),
(37.821637, 15, 463558801566202198479885.2069857662592280),
(37.821637, 16, 17532552720991931019508170.1002855156233684),
(37.821637, 17, 663109844696719094948877928.0672523682648687),
(37.821637, 18, 25079899837245684700124994552.6717306599041850),
(37.821637, 19, 948562867640665366544581398598.1275771806665398),
(37.821637, 20, 35876200451584291931921101974730.6901038166532866))
SELECT b, p, bc_result, b^p AS power, b^p - bc_result AS diff FROM t;
--
-- Tests for raising to non-integer powers
--
-- base less than 1
--
-- bc(1) results computed with a scale of 500 and truncated using the script
-- below, and then rounded by hand to match the precision of POW():
--
-- for n in {-20..20}
-- do
-- b="0.06933247"
-- p="$n.342987"
-- r=$(bc -ql <<< "scale=500 ; e($p*l($b))" | head -n 1)
-- echo "($b, $p, $r),"
-- done
WITH t(b, p, bc_result) AS (VALUES
(0.06933247, -20.342987, 379149253615977128356318.39406340),
(0.06933247, -19.342987, 26287354251852125772450.59436685),
(0.06933247, -18.342987, 1822567200045909954554.65766042),
(0.06933247, -17.342987, 126363085720167050546.86216560),
(0.06933247, -16.342987, 8761064849800910427.02880469),
(0.06933247, -15.342987, 607426265866876128.15466179),
(0.06933247, -14.342987, 42114363355427213.14899924),
(0.06933247, -13.342987, 2919892833909256.59283660),
(0.06933247, -12.342987, 202443382310228.51544515),
(0.06933247, -11.342987, 14035899730722.44924025),
(0.06933247, -10.342987, 973143597003.32229028),
(0.06933247, -9.342987, 67470449244.92493259),
(0.06933247, -8.342987, 4677892898.16028054),
(0.06933247, -7.342987, 324329869.02491071),
(0.06933247, -6.342987, 22486590.914273551),
(0.06933247, -5.342987, 1559050.8899661435),
(0.06933247, -4.342987, 108092.84905705095),
(0.06933247, -3.342987, 7494.3442144625131),
(0.06933247, -2.342987, 519.60139541889576),
(0.06933247, -1.342987, 36.025248159838727),
(0.06933247, 0.342987, .40036522320023350),
(0.06933247, 1.342987, .02775830982657349),
(0.06933247, 2.342987, .001924552183301612),
(0.06933247, 3.342987, .0001334339565121935),
(0.06933247, 4.342987, .000009251305786862961),
(0.06933247, 5.342987, .0000006414158809285026),
(0.06933247, 6.342987, .00000004447094732199898),
(0.06933247, 7.342987, .000000003083280621074075),
(0.06933247, 8.342987, .0000000002137714611621997),
(0.06933247, 9.342987, .00000000001482130341788437),
(0.06933247, 10.342987, .000000000001027597574581366),
(0.06933247, 11.342987, .00000000000007124587801173530),
(0.06933247, 12.342987, .000000000000004939652699872298),
(0.06933247, 13.342987, .0000000000000003424783226243151),
(0.06933247, 14.342987, .00000000000000002374486802900065),
(0.06933247, 15.342987, .000000000000000001646290350274646),
(0.06933247, 16.342987, .0000000000000000001141413763217064),
(0.06933247, 17.342987, .000000000000000000007913703549583420),
(0.06933247, 18.342987, .0000000000000000000005486766139403860),
(0.06933247, 19.342987, .00000000000000000000003804110487572339),
(0.06933247, 20.342987, .000000000000000000000002637483762562946))
SELECT b, p, bc_result, b^p AS power, b^p - bc_result AS diff FROM t;
-- base greater than 1
--
-- bc(1) results computed with a scale of 500 and truncated using the script
-- below, and then rounded by hand to match the precision of POW():
--
-- for n in {-20..20}
-- do
-- b="27.234987"
-- p="$n.230957"
-- r=$(bc -ql <<< "scale=500 ; e($p*l($b))" | head -n 1)
-- echo "($b, $p, $r),"
-- done
WITH t(b, p, bc_result) AS (VALUES
(27.234987, -20.230957, .000000000000000000000000000009247064512095633),
(27.234987, -19.230957, .0000000000000000000000000002518436817750859),
(27.234987, -18.230957, .000000000000000000000000006858959399176602),
(27.234987, -17.230957, .0000000000000000000000001868036700701026),
(27.234987, -16.230957, .000000000000000000000005087595525911532),
(27.234987, -15.230957, .0000000000000000000001385605980094587),
(27.234987, -14.230957, .000000000000000000003773696085499835),
(27.234987, -13.230957, .0000000000000000001027765638305389),
(27.234987, -12.230957, .000000000000000002799118379829397),
(27.234987, -11.230957, .00000000000000007623395268611469),
(27.234987, -10.230957, .000000000000002076230710364949),
(27.234987, -9.230957, .00000000000005654611640579014),
(27.234987, -8.230957, .000000000001540032745212181),
(27.234987, -7.230957, .00000000004194277179542807),
(27.234987, -6.230957, .000000001142310844592450),
(27.234987, -5.230957, .00000003111082100243440),
(27.234987, -4.230957, .0000008473028055606278),
(27.234987, -3.230957, .00002307628089450723),
(27.234987, -2.230957, .0006284822101702527),
(27.234987, -1.230957, .01711670482371810),
(27.234987, 0.230957, 2.1451253063142300),
(27.234987, 1.230957, 58.422459830839071),
(27.234987, 2.230957, 1591.1349340009243),
(27.234987, 3.230957, 43334.539242761031),
(27.234987, 4.230957, 1180215.6129275865),
(27.234987, 5.230957, 32143156.875279851),
(27.234987, 6.230957, 875418459.63720737),
(27.234987, 7.230957, 23842010367.779367),
(27.234987, 8.230957, 649336842420.336290),
(27.234987, 9.230957, 17684680461938.907402),
(27.234987, 10.230957, 481642042480060.137900),
(27.234987, 11.230957, 13117514765597885.614921),
(27.234987, 12.230957, 357255344113366461.949871),
(27.234987, 13.230957, 9729844652608062117.440722),
(27.234987, 14.230957, 264992192625800087863.690528),
(27.234987, 15.230957, 7217058921265161257566.469315),
(27.234987, 16.230957, 196556505898890690402726.443417),
(27.234987, 17.230957, 5353213882921711267539279.451015),
(27.234987, 18.230957, 145794710509592328389185797.837767),
(27.234987, 19.230957, 3970717045397510438979206144.696206),
(27.234987, 20.230957, 108142427112079606637962972621.121293))
SELECT b, p, bc_result, b^p AS power, b^p - bc_result AS diff FROM t;
--
-- Tests for EXP()
--
-- bc(1) results computed with a scale of 500 and truncated using the script
-- below, and then rounded by hand to match the precision of EXP():
--
-- for n in {-20..20}
-- do
-- x="$n.29837"
-- r=$(bc -ql <<< "scale=500 ; e($x)" | head -n 1)
-- echo "($x, $r),"
-- done
WITH t(x, bc_result) AS (VALUES
(-20.29837, .000000001529431101152222),
(-19.29837, .000000004157424770142192),
(-18.29837, .00000001130105220586304),
(-17.29837, .00000003071944485366452),
(-16.29837, .00000008350410872606600),
(-15.29837, .0000002269877013517336),
(-14.29837, .0000006170165438681061),
(-13.29837, .000001677224859055276),
(-12.29837, .000004559169856609741),
(-11.29837, .00001239310857408049),
(-10.29837, .00003368796183504298),
(-9.29837, .00009157337449401917),
(-8.29837, .0002489222398577673),
(-7.29837, .0006766408013046928),
(-6.29837, .001839300394580514),
(-5.29837, .004999736839665763),
(-4.29837, .01359069379834070),
(-3.29837, .03694333598818056),
(-2.29837, .1004223988993283),
(-1.29837, .2729763820983097),
(0.29837, 1.3476603299656679),
(1.29837, 3.6633205858807959),
(2.29837, 9.9579377804197108),
(3.29837, 27.068481317440698),
(4.29837, 73.579760889182206),
(5.29837, 200.01052696742555),
(6.29837, 543.68498095607070),
(7.29837, 1477.8890041389891),
(8.29837, 4017.3188244304487),
(9.29837, 10920.204759575742),
(10.29837, 29684.194161006717),
(11.29837, 80690.005580314652),
(12.29837, 219338.17590722828),
(13.29837, 596222.97785597218),
(14.29837, 1620702.0864156289),
(15.29837, 4405525.0308492653),
(16.29837, 11975458.636179032),
(17.29837, 32552671.598188404),
(18.29837, 88487335.673150406),
(19.29837, 240533516.60908059),
(20.29837, 653837887.33381570))
SELECT x, bc_result, exp(x), exp(x)-bc_result AS diff FROM t;
--
-- Tests for LN()
--
-- input very small
--
-- bc(1) results computed with a scale of 500 and truncated using the script
-- below, and then rounded by hand to match the precision of LN():
--
-- for p in {1..40}
-- do
-- l=$(bc -ql <<< "scale=500 ; l(10^-$p)" | head -n 1)
-- echo "('1.0e-$p', $l),"
-- done
WITH t(x, bc_result) AS (VALUES
('1.0e-1', -2.3025850929940457),
('1.0e-2', -4.6051701859880914),
('1.0e-3', -6.9077552789821371),
('1.0e-4', -9.2103403719761827),
('1.0e-5', -11.512925464970228),
('1.0e-6', -13.815510557964274),
('1.0e-7', -16.118095650958320),
('1.0e-8', -18.420680743952365),
('1.0e-9', -20.723265836946411),
('1.0e-10', -23.025850929940457),
('1.0e-11', -25.328436022934503),
('1.0e-12', -27.631021115928548),
('1.0e-13', -29.933606208922594),
('1.0e-14', -32.236191301916640),
('1.0e-15', -34.5387763949106853),
('1.0e-16', -36.84136148790473094),
('1.0e-17', -39.143946580898776628),
('1.0e-18', -41.4465316738928223123),
('1.0e-19', -43.74911676688686799634),
('1.0e-20', -46.051701859880913680360),
('1.0e-21', -48.3542869528749593643778),
('1.0e-22', -50.65687204586900504839581),
('1.0e-23', -52.959457138863050732413803),
('1.0e-24', -55.2620422318570964164317949),
('1.0e-25', -57.56462732485114210044978637),
('1.0e-26', -59.867212417845187784467777822),
('1.0e-27', -62.1697975108392334684857692765),
('1.0e-28', -64.47238260383327915250376073116),
('1.0e-29', -66.774967696827324836521752185847),
('1.0e-30', -69.0775527898213705205397436405309),
('1.0e-31', -71.38013788281541620455773509521529),
('1.0e-32', -73.682722975809461888575726549899655),
('1.0e-33', -75.9853080688035075725937180045840189),
('1.0e-34', -78.28789316179755325661170945926838306),
('1.0e-35', -80.590478254791598940629700913952747266),
('1.0e-36', -82.8930633477856446246476923686371114736),
('1.0e-37', -85.19564844077969030866568382332147568124),
('1.0e-38', -87.498233533773735992683675278005839888842),
('1.0e-39', -89.8008186267677816767016667326902040964430),
('1.0e-40', -92.10340371976182736071965818737456830404406))
SELECT x, bc_result, ln(x::numeric), ln(x::numeric)-bc_result AS diff FROM t;
-- input very close to but smaller than 1
--
-- bc(1) results computed with a scale of 500 and truncated using the script
-- below, and then rounded by hand to match the precision of LN():
--
-- for p in {1..40}
-- do
-- l=$(bc -ql <<< "scale=500 ; l(1-10^-$p)" | head -n 1)
-- echo "('1.0e-$p', $l),"
-- done
WITH t(x, bc_result) AS (VALUES
('1.0e-1', -.10536051565782630),
('1.0e-2', -.010050335853501441),
('1.0e-3', -.0010005003335835335),
('1.0e-4', -.00010000500033335834),
('1.0e-5', -.000010000050000333336),
('1.0e-6', -.0000010000005000003333),
('1.0e-7', -.00000010000000500000033),
('1.0e-8', -.000000010000000050000000),
('1.0e-9', -.0000000010000000005000000),
('1.0e-10', -.00000000010000000000500000),
('1.0e-11', -.000000000010000000000050000),
('1.0e-12', -.0000000000010000000000005000),
('1.0e-13', -.00000000000010000000000000500),
('1.0e-14', -.000000000000010000000000000050),
('1.0e-15', -.0000000000000010000000000000005),
('1.0e-16', -.00000000000000010000000000000001),
('1.0e-17', -.000000000000000010000000000000000),
('1.0e-18', -.0000000000000000010000000000000000),
('1.0e-19', -.00000000000000000010000000000000000),
('1.0e-20', -.000000000000000000010000000000000000),
('1.0e-21', -.0000000000000000000010000000000000000),
('1.0e-22', -.00000000000000000000010000000000000000),
('1.0e-23', -.000000000000000000000010000000000000000),
('1.0e-24', -.0000000000000000000000010000000000000000),
('1.0e-25', -.00000000000000000000000010000000000000000),
('1.0e-26', -.000000000000000000000000010000000000000000),
('1.0e-27', -.0000000000000000000000000010000000000000000),
('1.0e-28', -.00000000000000000000000000010000000000000000),
('1.0e-29', -.000000000000000000000000000010000000000000000),
('1.0e-30', -.0000000000000000000000000000010000000000000000),
('1.0e-31', -.00000000000000000000000000000010000000000000000),
('1.0e-32', -.000000000000000000000000000000010000000000000000),
('1.0e-33', -.0000000000000000000000000000000010000000000000000),
('1.0e-34', -.00000000000000000000000000000000010000000000000000),
('1.0e-35', -.000000000000000000000000000000000010000000000000000),
('1.0e-36', -.0000000000000000000000000000000000010000000000000000),
('1.0e-37', -.00000000000000000000000000000000000010000000000000000),
('1.0e-38', -.000000000000000000000000000000000000010000000000000000),
('1.0e-39', -.0000000000000000000000000000000000000010000000000000000),
('1.0e-40', -.00000000000000000000000000000000000000010000000000000000))
SELECT '1-'||x, bc_result, ln(1.0-x::numeric), ln(1.0-x::numeric)-bc_result AS diff FROM t;
-- input very close to but larger than 1
--
-- bc(1) results computed with a scale of 500 and truncated using the script
-- below, and then rounded by hand to match the precision of LN():
--
-- for p in {1..40}
-- do
-- l=$(bc -ql <<< "scale=500 ; l(1+10^-$p)" | head -n 1)
-- echo "('1.0e-$p', $l),"
-- done
WITH t(x, bc_result) AS (VALUES
('1.0e-1', .09531017980432486),
('1.0e-2', .009950330853168083),
('1.0e-3', .0009995003330835332),
('1.0e-4', .00009999500033330834),
('1.0e-5', .000009999950000333331),
('1.0e-6', .0000009999995000003333),
('1.0e-7', .00000009999999500000033),
('1.0e-8', .000000009999999950000000),
('1.0e-9', .0000000009999999995000000),
('1.0e-10', .00000000009999999999500000),
('1.0e-11', .000000000009999999999950000),
('1.0e-12', .0000000000009999999999995000),
('1.0e-13', .00000000000009999999999999500),
('1.0e-14', .000000000000009999999999999950),
('1.0e-15', .0000000000000009999999999999995),
('1.0e-16', .00000000000000010000000000000000),
('1.0e-17', .000000000000000010000000000000000),
('1.0e-18', .0000000000000000010000000000000000),
('1.0e-19', .00000000000000000010000000000000000),
('1.0e-20', .000000000000000000010000000000000000),
('1.0e-21', .0000000000000000000010000000000000000),
('1.0e-22', .00000000000000000000010000000000000000),
('1.0e-23', .000000000000000000000010000000000000000),
('1.0e-24', .0000000000000000000000010000000000000000),
('1.0e-25', .00000000000000000000000010000000000000000),
('1.0e-26', .000000000000000000000000010000000000000000),
('1.0e-27', .0000000000000000000000000010000000000000000),
('1.0e-28', .00000000000000000000000000010000000000000000),
('1.0e-29', .000000000000000000000000000010000000000000000),
('1.0e-30', .0000000000000000000000000000010000000000000000),
('1.0e-31', .00000000000000000000000000000010000000000000000),
('1.0e-32', .000000000000000000000000000000010000000000000000),
('1.0e-33', .0000000000000000000000000000000010000000000000000),
('1.0e-34', .00000000000000000000000000000000010000000000000000),
('1.0e-35', .000000000000000000000000000000000010000000000000000),
('1.0e-36', .0000000000000000000000000000000000010000000000000000),
('1.0e-37', .00000000000000000000000000000000000010000000000000000),
('1.0e-38', .000000000000000000000000000000000000010000000000000000),
('1.0e-39', .0000000000000000000000000000000000000010000000000000000),
('1.0e-40', .00000000000000000000000000000000000000010000000000000000))
SELECT '1+'||x, bc_result, ln(1.0+x::numeric), ln(1.0+x::numeric)-bc_result AS diff FROM t;
-- input very large
--
-- bc(1) results computed with a scale of 500 and truncated using the script
-- below, and then rounded by hand to match the precision of LN():
--
-- for p in {1..40}
-- do
-- l=$(bc -ql <<< "scale=500 ; l(10^$p)" | head -n 1)
-- echo "('1.0e$p', $l),"
-- done
WITH t(x, bc_result) AS (VALUES
('1.0e1', 2.3025850929940457),
('1.0e2', 4.6051701859880914),
('1.0e3', 6.9077552789821371),
('1.0e4', 9.2103403719761827),
('1.0e5', 11.512925464970228),
('1.0e6', 13.815510557964274),
('1.0e7', 16.118095650958320),
('1.0e8', 18.420680743952365),
('1.0e9', 20.723265836946411),
('1.0e10', 23.025850929940457),
('1.0e11', 25.328436022934503),
('1.0e12', 27.631021115928548),
('1.0e13', 29.933606208922594),
('1.0e14', 32.236191301916640),
('1.0e15', 34.538776394910685),
('1.0e16', 36.841361487904731),
('1.0e17', 39.143946580898777),
('1.0e18', 41.446531673892822),
('1.0e19', 43.749116766886868),
('1.0e20', 46.051701859880914),
('1.0e21', 48.354286952874959),
('1.0e22', 50.656872045869005),
('1.0e23', 52.959457138863051),
('1.0e24', 55.262042231857096),
('1.0e25', 57.564627324851142),
('1.0e26', 59.867212417845188),
('1.0e27', 62.169797510839233),
('1.0e28', 64.472382603833279),
('1.0e29', 66.774967696827325),
('1.0e30', 69.077552789821371),
('1.0e31', 71.380137882815416),
('1.0e32', 73.682722975809462),
('1.0e33', 75.985308068803508),
('1.0e34', 78.287893161797553),
('1.0e35', 80.590478254791599),
('1.0e36', 82.893063347785645),
('1.0e37', 85.195648440779690),
('1.0e38', 87.498233533773736),
('1.0e39', 89.800818626767782),
('1.0e40', 92.103403719761827))
SELECT x, bc_result, ln(x::numeric), ln(x::numeric)-bc_result AS diff FROM t;
-- input huge
--
-- bc(1) results computed with a scale of 1000 and truncated using the script
-- below, and then rounded by hand to match the precision of LN():
--
-- for p in {1..10}
-- do
-- l=$(bc -ql <<< "scale=1000 ; l(10^${p}00)" | head -n 1)
-- echo "('1.0e${p}00', $l),"
-- done
WITH t(x, bc_result) AS (VALUES
('1.0e100', 230.25850929940457),
('1.0e200', 460.51701859880914),
('1.0e300', 690.77552789821371),
('1.0e400', 921.03403719761827),
('1.0e500', 1151.2925464970228),
('1.0e600', 1381.5510557964274),
('1.0e700', 1611.8095650958320),
('1.0e800', 1842.0680743952365),
('1.0e900', 2072.3265836946411),
('1.0e1000', 2302.5850929940457))
SELECT x, bc_result, ln(x::numeric), ln(x::numeric)-bc_result AS diff FROM t;
--
-- Tests for LOG() (base 10)
--
-- input very small, exact result known
WITH t(x) AS (SELECT '1e-'||n FROM generate_series(1, 100) g(n))
SELECT x, log(x::numeric) FROM t;
-- input very small, non-exact results
--
-- bc(1) results computed with a scale of 500 and truncated using the script
-- below, and then rounded by hand to match the precision of LN():
--
-- for p in {1..50..7}
-- do
-- for d in {9..1..3}
-- do
-- l=$(bc -ql <<< "scale=500 ; l($d*10^-$p) / l(10)" | head -n 1)
-- echo "('${d}.0e-$p', $l),"
-- done
-- done
WITH t(x, bc_result) AS (VALUES
('9.0e-1', -.04575749056067513),
('6.0e-1', -.2218487496163564),
('3.0e-1', -.5228787452803376),
('9.0e-8', -7.045757490560675),
('6.0e-8', -7.221848749616356),
('3.0e-8', -7.522878745280338),
('9.0e-15', -14.0457574905606751),
('6.0e-15', -14.2218487496163564),
('3.0e-15', -14.5228787452803376),
('9.0e-22', -21.04575749056067512540994),
('6.0e-22', -21.22184874961635636749123),
('3.0e-22', -21.52287874528033756270497),
('9.0e-29', -28.045757490560675125409944193490),
('6.0e-29', -28.221848749616356367491233202020),
('3.0e-29', -28.522878745280337562704972096745),
('9.0e-36', -35.0457574905606751254099441934897693816),
('6.0e-36', -35.2218487496163563674912332020203916640),
('3.0e-36', -35.5228787452803375627049720967448846908),
('9.0e-43', -42.04575749056067512540994419348976938159974227),
('6.0e-43', -42.22184874961635636749123320202039166403168125),
('3.0e-43', -42.52287874528033756270497209674488469079987114),
('9.0e-50', -49.045757490560675125409944193489769381599742271618608),
('6.0e-50', -49.221848749616356367491233202020391664031681254347196),
('3.0e-50', -49.522878745280337562704972096744884690799871135809304))
SELECT x, bc_result, log(x::numeric), log(x::numeric)-bc_result AS diff FROM t;
-- input very close to but smaller than 1
--
-- bc(1) results computed with a scale of 500 and truncated using the script
-- below, and then rounded by hand to match the precision of LN():
--
-- for p in {1..40..7}
-- do
-- for d in {9..1..3}
-- do
-- l=$(bc -ql <<< "scale=500 ; l(1-$d*10^-$p) / l(10)" | head -n 1)
-- echo "('${d}.0e-$p', $l),"
-- done
-- done
WITH t(x, bc_result) AS (VALUES
('9.0e-1', -1.0000000000000000),
('6.0e-1', -.3979400086720376),
('3.0e-1', -.1549019599857432),
('9.0e-8', -.000000039086505130185422),
('6.0e-8', -.000000026057669695925208),
('3.0e-8', -.000000013028834652530076),
('9.0e-15', -.0000000000000039086503371292840),
('6.0e-15', -.0000000000000026057668914195188),
('3.0e-15', -.0000000000000013028834457097574),
('9.0e-22', -.00000000000000000000039086503371292664),
('6.0e-22', -.00000000000000000000026057668914195110),
('3.0e-22', -.00000000000000000000013028834457097555),
('9.0e-29', -.000000000000000000000000000039086503371292664),
('6.0e-29', -.000000000000000000000000000026057668914195110),
('3.0e-29', -.000000000000000000000000000013028834457097555),
('9.0e-36', -.0000000000000000000000000000000000039086503371292664),
('6.0e-36', -.0000000000000000000000000000000000026057668914195110),
('3.0e-36', -.0000000000000000000000000000000000013028834457097555))
SELECT '1-'||x, bc_result, log(1.0-x::numeric), log(1.0-x::numeric)-bc_result AS diff FROM t;
-- input very close to but larger than 1
--
-- bc(1) results computed with a scale of 500 and truncated using the script
-- below, and then rounded by hand to match the precision of LN():
--
-- for p in {1..40..7}
-- do
-- for d in {9..1..3}
-- do
-- l=$(bc -ql <<< "scale=500 ; l(1+$d*10^-$p) / l(10)" | head -n 1)
-- echo "('${d}.0e-$p', $l),"
-- done
-- done
WITH t(x, bc_result) AS (VALUES
('9.0e-1', .2787536009528290),
('6.0e-1', .2041199826559248),
('3.0e-1', .1139433523068368),
('9.0e-8', .000000039086501612400118),
('6.0e-8', .000000026057668132465074),
('3.0e-8', .000000013028834261665042),
('9.0e-15', .0000000000000039086503371292489),
('6.0e-15', .0000000000000026057668914195031),
('3.0e-15', .0000000000000013028834457097535),
('9.0e-22', .00000000000000000000039086503371292664),
('6.0e-22', .00000000000000000000026057668914195110),
('3.0e-22', .00000000000000000000013028834457097555),
('9.0e-29', .000000000000000000000000000039086503371292664),
('6.0e-29', .000000000000000000000000000026057668914195110),
('3.0e-29', .000000000000000000000000000013028834457097555),
('9.0e-36', .0000000000000000000000000000000000039086503371292664),
('6.0e-36', .0000000000000000000000000000000000026057668914195110),
('3.0e-36', .0000000000000000000000000000000000013028834457097555))
SELECT '1+'||x, bc_result, log(1.0+x::numeric), log(1.0+x::numeric)-bc_result AS diff FROM t;
-- input very large, exact result known
WITH t(x) AS (SELECT '1e'||n FROM generate_series(1, 100) g(n))
SELECT x, log(x::numeric) FROM t;
-- input very large, non-exact results
--
-- bc(1) results computed with a scale of 500 and truncated using the script
-- below, and then rounded by hand to match the precision of LN():
--
-- for p in {10..50..7}
-- do
-- for d in {2..9..3}
-- do
-- l=$(bc -ql <<< "scale=500 ; l($d*10^$p) / l(10)" | head -n 1)
-- echo "('${d}.0e$p', $l),"
-- done
-- done
WITH t(x, bc_result) AS (VALUES
('2.0e10', 10.301029995663981),
('5.0e10', 10.698970004336019),
('8.0e10', 10.903089986991944),
('2.0e17', 17.301029995663981),
('5.0e17', 17.698970004336019),
('8.0e17', 17.903089986991944),
('2.0e24', 24.301029995663981),
('5.0e24', 24.698970004336019),
('8.0e24', 24.903089986991944),
('2.0e31', 31.301029995663981),
('5.0e31', 31.698970004336019),
('8.0e31', 31.903089986991944),
('2.0e38', 38.301029995663981),
('5.0e38', 38.698970004336019),
('8.0e38', 38.903089986991944),
('2.0e45', 45.30102999566398),
('5.0e45', 45.69897000433602),
('8.0e45', 45.90308998699194))
SELECT x, bc_result, log(x::numeric), log(x::numeric)-bc_result AS diff FROM t;

Write Preview
Loading…
Cancel
Save