-/* ----------
+/*-------------------------------------------------------------------------
+ *
* numeric.c
+ * An exact numeric data type for the Postgres database system
*
- * An exact numeric data type for the Postgres database system
+ * Original coding 1998, Jan Wieck. Heavily revised 2003, Tom Lane.
*
- * 1998 Jan Wieck
+ * Many of the algorithmic ideas are borrowed from David M. Smith's "FM"
+ * multiple-precision math library, most recently published as Algorithm
+ * 786: Multiple-Precision Complex Arithmetic and Functions, ACM
+ * Transactions on Mathematical Software, Vol. 24, No. 4, December 1998,
+ * pages 359-367.
*
- * $Header: /cvsroot/pgsql/src/backend/utils/adt/numeric.c,v 1.58 2003/03/14 00:15:32 tgl Exp $
+ * Copyright (c) 1998-2003, PostgreSQL Global Development Group
*
- * ----------
+ * IDENTIFICATION
+ * $Header: /cvsroot/pgsql/src/backend/utils/adt/numeric.c,v 1.59 2003/03/21 01:58:04 tgl Exp $
+ *
+ *-------------------------------------------------------------------------
*/
#include "postgres.h"
#include <ctype.h>
#include <float.h>
+#include <limits.h>
#include <math.h>
#include <errno.h>
/* ----------
* Local data types
*
- * Note: the first digit of a NumericVar's value is assumed to be multiplied
- * by 10 ** weight. Another way to say it is that there are weight+1 digits
- * before the decimal point. It is possible to have weight < 0.
- *
+ * Numeric values are represented in a base-NBASE floating point format.
+ * Each "digit" ranges from 0 to NBASE-1. The type NumericDigit is signed
+ * and wide enough to store a digit. We assume that NBASE*NBASE can fit in
+ * an int. Although the purely calculational routines could handle any even
+ * NBASE that's less than sqrt(INT_MAX), in practice we are only interested
+ * in NBASE a power of ten, so that I/O conversions and decimal rounding
+ * are easy. Also, it's actually more efficient if NBASE is rather less than
+ * sqrt(INT_MAX), so that there is "headroom" for mul_var and div_var to
+ * postpone processing carries.
+ * ----------
+ */
+
+#if 0
+#define NBASE 10
+#define HALF_NBASE 5
+#define DEC_DIGITS 1 /* decimal digits per NBASE digit */
+#define MUL_GUARD_DIGITS 4 /* these are measured in NBASE digits */
+#define DIV_GUARD_DIGITS 8
+
+typedef signed char NumericDigit;
+#endif
+
+#if 0
+#define NBASE 100
+#define HALF_NBASE 50
+#define DEC_DIGITS 2 /* decimal digits per NBASE digit */
+#define MUL_GUARD_DIGITS 3 /* these are measured in NBASE digits */
+#define DIV_GUARD_DIGITS 6
+
+typedef signed char NumericDigit;
+#endif
+
+#if 1
+#define NBASE 10000
+#define HALF_NBASE 5000
+#define DEC_DIGITS 4 /* decimal digits per NBASE digit */
+#define MUL_GUARD_DIGITS 2 /* these are measured in NBASE digits */
+#define DIV_GUARD_DIGITS 4
+
+typedef int16 NumericDigit;
+#endif
+
+
+/* ----------
* The value represented by a NumericVar is determined by the sign, weight,
- * ndigits, and digits[] array. The rscale and dscale are carried along,
- * but they are just auxiliary information until rounding is done before
- * final storage or display. (Scales are the number of digits wanted
- * *after* the decimal point. Scales are always >= 0.)
+ * 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.
*
* buf points at the physical start of the palloc'd digit buffer for the
* NumericVar. digits points at the first digit in actual use (the one
- * with the specified weight). We normally leave an unused byte or two
+ * with the specified weight). We normally leave an unused digit or two
* (preset to zeroes) between buf and digits, so that there is room to store
* a carry out of the top digit without special pushups. We just need to
* decrement digits (and increment weight) to make room for the carry digit.
+ * (There is no such extra space in a numeric value stored in the database,
+ * only in a NumericVar in memory.)
*
* If buf is NULL then the digit buffer isn't actually palloc'd and should
* not be freed --- see the constants below for an example.
*
+ * dscale, or display scale, is the nominal precision expressed as number
+ * of digits after the decimal point (it must always be >= 0 at present).
+ * dscale may be more than the number of physically stored fractional digits,
+ * implying that we have suppressed storage of significant trailing zeroes.
+ * It should never be less than the number of stored digits, since that would
+ * imply hiding digits that are present. NOTE that dscale is always expressed
+ * in *decimal* digits, and so it may correspond to a fractional number of
+ * base-NBASE digits --- divide by DEC_DIGITS to convert to NBASE digits.
+ *
+ * rscale, or result scale, is the target precision for a computation.
+ * Like dscale it is expressed as number of *decimal* digits after the decimal
+ * point, and is always >= 0 at present.
+ * Note that rscale is not stored in variables --- it's figured on-the-fly
+ * from the dscales of the inputs.
+ *
* 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.
* ----------
*/
-typedef unsigned char NumericDigit;
-
typedef struct NumericVar
{
- int ndigits; /* number of digits in digits[] - can be
- * 0! */
+ int ndigits; /* # of digits in digits[] - can be 0! */
int weight; /* weight of first digit */
- int rscale; /* result scale */
- int dscale; /* display scale */
int sign; /* NUMERIC_POS, NUMERIC_NEG, or
* NUMERIC_NAN */
+ int dscale; /* display scale */
NumericDigit *buf; /* start of palloc'd space for digits[] */
- NumericDigit *digits; /* decimal digits */
+ NumericDigit *digits; /* base-NBASE digits */
} NumericVar;
/* ----------
- * Local data
- * ----------
- */
-static int global_rscale = 0;
-
-/* ----------
- * Some preinitialized variables we need often
+ * Some preinitialized constants
* ----------
*/
static NumericDigit const_zero_data[1] = {0};
static NumericVar const_zero =
-{0, 0, 0, 0, NUMERIC_POS, NULL, const_zero_data};
+{0, 0, NUMERIC_POS, 0, NULL, const_zero_data};
static NumericDigit const_one_data[1] = {1};
static NumericVar const_one =
-{1, 0, 0, 0, NUMERIC_POS, NULL, const_one_data};
+{1, 0, NUMERIC_POS, 0, NULL, const_one_data};
static NumericDigit const_two_data[1] = {2};
static NumericVar const_two =
-{1, 0, 0, 0, NUMERIC_POS, NULL, const_two_data};
-
-static NumericDigit const_zero_point_one_data[1] = {1};
-static NumericVar const_zero_point_one =
-{1, -1, 1, 1, NUMERIC_POS, NULL, const_zero_point_one_data};
-
+{1, 0, NUMERIC_POS, 0, NULL, const_two_data};
+
+#if DEC_DIGITS == 4
+static NumericDigit const_zero_point_five_data[1] = {5000};
+#elif DEC_DIGITS == 2
+static NumericDigit const_zero_point_five_data[1] = {50};
+#elif DEC_DIGITS == 1
+static NumericDigit const_zero_point_five_data[1] = {5};
+#endif
+static NumericVar const_zero_point_five =
+{1, -1, NUMERIC_POS, 1, NULL, const_zero_point_five_data};
+
+#if DEC_DIGITS == 4
+static NumericDigit const_zero_point_nine_data[1] = {9000};
+#elif DEC_DIGITS == 2
+static NumericDigit const_zero_point_nine_data[1] = {90};
+#elif DEC_DIGITS == 1
static NumericDigit const_zero_point_nine_data[1] = {9};
+#endif
static NumericVar const_zero_point_nine =
-{1, -1, 1, 1, NUMERIC_POS, NULL, const_zero_point_nine_data};
+{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};
+#elif DEC_DIGITS == 1
static NumericDigit const_one_point_one_data[2] = {1, 1};
+#endif
static NumericVar const_one_point_one =
-{2, 0, 1, 1, NUMERIC_POS, NULL, const_one_point_one_data};
+{2, 0, NUMERIC_POS, 1, NULL, const_one_point_one_data};
static NumericVar const_nan =
-{0, 0, 0, 0, NUMERIC_NAN, NULL, NULL};
+{0, 0, NUMERIC_NAN, 0, NULL, NULL};
+#if DEC_DIGITS == 4
+static const int round_powers[4] = { 0, 1000, 100, 10 };
+#endif
/* ----------
*/
#ifdef NUMERIC_DEBUG
-static void dump_numeric(char *str, Numeric num);
-static void dump_var(char *str, NumericVar *var);
+static void dump_numeric(const char *str, Numeric num);
+static void dump_var(const char *str, NumericVar *var);
#else
#define dump_numeric(s,n)
#define dump_var(s,v)
#endif
-#define digitbuf_alloc(size) ((NumericDigit *) palloc(size))
+#define digitbuf_alloc(ndigits) \
+ ((NumericDigit *) palloc((ndigits) * sizeof(NumericDigit)))
#define digitbuf_free(buf) \
do { \
if ((buf) != NULL) \
pfree(buf); \
} while (0)
-#define init_var(v) memset(v,0,sizeof(NumericVar))
+#define init_var(v) MemSetAligned(v, 0, sizeof(NumericVar))
+
static void alloc_var(NumericVar *var, int ndigits);
static void free_var(NumericVar *var);
static void zero_var(NumericVar *var);
-static void set_var_from_str(char *str, NumericVar *dest);
+static void set_var_from_str(const char *str, NumericVar *dest);
static void set_var_from_num(Numeric value, NumericVar *dest);
static void set_var_from_var(NumericVar *value, NumericVar *dest);
static char *get_str_from_var(NumericVar *var, int dscale);
static void apply_typmod(NumericVar *var, int32 typmod);
+static bool numericvar_to_int8(NumericVar *var, int64 *result);
+static void int8_to_numericvar(int64 val, NumericVar *var);
static double numeric_to_double_no_overflow(Numeric num);
static double numericvar_to_double_no_overflow(NumericVar *var);
static int cmp_var(NumericVar *var1, NumericVar *var2);
static void add_var(NumericVar *var1, NumericVar *var2, NumericVar *result);
static void sub_var(NumericVar *var1, NumericVar *var2, NumericVar *result);
-static void mul_var(NumericVar *var1, NumericVar *var2, NumericVar *result);
-static void div_var(NumericVar *var1, NumericVar *var2, NumericVar *result);
+static void mul_var(NumericVar *var1, NumericVar *var2, NumericVar *result,
+ int rscale);
+static void div_var(NumericVar *var1, NumericVar *var2, NumericVar *result,
+ int rscale);
static int select_div_scale(NumericVar *var1, NumericVar *var2);
static void mod_var(NumericVar *var1, NumericVar *var2, NumericVar *result);
static void ceil_var(NumericVar *var, NumericVar *result);
static void floor_var(NumericVar *var, NumericVar *result);
-static void sqrt_var(NumericVar *arg, NumericVar *result);
-static void exp_var(NumericVar *arg, NumericVar *result);
-static void ln_var(NumericVar *arg, 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 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);
+static void power_var_int(NumericVar *base, int exp, NumericVar *result,
+ int rscale);
static int cmp_abs(NumericVar *var1, NumericVar *var2);
static void add_abs(NumericVar *var1, NumericVar *var2, NumericVar *result);
static void sub_abs(NumericVar *var1, NumericVar *var2, NumericVar *result);
-
+static void round_var(NumericVar *var, int rscale);
+static void trunc_var(NumericVar *var, int rscale);
+static void strip_var(NumericVar *var);
/* ----------------------------------------------------------------------
*/
-/* ----------
+/*
* numeric_in() -
*
* Input function for numeric data type
- * ----------
*/
Datum
numeric_in(PG_FUNCTION_ARGS)
}
-/* ----------
+/*
* numeric_out() -
*
* Output function for numeric data type
- * ----------
*/
Datum
numeric_out(PG_FUNCTION_ARGS)
}
-/* ----------
+/*
* numeric() -
*
* This is a special function called by the Postgres database system
- * before a value is stored in a tuples attribute. The precision and
+ * before a value is stored in a tuple's attribute. The precision and
* scale of the attribute have to be applied on the value.
- * ----------
*/
Datum
numeric(PG_FUNCTION_ARGS)
int32 tmp_typmod;
int precision;
int scale;
- int maxweight;
+ int ddigits;
+ int maxdigits;
NumericVar var;
/*
tmp_typmod = typmod - VARHDRSZ;
precision = (tmp_typmod >> 16) & 0xffff;
scale = tmp_typmod & 0xffff;
- maxweight = precision - scale;
+ maxdigits = precision - scale;
/*
- * If the number is in bounds and due to the present result scale no
+ * If the number is certainly in bounds and due to the target scale no
* rounding could be necessary, just make a copy of the input and
- * modify its scale fields.
+ * modify its scale fields. (Note we assume the existing dscale is
+ * honest...)
*/
- if (num->n_weight < maxweight && scale >= num->n_rscale)
+ ddigits = (num->n_weight + 1) * DEC_DIGITS;
+ if (ddigits <= maxdigits && scale >= NUMERIC_DSCALE(num))
{
new = (Numeric) palloc(num->varlen);
memcpy(new, num, num->varlen);
- new->n_rscale = scale;
new->n_sign_dscale = NUMERIC_SIGN(new) |
((uint16) scale & NUMERIC_DSCALE_MASK);
PG_RETURN_NUMERIC(new);
PG_RETURN_NUMERIC(res);
}
-/* ----------
+/*
* numeric_sign() -
*
* returns -1 if the argument is less than 0, 0 if the argument is equal
* to 0, and 1 if the argument is greater than zero.
- * ----------
*/
Datum
numeric_sign(PG_FUNCTION_ARGS)
}
-/* ----------
+/*
* numeric_round() -
*
* Round a value to have 'scale' digits after the decimal point.
* We allow negative 'scale', implying rounding before the decimal
* point --- Oracle interprets rounding that way.
- * ----------
*/
Datum
numeric_round(PG_FUNCTION_ARGS)
int32 scale = PG_GETARG_INT32(1);
Numeric res;
NumericVar arg;
- int i;
/*
* Handle NaN
/*
* Limit the scale value to avoid possible overflow in calculations
- * below.
*/
scale = Max(scale, -NUMERIC_MAX_RESULT_SCALE);
scale = Min(scale, NUMERIC_MAX_RESULT_SCALE);
init_var(&arg);
set_var_from_num(num, &arg);
- i = arg.weight + scale + 1;
-
- if (i < arg.ndigits)
- {
- /*
- * If i = 0, the value loses all digits, but could round up if its
- * first digit is more than 4. If i < 0 the result must be 0.
- */
- if (i < 0)
- arg.ndigits = 0;
- else
- {
- int carry = (arg.digits[i] > 4) ? 1 : 0;
-
- arg.ndigits = i;
-
- while (carry)
- {
- carry += arg.digits[--i];
- arg.digits[i] = carry % 10;
- carry /= 10;
- }
-
- if (i < 0)
- {
- Assert(i == -1); /* better not have added more than 1 digit */
- Assert(arg.digits > arg.buf);
- arg.digits--;
- arg.ndigits++;
- arg.weight++;
- }
- }
- }
+ round_var(&arg, scale);
- /*
- * Set result's scale to something reasonable.
- */
- scale = Min(NUMERIC_MAX_DISPLAY_SCALE, Max(0, scale));
- arg.rscale = scale;
- arg.dscale = scale;
+ /* We don't allow negative output dscale */
+ if (scale < 0)
+ arg.dscale = 0;
/*
* Return the rounded result
}
-/* ----------
+/*
* numeric_trunc() -
*
* Truncate a value to have 'scale' digits after the decimal point.
* We allow negative 'scale', implying a truncation before the decimal
* point --- Oracle interprets truncation that way.
- * ----------
*/
Datum
numeric_trunc(PG_FUNCTION_ARGS)
/*
* Limit the scale value to avoid possible overflow in calculations
- * below.
*/
scale = Max(scale, -NUMERIC_MAX_RESULT_SCALE);
scale = Min(scale, NUMERIC_MAX_RESULT_SCALE);
init_var(&arg);
set_var_from_num(num, &arg);
- arg.ndigits = Min(arg.ndigits, Max(0, arg.weight + scale + 1));
+ trunc_var(&arg, scale);
- /*
- * Set result's scale to something reasonable.
- */
- scale = Min(NUMERIC_MAX_DISPLAY_SCALE, Max(0, scale));
- arg.rscale = scale;
- arg.dscale = scale;
+ /* We don't allow negative output dscale */
+ if (scale < 0)
+ arg.dscale = 0;
/*
* Return the truncated result
}
-/* ----------
+/*
* numeric_ceil() -
*
* Return the smallest integer greater than or equal to the argument
- * ----------
*/
Datum
numeric_ceil(PG_FUNCTION_ARGS)
set_var_from_num(num, &result);
ceil_var(&result, &result);
- result.dscale = 0;
-
res = make_result(&result);
free_var(&result);
}
-/* ----------
+/*
* numeric_floor() -
*
* Return the largest integer equal to or less than the argument
- * ----------
*/
Datum
numeric_floor(PG_FUNCTION_ARGS)
set_var_from_num(num, &result);
floor_var(&result, &result);
- result.dscale = 0;
-
res = make_result(&result);
free_var(&result);
/* ----------------------------------------------------------------------
*
- * Arithmetic base functions
+ * Basic arithmetic functions
*
* ----------------------------------------------------------------------
*/
-/* ----------
+/*
* numeric_add() -
*
* Add two numerics
- * ----------
*/
Datum
numeric_add(PG_FUNCTION_ARGS)
/*
* Unpack the values, let add_var() compute the result and return it.
- * The internals of add_var() will automatically set the correct
- * result and display scales in the result.
*/
init_var(&arg1);
init_var(&arg2);
set_var_from_num(num2, &arg2);
add_var(&arg1, &arg2, &result);
+
res = make_result(&result);
free_var(&arg1);
}
-/* ----------
+/*
* numeric_sub() -
*
* Subtract one numeric from another
- * ----------
*/
Datum
numeric_sub(PG_FUNCTION_ARGS)
PG_RETURN_NUMERIC(make_result(&const_nan));
/*
- * Unpack the two arguments, let sub_var() compute the result and
- * return it.
+ * Unpack the values, let sub_var() compute the result and return it.
*/
init_var(&arg1);
init_var(&arg2);
set_var_from_num(num2, &arg2);
sub_var(&arg1, &arg2, &result);
+
res = make_result(&result);
free_var(&arg1);
}
-/* ----------
+/*
* numeric_mul() -
*
* Calculate the product of two numerics
- * ----------
*/
Datum
numeric_mul(PG_FUNCTION_ARGS)
PG_RETURN_NUMERIC(make_result(&const_nan));
/*
- * Unpack the arguments, let mul_var() compute the result and return
- * it. Unlike add_var() and sub_var(), mul_var() will round the result
- * to the scale stored in global_rscale. In the case of numeric_mul(),
- * which is invoked for the * operator on numerics, we set it to the
- * exact representation for the product (rscale = sum(rscale of arg1,
- * rscale of arg2) and the same for the dscale).
+ * Unpack the values, let mul_var() compute the result and return it.
+ * Unlike add_var() and sub_var(), mul_var() will round its result.
+ * In the case of numeric_mul(), which is invoked for the * operator on
+ * numerics, we request exact representation for the product (rscale =
+ * sum(dscale of arg1, dscale of arg2)).
*/
init_var(&arg1);
init_var(&arg2);
set_var_from_num(num1, &arg1);
set_var_from_num(num2, &arg2);
- global_rscale = arg1.rscale + arg2.rscale;
-
- mul_var(&arg1, &arg2, &result);
-
- result.dscale = arg1.dscale + arg2.dscale;
+ mul_var(&arg1, &arg2, &result, arg1.dscale + arg2.dscale);
res = make_result(&result);
}
-/* ----------
+/*
* numeric_div() -
*
* Divide one numeric into another
- * ----------
*/
Datum
numeric_div(PG_FUNCTION_ARGS)
NumericVar arg2;
NumericVar result;
Numeric res;
- int res_dscale;
+ int rscale;
/*
* Handle NaN
set_var_from_num(num1, &arg1);
set_var_from_num(num2, &arg2);
- res_dscale = select_div_scale(&arg1, &arg2);
-
/*
- * Do the divide, set the display scale and return the result
+ * Select scale for division result
*/
- div_var(&arg1, &arg2, &result);
+ rscale = select_div_scale(&arg1, &arg2);
- result.dscale = res_dscale;
+ /*
+ * Do the divide and return the result
+ */
+ div_var(&arg1, &arg2, &result, rscale);
res = make_result(&result);
}
-/* ----------
+/*
* numeric_mod() -
*
* Calculate the modulo of two numerics
- * ----------
*/
Datum
numeric_mod(PG_FUNCTION_ARGS)
}
-/* ----------
+/*
* numeric_inc() -
*
* Increment a number by one
- * ----------
*/
Datum
numeric_inc(PG_FUNCTION_ARGS)
set_var_from_num(num, &arg);
add_var(&arg, &const_one, &arg);
+
res = make_result(&arg);
free_var(&arg);
}
-/* ----------
+/*
* numeric_smaller() -
*
* Return the smaller of two numbers
- * ----------
*/
Datum
numeric_smaller(PG_FUNCTION_ARGS)
}
-/* ----------
+/*
* numeric_larger() -
*
* Return the larger of two numbers
- * ----------
*/
Datum
numeric_larger(PG_FUNCTION_ARGS)
/* ----------------------------------------------------------------------
*
- * Complex math functions
+ * Advanced math functions
*
* ----------------------------------------------------------------------
*/
-/* ----------
+/*
* numeric_sqrt() -
*
* Compute the square root of a numeric.
- * ----------
*/
Datum
numeric_sqrt(PG_FUNCTION_ARGS)
NumericVar arg;
NumericVar result;
int sweight;
- int res_dscale;
+ int rscale;
/*
* Handle NaN
PG_RETURN_NUMERIC(make_result(&const_nan));
/*
- * Unpack the argument and determine the scales. We choose a display
+ * Unpack the argument and determine the result scale. We choose a
* scale to give at least NUMERIC_MIN_SIG_DIGITS significant digits;
* but in any case not less than the input's dscale.
*/
set_var_from_num(num, &arg);
/* Assume the input was normalized, so arg.weight is accurate */
- sweight = (arg.weight / 2) - 1;
-
- res_dscale = NUMERIC_MIN_SIG_DIGITS - sweight;
- res_dscale = Max(res_dscale, arg.dscale);
- res_dscale = Max(res_dscale, NUMERIC_MIN_DISPLAY_SCALE);
- res_dscale = Min(res_dscale, NUMERIC_MAX_DISPLAY_SCALE);
+ sweight = (arg.weight + 1) * DEC_DIGITS / 2 - 1;
- global_rscale = res_dscale + 8;
+ rscale = NUMERIC_MIN_SIG_DIGITS - sweight;
+ rscale = Max(rscale, arg.dscale);
+ rscale = Max(rscale, NUMERIC_MIN_DISPLAY_SCALE);
+ rscale = Min(rscale, NUMERIC_MAX_DISPLAY_SCALE);
/*
* Let sqrt_var() do the calculation and return the result.
*/
- sqrt_var(&arg, &result);
-
- result.dscale = res_dscale;
+ sqrt_var(&arg, &result, rscale);
res = make_result(&result);
}
-/* ----------
+/*
* numeric_exp() -
*
* Raise e to the power of x
- * ----------
*/
Datum
numeric_exp(PG_FUNCTION_ARGS)
Numeric res;
NumericVar arg;
NumericVar result;
- int res_dscale;
+ int rscale;
double val;
/*
PG_RETURN_NUMERIC(make_result(&const_nan));
/*
- * Unpack the argument and determine the scales. We choose a display
+ * Unpack the argument and determine the result scale. We choose a
* scale to give at least NUMERIC_MIN_SIG_DIGITS significant digits;
* but in any case not less than the input's dscale.
*/
set_var_from_num(num, &arg);
/* convert input to float8, ignoring overflow */
- val = numeric_to_double_no_overflow(num);
+ val = numericvar_to_double_no_overflow(&arg);
- /* log10(result) = num * log10(e), so this is approximately the weight: */
+ /*
+ * log10(result) = num * log10(e), so this is approximately the decimal
+ * weight of the result:
+ */
val *= 0.434294481903252;
/* limit to something that won't cause integer overflow */
val = Max(val, -NUMERIC_MAX_RESULT_SCALE);
val = Min(val, NUMERIC_MAX_RESULT_SCALE);
- res_dscale = NUMERIC_MIN_SIG_DIGITS - (int) val;
- res_dscale = Max(res_dscale, arg.dscale);
- res_dscale = Max(res_dscale, NUMERIC_MIN_DISPLAY_SCALE);
- res_dscale = Min(res_dscale, NUMERIC_MAX_DISPLAY_SCALE);
-
- global_rscale = res_dscale + NUMERIC_EXTRA_DIGITS;
+ rscale = NUMERIC_MIN_SIG_DIGITS - (int) val;
+ rscale = Max(rscale, arg.dscale);
+ rscale = Max(rscale, NUMERIC_MIN_DISPLAY_SCALE);
+ rscale = Min(rscale, NUMERIC_MAX_DISPLAY_SCALE);
/*
* Let exp_var() do the calculation and return the result.
*/
- exp_var(&arg, &result);
-
- result.dscale = res_dscale;
+ exp_var(&arg, &result, rscale);
res = make_result(&result);
}
-/* ----------
+/*
* numeric_ln() -
*
* Compute the natural logarithm of x
- * ----------
*/
Datum
numeric_ln(PG_FUNCTION_ARGS)
Numeric res;
NumericVar arg;
NumericVar result;
- int res_dscale;
+ int dec_digits;
+ int rscale;
/*
* Handle NaN
set_var_from_num(num, &arg);
- if (arg.weight > 0)
- res_dscale = NUMERIC_MIN_SIG_DIGITS - (int) log10(arg.weight);
- else if (arg.weight < 0)
- res_dscale = NUMERIC_MIN_SIG_DIGITS - (int) log10(- arg.weight);
- else
- res_dscale = NUMERIC_MIN_SIG_DIGITS;
-
- res_dscale = Max(res_dscale, arg.dscale);
- res_dscale = Max(res_dscale, NUMERIC_MIN_DISPLAY_SCALE);
- res_dscale = Min(res_dscale, NUMERIC_MAX_DISPLAY_SCALE);
+ /* Approx decimal digits before decimal point */
+ dec_digits = (arg.weight + 1) * DEC_DIGITS;
- global_rscale = res_dscale + NUMERIC_EXTRA_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;
- ln_var(&arg, &result);
+ rscale = Max(rscale, arg.dscale);
+ rscale = Max(rscale, NUMERIC_MIN_DISPLAY_SCALE);
+ rscale = Min(rscale, NUMERIC_MAX_DISPLAY_SCALE);
- result.dscale = res_dscale;
+ ln_var(&arg, &result, rscale);
res = make_result(&result);
}
-/* ----------
+/*
* numeric_log() -
*
* Compute the logarithm of x in a given base
- * ----------
*/
Datum
numeric_log(PG_FUNCTION_ARGS)
/*
* Call log_var() to compute and return the result; note it handles
- * rscale/dscale itself.
+ * scale selection itself.
*/
log_var(&arg1, &arg2, &result);
}
-/* ----------
+/*
* numeric_power() -
*
* Raise b to the power of x
- * ----------
*/
Datum
numeric_power(PG_FUNCTION_ARGS)
/*
* Call power_var() to compute and return the result; note it handles
- * rscale/dscale itself.
+ * scale selection itself.
*/
power_var(&arg1, &arg2, &result);
int32 val = PG_GETARG_INT32(0);
Numeric res;
NumericVar result;
- char *tmp;
init_var(&result);
- tmp = DatumGetCString(DirectFunctionCall1(int4out,
- Int32GetDatum(val)));
- set_var_from_str(tmp, &result);
+ int8_to_numericvar((int64) val, &result);
+
res = make_result(&result);
free_var(&result);
- pfree(tmp);
PG_RETURN_NUMERIC(res);
}
{
Numeric num = PG_GETARG_NUMERIC(0);
NumericVar x;
- char *str;
- Datum result;
+ int64 val;
+ int32 result;
/* XXX would it be better to return NULL? */
if (NUMERIC_IS_NAN(num))
elog(ERROR, "Cannot convert NaN to int4");
- /*
- * Get the number in the variable format so we can round to integer.
- */
+ /* Convert to variable format and thence to int8 */
init_var(&x);
set_var_from_num(num, &x);
- str = get_str_from_var(&x, 0); /* dscale = 0 produces rounding */
+ if (!numericvar_to_int8(&x, &val))
+ elog(ERROR, "numeric conversion to int4 is out of range");
free_var(&x);
- result = DirectFunctionCall1(int4in, CStringGetDatum(str));
- pfree(str);
+ /* Down-convert to int4 */
+ result = (int32) val;
- PG_RETURN_DATUM(result);
+ /* Test for overflow by reverse-conversion. */
+ if ((int64) result != val)
+ elog(ERROR, "numeric conversion to int4 is out of range");
+
+ PG_RETURN_INT32(result);
}
Datum
int8_numeric(PG_FUNCTION_ARGS)
{
- Datum val = PG_GETARG_DATUM(0);
+ int64 val = PG_GETARG_INT64(0);
Numeric res;
NumericVar result;
- char *tmp;
init_var(&result);
- tmp = DatumGetCString(DirectFunctionCall1(int8out, val));
- set_var_from_str(tmp, &result);
+ int8_to_numericvar(val, &result);
+
res = make_result(&result);
free_var(&result);
- pfree(tmp);
PG_RETURN_NUMERIC(res);
}
{
Numeric num = PG_GETARG_NUMERIC(0);
NumericVar x;
- char *str;
- Datum result;
+ int64 result;
/* XXX would it be better to return NULL? */
if (NUMERIC_IS_NAN(num))
elog(ERROR, "Cannot convert NaN to int8");
- /*
- * Get the number in the variable format so we can round to integer.
- */
+ /* Convert to variable format and thence to int8 */
init_var(&x);
set_var_from_num(num, &x);
- str = get_str_from_var(&x, 0); /* dscale = 0 produces rounding */
+ if (!numericvar_to_int8(&x, &result))
+ elog(ERROR, "numeric conversion to int8 is out of range");
free_var(&x);
- result = DirectFunctionCall1(int8in, CStringGetDatum(str));
- pfree(str);
-
- PG_RETURN_DATUM(result);
+ PG_RETURN_INT64(result);
}
int16 val = PG_GETARG_INT16(0);
Numeric res;
NumericVar result;
- char *tmp;
init_var(&result);
- tmp = DatumGetCString(DirectFunctionCall1(int2out,
- Int16GetDatum(val)));
- set_var_from_str(tmp, &result);
+ int8_to_numericvar((int64) val, &result);
+
res = make_result(&result);
free_var(&result);
- pfree(tmp);
PG_RETURN_NUMERIC(res);
}
{
Numeric num = PG_GETARG_NUMERIC(0);
NumericVar x;
- char *str;
- Datum result;
+ int64 val;
+ int16 result;
/* XXX would it be better to return NULL? */
if (NUMERIC_IS_NAN(num))
elog(ERROR, "Cannot convert NaN to int2");
- /*
- * Get the number in the variable format so we can round to integer.
- */
+ /* Convert to variable format and thence to int8 */
init_var(&x);
set_var_from_num(num, &x);
- str = get_str_from_var(&x, 0); /* dscale = 0 produces rounding */
+ if (!numericvar_to_int8(&x, &val))
+ elog(ERROR, "numeric conversion to int2 is out of range");
free_var(&x);
- result = DirectFunctionCall1(int2in, CStringGetDatum(str));
- pfree(str);
+ /* Down-convert to int2 */
+ result = (int16) val;
- return result;
+ /* Test for overflow by reverse-conversion. */
+ if ((int64) result != val)
+ elog(ERROR, "numeric conversion to int2 is out of range");
+
+ PG_RETURN_INT16(result);
}
vsumX,
vsumX2,
vNminus1;
- int div_dscale;
+ int rscale;
/* We assume the input is array of numeric */
deconstruct_array(transarray,
init_var(&vsumX2);
set_var_from_num(sumX2, &vsumX2);
- /* set rscale for mul_var calls */
- global_rscale = vsumX.rscale * 2;
+ /* compute rscale for mul_var calls */
+ rscale = vsumX.dscale * 2;
- mul_var(&vsumX, &vsumX, &vsumX); /* now vsumX contains sumX * sumX */
- mul_var(&vN, &vsumX2, &vsumX2); /* now vsumX2 contains N * sumX2 */
+ mul_var(&vsumX, &vsumX, &vsumX, rscale); /* vsumX = sumX * sumX */
+ mul_var(&vN, &vsumX2, &vsumX2, rscale); /* vsumX2 = N * sumX2 */
sub_var(&vsumX2, &vsumX, &vsumX2); /* N * sumX2 - sumX * sumX */
if (cmp_var(&vsumX2, &const_zero) <= 0)
}
else
{
- mul_var(&vN, &vNminus1, &vNminus1); /* N * (N - 1) */
- div_dscale = select_div_scale(&vsumX2, &vNminus1);
- div_var(&vsumX2, &vNminus1, &vsumX); /* variance */
- vsumX.dscale = div_dscale;
+ mul_var(&vN, &vNminus1, &vNminus1, 0); /* N * (N - 1) */
+ rscale = select_div_scale(&vsumX2, &vNminus1);
+ div_var(&vsumX2, &vNminus1, &vsumX, rscale); /* variance */
res = make_result(&vsumX);
}
vsumX,
vsumX2,
vNminus1;
- int div_dscale;
+ int rscale;
/* We assume the input is array of numeric */
deconstruct_array(transarray,
init_var(&vsumX2);
set_var_from_num(sumX2, &vsumX2);
- /* set rscale for mul_var calls */
- global_rscale = vsumX.rscale * 2;
+ /* compute rscale for mul_var calls */
+ rscale = vsumX.dscale * 2;
- mul_var(&vsumX, &vsumX, &vsumX); /* now vsumX contains sumX * sumX */
- mul_var(&vN, &vsumX2, &vsumX2); /* now vsumX2 contains N * sumX2 */
+ mul_var(&vsumX, &vsumX, &vsumX, rscale); /* vsumX = sumX * sumX */
+ mul_var(&vN, &vsumX2, &vsumX2, rscale); /* vsumX2 = N * sumX2 */
sub_var(&vsumX2, &vsumX, &vsumX2); /* N * sumX2 - sumX * sumX */
if (cmp_var(&vsumX2, &const_zero) <= 0)
}
else
{
- mul_var(&vN, &vNminus1, &vNminus1); /* N * (N - 1) */
- div_dscale = select_div_scale(&vsumX2, &vNminus1);
- div_var(&vsumX2, &vNminus1, &vsumX); /* variance */
- vsumX.dscale = div_dscale;
- sqrt_var(&vsumX, &vsumX); /* stddev */
+ mul_var(&vN, &vNminus1, &vNminus1, 0); /* N * (N - 1) */
+ rscale = select_div_scale(&vsumX2, &vNminus1);
+ div_var(&vsumX2, &vNminus1, &vsumX, rscale); /* variance */
+ sqrt_var(&vsumX, &vsumX, rscale); /* stddev */
res = make_result(&vsumX);
}
/* ----------------------------------------------------------------------
*
- * Local functions follow
+ * Debug support
*
* ----------------------------------------------------------------------
*/
-
#ifdef NUMERIC_DEBUG
-/* ----------
+/*
* dump_numeric() - Dump a value in the db storage format for debugging
- * ----------
*/
static void
-dump_numeric(char *str, Numeric num)
+dump_numeric(const char *str, Numeric num)
{
+ NumericDigit *digits = (NumericDigit *) num->n_data;
+ int ndigits;
int i;
- printf("%s: NUMERIC w=%d r=%d d=%d ", str, num->n_weight, num->n_rscale,
- NUMERIC_DSCALE(num));
+ ndigits = (num->varlen - NUMERIC_HDRSZ) / sizeof(NumericDigit);
+
+ printf("%s: NUMERIC w=%d d=%d ", str, num->n_weight, NUMERIC_DSCALE(num));
switch (NUMERIC_SIGN(num))
{
case NUMERIC_POS:
break;
}
- for (i = 0; i < num->varlen - NUMERIC_HDRSZ; i++)
- printf(" %d %d", (num->n_data[i] >> 4) & 0x0f, num->n_data[i] & 0x0f);
+ for (i = 0; i < ndigits; i++)
+ printf(" %0*d", DEC_DIGITS, digits[i]);
printf("\n");
}
-/* ----------
+/*
* dump_var() - Dump a value in the variable format for debugging
- * ----------
*/
static void
-dump_var(char *str, NumericVar *var)
+dump_var(const char *str, NumericVar *var)
{
int i;
- printf("%s: VAR w=%d r=%d d=%d ", str, var->weight, var->rscale,
- var->dscale);
+ printf("%s: VAR w=%d d=%d ", str, var->weight, var->dscale);
switch (var->sign)
{
case NUMERIC_POS:
}
for (i = 0; i < var->ndigits; i++)
- printf(" %d", var->digits[i]);
+ printf(" %0*d", DEC_DIGITS, var->digits[i]);
printf("\n");
}
+
#endif /* NUMERIC_DEBUG */
-/* ----------
+/* ----------------------------------------------------------------------
+ *
+ * Local functions follow
+ *
+ * In general, these do not support NaNs --- callers must eliminate
+ * the possibility of NaN first. (make_result() is an exception.)
+ *
+ * ----------------------------------------------------------------------
+ */
+
+
+/*
* alloc_var() -
*
* Allocate a digit buffer of ndigits digits (plus a spare digit for rounding)
- * ----------
*/
static void
alloc_var(NumericVar *var, int ndigits)
{
digitbuf_free(var->buf);
var->buf = digitbuf_alloc(ndigits + 1);
- var->buf[0] = 0;
+ var->buf[0] = 0; /* spare digit for rounding */
var->digits = var->buf + 1;
var->ndigits = ndigits;
}
-/* ----------
+/*
* free_var() -
*
* Return the digit buffer of a variable to the free pool
- * ----------
*/
static void
free_var(NumericVar *var)
}
-/* ----------
+/*
* zero_var() -
*
* Set a variable to ZERO.
- * Note: rscale and dscale are not touched.
- * ----------
+ * Note: its dscale is not touched.
*/
static void
zero_var(NumericVar *var)
}
-/* ----------
+/*
* set_var_from_str()
*
* Parse a string and put the number into a variable
- * ----------
*/
static void
-set_var_from_str(char *str, NumericVar *dest)
+set_var_from_str(const char *str, NumericVar *dest)
{
- char *cp = str;
+ const char *cp = str;
bool have_dp = FALSE;
- int i = 0;
+ int i;
+ unsigned char *decdigits;
+ int sign = NUMERIC_POS;
+ int dweight = -1;
+ int ddigits;
+ int dscale = 0;
+ int weight;
+ int ndigits;
+ int offset;
+ NumericDigit *digits;
+
+ /*
+ * We first parse the string to extract decimal digits and determine the
+ * correct decimal weight. Then convert to NBASE representation.
+ */
+ /* skip leading spaces */
while (*cp)
{
if (!isspace((unsigned char) *cp))
cp++;
}
- alloc_var(dest, strlen(cp));
- dest->weight = -1;
- dest->dscale = 0;
- dest->sign = NUMERIC_POS;
-
switch (*cp)
{
case '+':
- dest->sign = NUMERIC_POS;
+ sign = NUMERIC_POS;
cp++;
break;
case '-':
- dest->sign = NUMERIC_NEG;
+ sign = NUMERIC_NEG;
cp++;
break;
}
if (!isdigit((unsigned char) *cp))
elog(ERROR, "Bad numeric input format '%s'", str);
+ decdigits = (unsigned char *) palloc(strlen(cp) + DEC_DIGITS*2);
+
+ /* leading padding for digit alignment later */
+ memset(decdigits, 0, DEC_DIGITS);
+ i = DEC_DIGITS;
+
while (*cp)
{
if (isdigit((unsigned char) *cp))
{
- dest->digits[i++] = *cp++ - '0';
+ decdigits[i++] = *cp++ - '0';
if (!have_dp)
- dest->weight++;
+ dweight++;
else
- dest->dscale++;
+ dscale++;
}
else if (*cp == '.')
{
else
break;
}
- dest->ndigits = i;
+
+ ddigits = i - DEC_DIGITS;
+ /* trailing padding for digit alignment later */
+ memset(decdigits + i, 0, DEC_DIGITS-1);
/* Handle exponent, if any */
if (*cp == 'e' || *cp == 'E')
if (exponent > NUMERIC_MAX_PRECISION ||
exponent < -NUMERIC_MAX_PRECISION)
elog(ERROR, "Bad numeric input format '%s'", str);
- dest->weight += (int) exponent;
- dest->dscale -= (int) exponent;
- if (dest->dscale < 0)
- dest->dscale = 0;
+ dweight += (int) exponent;
+ dscale -= (int) exponent;
+ if (dscale < 0)
+ dscale = 0;
}
/* Should be nothing left but spaces */
cp++;
}
- /* Strip any leading zeroes */
- while (dest->ndigits > 0 && *(dest->digits) == 0)
+ /*
+ * Okay, convert pure-decimal representation to base NBASE. First we
+ * need to determine the converted weight and ndigits. offset is the
+ * number of decimal zeroes to insert before the first given digit to
+ * have a correctly aligned first NBASE digit.
+ */
+ if (dweight >= 0)
+ weight = (dweight + 1 + DEC_DIGITS-1) / DEC_DIGITS - 1;
+ else
+ weight = - ((-dweight - 1) / DEC_DIGITS + 1);
+ offset = (weight + 1) * DEC_DIGITS - (dweight + 1);
+ ndigits = (ddigits + offset + DEC_DIGITS-1) / DEC_DIGITS;
+
+ alloc_var(dest, ndigits);
+ dest->sign = sign;
+ dest->weight = weight;
+ dest->dscale = dscale;
+
+ i = DEC_DIGITS - offset;
+ digits = dest->digits;
+
+ while (ndigits-- > 0)
{
- (dest->digits)++;
- (dest->weight)--;
- (dest->ndigits)--;
+#if DEC_DIGITS == 4
+ *digits++ = ((decdigits[i] * 10 + decdigits[i+1]) * 10 +
+ decdigits[i+2]) * 10 + decdigits[i+3];
+#elif DEC_DIGITS == 2
+ *digits++ = decdigits[i] * 10 + decdigits[i+1];
+#elif DEC_DIGITS == 1
+ *digits++ = decdigits[i];
+#else
+#error unsupported NBASE
+#endif
+ i += DEC_DIGITS;
}
- if (dest->ndigits == 0)
- dest->weight = 0;
- dest->rscale = dest->dscale;
+ pfree(decdigits);
+
+ /* Strip any leading/trailing zeroes, and normalize weight if zero */
+ strip_var(dest);
}
/*
* set_var_from_num() -
*
- * Parse back the packed db format into a variable
- *
+ * Convert the packed db format into a variable
*/
static void
set_var_from_num(Numeric num, NumericVar *dest)
{
- NumericDigit *digit;
- int i;
- int n;
+ int ndigits;
- n = num->varlen - NUMERIC_HDRSZ; /* number of digit-pairs in packed
- * fmt */
+ ndigits = (num->varlen - NUMERIC_HDRSZ) / sizeof(NumericDigit);
- alloc_var(dest, n * 2);
+ alloc_var(dest, ndigits);
dest->weight = num->n_weight;
- dest->rscale = num->n_rscale;
- dest->dscale = NUMERIC_DSCALE(num);
dest->sign = NUMERIC_SIGN(num);
+ dest->dscale = NUMERIC_DSCALE(num);
- digit = dest->digits;
-
- for (i = 0; i < n; i++)
- {
- unsigned char digitpair = num->n_data[i];
-
- *digit++ = (digitpair >> 4) & 0x0f;
- *digit++ = digitpair & 0x0f;
- }
+ memcpy(dest->digits, num->n_data, ndigits * sizeof(NumericDigit));
}
-/* ----------
+/*
* set_var_from_var() -
*
* Copy one variable into another
- * ----------
*/
static void
set_var_from_var(NumericVar *value, NumericVar *dest)
newbuf = digitbuf_alloc(value->ndigits + 1);
newbuf[0] = 0; /* spare digit for rounding */
- memcpy(newbuf + 1, value->digits, value->ndigits);
+ memcpy(newbuf + 1, value->digits, value->ndigits * sizeof(NumericDigit));
digitbuf_free(dest->buf);
}
-/* ----------
+/*
* get_str_from_var() -
*
* Convert a var to text representation (guts of numeric_out).
* CAUTION: var's contents may be modified by rounding!
- * Caller must have checked for NaN case.
* Returns a palloc'd string.
- * ----------
*/
static char *
get_str_from_var(NumericVar *var, int dscale)
{
char *str;
char *cp;
+ char *endcp;
int i;
int d;
+ NumericDigit dig;
+#if DEC_DIGITS > 1
+ NumericDigit d1;
+#endif
+
+ if (dscale < 0)
+ dscale = 0;
/*
* Check if we must round up before printing the value and do so.
*/
- i = dscale + var->weight + 1;
- if (i >= 0 && var->ndigits > i)
- {
- int carry = (var->digits[i] > 4) ? 1 : 0;
-
- var->ndigits = i;
-
- while (carry)
- {
- carry += var->digits[--i];
- var->digits[i] = carry % 10;
- carry /= 10;
- }
-
- if (i < 0)
- {
- Assert(i == -1); /* better not have added more than 1 digit */
- Assert(var->digits > var->buf);
- var->digits--;
- var->ndigits++;
- var->weight++;
- }
- }
- else
- var->ndigits = Max(0, Min(i, var->ndigits));
+ round_var(var, dscale);
/*
- * Allocate space for the result
+ * Allocate space for the result.
+ *
+ * i is set to to # of decimal digits before decimal point.
+ * dscale is the # of decimal digits we will print after decimal point.
+ * We may generate as many as DEC_DIGITS-1 excess digits at the end,
+ * and in addition we need room for sign, decimal point, null terminator.
*/
- str = palloc(Max(0, dscale) + Max(0, var->weight) + 4);
+ i = (var->weight + 1) * DEC_DIGITS;
+ if (i <= 0)
+ i = 1;
+
+ str = palloc(i + dscale + DEC_DIGITS + 2);
cp = str;
/*
/*
* Output all digits before the decimal point
*/
- i = Max(var->weight, 0);
- d = 0;
-
- while (i >= 0)
+ if (var->weight < 0)
{
- if (i <= var->weight && d < var->ndigits)
- *cp++ = var->digits[d++] + '0';
- else
- *cp++ = '0';
- i--;
+ d = var->weight + 1;
+ *cp++ = '0';
+ }
+ else
+ {
+ for (d = 0; d <= var->weight; d++)
+ {
+ dig = (d < var->ndigits) ? var->digits[d] : 0;
+ /* In the first digit, suppress extra leading decimal zeroes */
+#if DEC_DIGITS == 4
+ {
+ bool putit = (d > 0);
+
+ d1 = dig / 1000;
+ dig -= d1 * 1000;
+ putit |= (d1 > 0);
+ if (putit)
+ *cp++ = d1 + '0';
+ d1 = dig / 100;
+ dig -= d1 * 100;
+ putit |= (d1 > 0);
+ if (putit)
+ *cp++ = d1 + '0';
+ d1 = dig / 10;
+ dig -= d1 * 10;
+ putit |= (d1 > 0);
+ if (putit)
+ *cp++ = d1 + '0';
+ *cp++ = dig + '0';
+ }
+#elif DEC_DIGITS == 2
+ d1 = dig / 10;
+ dig -= d1 * 10;
+ if (d1 > 0 || d > 0)
+ *cp++ = d1 + '0';
+ *cp++ = dig + '0';
+#elif DEC_DIGITS == 1
+ *cp++ = dig + '0';
+#else
+#error unsupported NBASE
+#endif
+ }
}
/*
* If requested, output a decimal point and all the digits that follow
- * it.
+ * it. We initially put out a multiple of DEC_DIGITS digits, then
+ * truncate if needed.
*/
if (dscale > 0)
{
*cp++ = '.';
- while (i >= -dscale)
+ endcp = cp + dscale;
+ for (i = 0; i < dscale; d++, i += DEC_DIGITS)
{
- if (i <= var->weight && d < var->ndigits)
- *cp++ = var->digits[d++] + '0';
- else
- *cp++ = '0';
- i--;
+ dig = (d >= 0 && d < var->ndigits) ? var->digits[d] : 0;
+#if DEC_DIGITS == 4
+ d1 = dig / 1000;
+ dig -= d1 * 1000;
+ *cp++ = d1 + '0';
+ d1 = dig / 100;
+ dig -= d1 * 100;
+ *cp++ = d1 + '0';
+ d1 = dig / 10;
+ dig -= d1 * 10;
+ *cp++ = d1 + '0';
+ *cp++ = dig + '0';
+#elif DEC_DIGITS == 2
+ d1 = dig / 10;
+ dig -= d1 * 10;
+ *cp++ = d1 + '0';
+ *cp++ = dig + '0';
+#elif DEC_DIGITS == 1
+ *cp++ = dig + '0';
+#else
+#error unsupported NBASE
+#endif
}
+ cp = endcp;
}
/*
}
-/* ----------
+/*
* make_result() -
*
* Create the packed db numeric format in palloc()'d memory from
- * a variable. The var's rscale determines the number of digits kept.
- * ----------
+ * a variable.
*/
static Numeric
make_result(NumericVar *var)
{
Numeric result;
- NumericDigit *digit = var->digits;
+ NumericDigit *digits = var->digits;
int weight = var->weight;
int sign = var->sign;
int n;
- int i,
- j;
+ Size len;
if (sign == NUMERIC_NAN)
{
result->varlen = NUMERIC_HDRSZ;
result->n_weight = 0;
- result->n_rscale = 0;
result->n_sign_dscale = NUMERIC_NAN;
dump_numeric("make_result()", result);
return result;
}
- n = Max(0, Min(var->ndigits, var->weight + var->rscale + 1));
+ n = var->ndigits;
/* truncate leading zeroes */
- while (n > 0 && *digit == 0)
+ while (n > 0 && *digits == 0)
{
- digit++;
+ digits++;
weight--;
n--;
}
/* truncate trailing zeroes */
- while (n > 0 && digit[n - 1] == 0)
+ while (n > 0 && digits[n - 1] == 0)
n--;
/* If zero result, force to weight=0 and positive sign */
sign = NUMERIC_POS;
}
- result = (Numeric) palloc(NUMERIC_HDRSZ + (n + 1) / 2);
- result->varlen = NUMERIC_HDRSZ + (n + 1) / 2;
+ /* Build the result */
+ len = NUMERIC_HDRSZ + n * sizeof(NumericDigit);
+ result = (Numeric) palloc(len);
+ result->varlen = len;
result->n_weight = weight;
- result->n_rscale = var->rscale;
- result->n_sign_dscale = sign |
- ((uint16) var->dscale & NUMERIC_DSCALE_MASK);
+ result->n_sign_dscale = sign | (var->dscale & NUMERIC_DSCALE_MASK);
- i = 0;
- j = 0;
- while (j < n)
- {
- unsigned char digitpair = digit[j++] << 4;
+ memcpy(result->n_data, digits, n * sizeof(NumericDigit));
- if (j < n)
- digitpair |= digit[j++];
- result->n_data[i++] = digitpair;
- }
+ /* Check for overflow of int16 fields */
+ if (result->n_weight != weight ||
+ NUMERIC_DSCALE(result) != var->dscale)
+ elog(ERROR, "Value overflows numeric format");
dump_numeric("make_result()", result);
return result;
}
-/* ----------
+/*
* apply_typmod() -
*
* Do bounds checking and rounding according to the attributes
* typmod field.
- * ----------
*/
static void
apply_typmod(NumericVar *var, int32 typmod)
{
int precision;
int scale;
- int maxweight;
+ int maxdigits;
+ int ddigits;
int i;
/* Do nothing if we have a default typmod (-1) */
typmod -= VARHDRSZ;
precision = (typmod >> 16) & 0xffff;
scale = typmod & 0xffff;
- maxweight = precision - scale;
-
- /* Round to target scale */
- i = scale + var->weight + 1;
- if (i >= 0 && var->ndigits > i)
- {
- int carry = (var->digits[i] > 4) ? 1 : 0;
-
- var->ndigits = i;
+ maxdigits = precision - scale;
- while (carry)
- {
- carry += var->digits[--i];
- var->digits[i] = carry % 10;
- carry /= 10;
- }
-
- if (i < 0)
- {
- Assert(i == -1); /* better not have added more than 1 digit */
- Assert(var->digits > var->buf);
- var->digits--;
- var->ndigits++;
- var->weight++;
- }
- }
- else
- var->ndigits = Max(0, Min(i, var->ndigits));
+ /* Round to target scale (and set var->dscale) */
+ round_var(var, scale);
/*
* Check for overflow - note we can't do this before rounding, because
* storage but perhaps might not have been yet. In any case, we must
* recognize a true zero, whose weight doesn't mean anything.
*/
- if (var->weight >= maxweight)
+ ddigits = (var->weight + 1) * DEC_DIGITS;
+ if (ddigits > maxdigits)
{
/* Determine true weight; and check for all-zero result */
- int tweight = var->weight;
-
for (i = 0; i < var->ndigits; i++)
{
- if (var->digits[i])
+ NumericDigit dig = var->digits[i];
+
+ if (dig)
+ {
+ /* Adjust for any high-order decimal zero digits */
+#if DEC_DIGITS == 4
+ if (dig < 10)
+ ddigits -= 3;
+ else if (dig < 100)
+ ddigits -= 2;
+ else if (dig < 1000)
+ ddigits -= 1;
+#elif DEC_DIGITS == 2
+ if (dig < 10)
+ ddigits -= 1;
+#elif DEC_DIGITS == 1
+ /* no adjustment */
+#else
+#error unsupported NBASE
+#endif
+ if (ddigits > maxdigits)
+ elog(ERROR, "overflow on numeric "
+ "ABS(value) >= 10^%d for field with precision %d scale %d",
+ ddigits-1, precision, scale);
break;
- tweight--;
+ }
+ ddigits -= DEC_DIGITS;
}
+ }
+}
+
+/*
+ * Convert numeric to int8, rounding if needed.
+ *
+ * If overflow, return FALSE (no error is raised). Return TRUE if okay.
+ *
+ * CAUTION: var's contents may be modified by rounding!
+ */
+static bool
+numericvar_to_int8(NumericVar *var, int64 *result)
+{
+ NumericDigit *digits;
+ int ndigits;
+ int i;
+ int64 val,
+ oldval;
+ bool neg;
+
+ /* Round to nearest integer */
+ round_var(var, 0);
- if (tweight >= maxweight && i < var->ndigits)
- elog(ERROR, "overflow on numeric "
- "ABS(value) >= 10^%d for field with precision %d scale %d",
- tweight, precision, scale);
+ /* Check for zero input */
+ strip_var(var);
+ ndigits = var->ndigits;
+ if (ndigits == 0)
+ {
+ *result = 0;
+ return true;
}
- var->rscale = scale;
- var->dscale = scale;
+ /* Construct the result */
+ digits = var->digits;
+ neg = (var->sign == NUMERIC_NEG);
+ val = digits[0];
+ for (i = 1; i < ndigits; i++)
+ {
+ oldval = val;
+ val *= NBASE;
+ val += digits[i];
+ /*
+ * The overflow check is a bit tricky because we want to accept
+ * INT64_MIN, which will overflow the positive accumulator. We
+ * can detect this case easily though because INT64_MIN is the
+ * only nonzero value for which -val == val (on a two's complement
+ * machine, anyway).
+ */
+ if ((val / NBASE) != oldval) /* possible overflow? */
+ {
+ if (!neg || (-val) != val || val == 0 || oldval < 0)
+ return false;
+ }
+ }
+
+ *result = neg ? -val : val;
+ return true;
}
-/* Convert numeric to float8; if out of range, return +/- HUGE_VAL */
-/* Caller should have eliminated the possibility of NAN */
+/*
+ * Convert int8 value to numeric.
+ */
+static void
+int8_to_numericvar(int64 val, NumericVar *var)
+{
+ uint64 uval,
+ newuval;
+ NumericDigit *ptr;
+ int ndigits;
+
+ /* int8 can require at most 19 decimal digits; add one for safety */
+ alloc_var(var, 20/DEC_DIGITS);
+ if (val < 0)
+ {
+ var->sign = NUMERIC_NEG;
+ uval = -val;
+ }
+ else
+ {
+ var->sign = NUMERIC_POS;
+ uval = val;
+ }
+ var->dscale = 0;
+ if (val == 0)
+ {
+ var->ndigits = 0;
+ var->weight = 0;
+ return;
+ }
+ ptr = var->digits + var->ndigits;
+ ndigits = 0;
+ do {
+ ptr--;
+ ndigits++;
+ newuval = uval / NBASE;
+ *ptr = uval - newuval * NBASE;
+ uval = newuval;
+ } while (uval);
+ var->digits = ptr;
+ var->ndigits = ndigits;
+ var->weight = ndigits - 1;
+}
+
+/*
+ * Convert numeric to float8; if out of range, return +/- HUGE_VAL
+ */
static double
numeric_to_double_no_overflow(Numeric num)
{
}
-/* ----------
+/*
* cmp_var() -
*
- * Compare two values on variable level
- * ----------
+ * Compare two values on variable level. We assume zeroes have been
+ * truncated to no digits.
*/
static int
cmp_var(NumericVar *var1, NumericVar *var2)
}
-/* ----------
+/*
* add_var() -
*
* Full version of add functionality on variable level (handling signs).
* result might point to one of the operands too without danger.
- * ----------
*/
static void
add_var(NumericVar *var1, NumericVar *var2, NumericVar *result)
* ----------
*/
zero_var(result);
- result->rscale = Max(var1->rscale, var2->rscale);
result->dscale = Max(var1->dscale, var2->dscale);
break;
* ----------
*/
zero_var(result);
- result->rscale = Max(var1->rscale, var2->rscale);
result->dscale = Max(var1->dscale, var2->dscale);
break;
}
-/* ----------
+/*
* sub_var() -
*
* Full version of sub functionality on variable level (handling signs).
* result might point to one of the operands too without danger.
- * ----------
*/
static void
sub_var(NumericVar *var1, NumericVar *var2, NumericVar *result)
* ----------
*/
zero_var(result);
- result->rscale = Max(var1->rscale, var2->rscale);
result->dscale = Max(var1->dscale, var2->dscale);
break;
* ----------
*/
zero_var(result);
- result->rscale = Max(var1->rscale, var2->rscale);
result->dscale = Max(var1->dscale, var2->dscale);
break;
}
-/* ----------
+/*
* mul_var() -
*
* Multiplication on variable level. Product of var1 * var2 is stored
- * in result. Accuracy of result is determined by global_rscale.
- * ----------
+ * in result. Result is rounded to no more than rscale fractional digits.
*/
static void
-mul_var(NumericVar *var1, NumericVar *var2, NumericVar *result)
+mul_var(NumericVar *var1, NumericVar *var2, NumericVar *result,
+ int rscale)
{
- NumericDigit *res_buf;
- NumericDigit *res_digits;
int res_ndigits;
- int res_weight;
int res_sign;
+ int res_weight;
+ int maxdigits;
+ int *dig;
+ int carry;
+ int maxdig;
+ int newdig;
+ NumericDigit *res_digits;
int i,
ri,
i1,
i2;
- long sum = 0;
+ /* 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;
- res_weight = var1->weight + var2->weight + 2;
- res_ndigits = var1->ndigits + var2->ndigits + 1;
+ if (var1ndigits == 0 || var2ndigits == 0)
+ {
+ /* one or both inputs is zero; so is result */
+ zero_var(result);
+ result->dscale = rscale;
+ return;
+ }
+
+ /* Determine result sign and (maximum possible) weight */
if (var1->sign == var2->sign)
res_sign = NUMERIC_POS;
else
res_sign = NUMERIC_NEG;
+ res_weight = var1->weight + var2->weight + 2;
- res_buf = digitbuf_alloc(res_ndigits);
- res_digits = res_buf;
- memset(res_digits, 0, res_ndigits);
-
- ri = res_ndigits;
- for (i1 = var1->ndigits - 1; i1 >= 0; i1--)
+ /*
+ * Determine 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.
+ */
+ res_ndigits = var1ndigits + var2ndigits + 1;
+ maxdigits = res_weight + 1 + (rscale * DEC_DIGITS) + MUL_GUARD_DIGITS;
+ if (res_ndigits > maxdigits)
{
- sum = 0;
- i = --ri;
-
- for (i2 = var2->ndigits - 1; i2 >= 0; i2--)
+ 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
{
- sum += res_digits[i] + var1->digits[i1] * var2->digits[i2];
- res_digits[i--] = sum % 10;
- sum /= 10;
+ var2ndigits -= res_ndigits - maxdigits;
+ if (var2ndigits < var1ndigits)
+ var1ndigits = var2ndigits = (var1ndigits + var2ndigits) / 2;
}
- res_digits[i] = sum;
+ res_ndigits = maxdigits;
+ Assert(res_ndigits == var1ndigits + var2ndigits + 1);
}
- i = res_weight + global_rscale + 2;
- if (i >= 0 && i < res_ndigits)
+ /*
+ * We do the arithmetic in an array "dig[]" of signed int's. Since
+ * INT_MAX is noticeably larger than NBASE*NBASE, this gives us headroom
+ * to avoid normalizing carries immediately.
+ *
+ * maxdig tracks the maximum possible value of any dig[] entry;
+ * when this threatens to exceed INT_MAX, we take the time to propagate
+ * carries. To avoid overflow in maxdig itself, it actually represents
+ * the max possible value divided by NBASE-1.
+ */
+ dig = (int *) palloc0(res_ndigits * sizeof(int));
+ maxdig = 0;
+
+ ri = res_ndigits - 1;
+ for (i1 = var1ndigits - 1; i1 >= 0; ri--, i1--)
{
- sum = (res_digits[i] > 4) ? 1 : 0;
- res_ndigits = i;
- i--;
- while (sum)
+ int var1digit = var1digits[i1];
+
+ if (var1digit == 0)
+ continue;
+
+ /* Time to normalize? */
+ maxdig += var1digit;
+ if (maxdig > INT_MAX/(NBASE-1))
+ {
+ /* Yes, do it */
+ carry = 0;
+ for (i = res_ndigits-1; i >= 0; i--)
+ {
+ newdig = dig[i] + carry;
+ if (newdig >= NBASE)
+ {
+ carry = newdig/NBASE;
+ newdig -= carry*NBASE;
+ }
+ else
+ carry = 0;
+ dig[i] = newdig;
+ }
+ Assert(carry == 0);
+ /* Reset maxdig to indicate new worst-case */
+ maxdig = 1 + var1digit;
+ }
+
+ /* Add appropriate multiple of var2 into the accumulator */
+ i = ri;
+ for (i2 = var2ndigits - 1; i2 >= 0; i2--)
{
- sum += res_digits[i];
- res_digits[i--] = sum % 10;
- sum /= 10;
+ dig[i--] += var1digit * var2digits[i2];
}
}
- while (res_ndigits > 0 && *res_digits == 0)
+ /*
+ * Now we do a final carry propagation pass to normalize the result,
+ * which we combine with storing the result digits into the output.
+ * Note that this is still done at full precision w/guard digits.
+ */
+ alloc_var(result, res_ndigits);
+ res_digits = result->digits;
+ carry = 0;
+ for (i = res_ndigits-1; i >= 0; i--)
{
- res_digits++;
- res_weight--;
- res_ndigits--;
+ newdig = dig[i] + carry;
+ if (newdig >= NBASE)
+ {
+ carry = newdig/NBASE;
+ newdig -= carry*NBASE;
+ }
+ else
+ carry = 0;
+ res_digits[i] = newdig;
}
- while (res_ndigits > 0 && res_digits[res_ndigits - 1] == 0)
- res_ndigits--;
+ Assert(carry == 0);
- if (res_ndigits == 0)
- {
- res_sign = NUMERIC_POS;
- res_weight = 0;
- }
+ pfree(dig);
- digitbuf_free(result->buf);
- result->buf = res_buf;
- result->digits = res_digits;
- result->ndigits = res_ndigits;
+ /*
+ * Finally, round the result to the requested precision.
+ */
result->weight = res_weight;
- result->rscale = global_rscale;
result->sign = res_sign;
+
+ /* Round to target rscale (and set result->dscale) */
+ round_var(result, rscale);
+
+ /* Strip leading and trailing zeroes */
+ strip_var(result);
}
-/* ----------
+/*
* div_var() -
*
- * Division on variable level. Accuracy of result is determined by
- * global_rscale.
- * ----------
+ * Division on variable level. Quotient of var1 / var2 is stored
+ * in result. Result is rounded to no more than rscale fractional digits.
*/
static void
-div_var(NumericVar *var1, NumericVar *var2, NumericVar *result)
+div_var(NumericVar *var1, NumericVar *var2, NumericVar *result,
+ int rscale)
{
- NumericDigit *res_digits;
- int res_ndigits;
+ int div_ndigits;
int res_sign;
int res_weight;
- NumericVar dividend;
- NumericVar divisor[10];
- int ndigits_tmp;
- int weight_tmp;
- int rscale_tmp;
- int ri;
+ int *div;
+ int qdigit;
+ int carry;
+ int maxdiv;
+ int newdig;
+ NumericDigit *res_digits;
+ double fdividend,
+ fdivisor,
+ fdivisorinverse,
+ fquotient;
+ int qi;
int i;
- long guess;
- long first_have;
- long first_div;
- int first_nextdigit;
- int stat = 0;
+ /* 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;
/*
- * First of all division by zero check
+ * First of all division by zero check; we must not be handed an
+ * unnormalized divisor.
*/
- ndigits_tmp = var2->ndigits + 1;
- if (ndigits_tmp == 1)
+ if (var2ndigits == 0 || var2digits[0] == 0)
elog(ERROR, "division by zero");
/*
- * Determine the result sign, weight and number of digits to calculate
- */
- if (var1->sign == var2->sign)
- res_sign = NUMERIC_POS;
- else
- res_sign = NUMERIC_NEG;
- res_weight = var1->weight - var2->weight + 1;
- res_ndigits = global_rscale + res_weight;
- if (res_ndigits <= 0)
- res_ndigits = 1;
-
- /*
* Now result zero check
*/
- if (var1->ndigits == 0)
+ if (var1ndigits == 0)
{
zero_var(result);
- result->rscale = global_rscale;
+ result->dscale = rscale;
return;
}
/*
- * Initialize local variables
+ * Determine the result sign, weight and number of digits to calculate
*/
- init_var(÷nd);
- for (i = 1; i < 10; i++)
- init_var(&divisor[i]);
+ if (var1->sign == var2->sign)
+ res_sign = NUMERIC_POS;
+ else
+ res_sign = NUMERIC_NEG;
+ res_weight = var1->weight - var2->weight + 1;
+ /* The number of accurate result digits we need to produce: */
+ div_ndigits = res_weight + 1 + (rscale + DEC_DIGITS-1)/DEC_DIGITS;
+ /* Add guard digits for roundoff error */
+ div_ndigits += DIV_GUARD_DIGITS;
+ if (div_ndigits < DIV_GUARD_DIGITS)
+ div_ndigits = DIV_GUARD_DIGITS;
+ /* Must be at least var1ndigits, too, to simplify data-loading loop */
+ if (div_ndigits < var1ndigits)
+ div_ndigits = var1ndigits;
/*
- * Make a copy of the divisor which has one leading zero digit
+ * We do the arithmetic in an array "div[]" of signed int's. Since
+ * INT_MAX is noticeably larger than NBASE*NBASE, this gives us headroom
+ * to avoid normalizing carries immediately.
+ *
+ * We start with div[] containing one zero digit followed by the
+ * dividend's digits (plus appended zeroes to reach the desired
+ * precision including guard digits). Each step of the main loop
+ * computes an (approximate) quotient digit and stores it into div[],
+ * removing one position of dividend space. A final pass of carry
+ * propagation takes care of any mistaken quotient digits.
*/
- divisor[1].ndigits = ndigits_tmp;
- divisor[1].rscale = var2->ndigits;
- divisor[1].sign = NUMERIC_POS;
- divisor[1].buf = digitbuf_alloc(ndigits_tmp);
- divisor[1].digits = divisor[1].buf;
- divisor[1].digits[0] = 0;
- memcpy(&(divisor[1].digits[1]), var2->digits, ndigits_tmp - 1);
+ div = (int *) palloc0((div_ndigits + 1) * sizeof(int));
+ for (i = 0; i < var1ndigits; i++)
+ div[i+1] = var1digits[i];
/*
- * Make a copy of the dividend
+ * We estimate each quotient digit using floating-point arithmetic,
+ * taking the first four digits of the (current) dividend and divisor.
+ * This must be float to avoid overflow.
*/
- dividend.ndigits = var1->ndigits;
- dividend.weight = 0;
- dividend.rscale = var1->ndigits;
- dividend.sign = NUMERIC_POS;
- dividend.buf = digitbuf_alloc(var1->ndigits);
- dividend.digits = dividend.buf;
- memcpy(dividend.digits, var1->digits, var1->ndigits);
+ fdivisor = (double) var2digits[0];
+ for (i = 1; i < 4; i++)
+ {
+ fdivisor *= NBASE;
+ if (i < var2ndigits)
+ fdivisor += (double) var2digits[i];
+ }
+ fdivisorinverse = 1.0 / fdivisor;
/*
- * Setup the result
+ * maxdiv tracks the maximum possible absolute value of any div[] entry;
+ * when this threatens to exceed INT_MAX, we take the time to propagate
+ * carries. To avoid overflow in maxdiv itself, it actually represents
+ * the max possible abs. value divided by NBASE-1.
*/
- digitbuf_free(result->buf);
- result->buf = digitbuf_alloc(res_ndigits + 2);
- res_digits = result->buf;
- result->digits = res_digits;
- result->ndigits = res_ndigits;
- result->weight = res_weight;
- result->rscale = global_rscale;
- result->sign = res_sign;
- res_digits[0] = 0;
-
- first_div = divisor[1].digits[1] * 10;
- if (ndigits_tmp > 2)
- first_div += divisor[1].digits[2];
+ maxdiv = 1;
- first_have = 0;
- first_nextdigit = 0;
-
- weight_tmp = 1;
- rscale_tmp = divisor[1].rscale;
-
- for (ri = 0; ri <= res_ndigits; ri++)
+ /*
+ * Outer loop computes next quotient digit, which will go into div[qi]
+ */
+ for (qi = 0; qi < div_ndigits; qi++)
{
- first_have = first_have * 10;
- if (first_nextdigit >= 0 && first_nextdigit < dividend.ndigits)
- first_have += dividend.digits[first_nextdigit];
- first_nextdigit++;
-
- guess = (first_have * 10) / first_div + 1;
- if (guess > 9)
- guess = 9;
+ /* Approximate the current dividend value */
+ fdividend = (double) div[qi];
+ for (i = 1; i < 4; i++)
+ {
+ fdividend *= NBASE;
+ if (qi+i <= div_ndigits)
+ fdividend += (double) div[qi+i];
+ }
+ /* Compute the (approximate) quotient digit */
+ fquotient = fdividend * fdivisorinverse;
+ qdigit = (fquotient >= 0.0) ? ((int) fquotient) :
+ (((int) fquotient) - 1); /* truncate towards -infinity */
- while (guess > 0)
+ if (qdigit != 0)
{
- if (divisor[guess].buf == NULL)
+ /* Do we need to normalize now? */
+ maxdiv += Abs(qdigit);
+ if (maxdiv > INT_MAX/(NBASE-1))
{
- int i;
- long sum = 0;
-
- memcpy(&divisor[guess], &divisor[1], sizeof(NumericVar));
- divisor[guess].buf = digitbuf_alloc(divisor[guess].ndigits);
- divisor[guess].digits = divisor[guess].buf;
- for (i = divisor[1].ndigits - 1; i >= 0; i--)
+ /* Yes, do it */
+ carry = 0;
+ for (i = div_ndigits; i > qi; i--)
{
- sum += divisor[1].digits[i] * guess;
- divisor[guess].digits[i] = sum % 10;
- sum /= 10;
+ newdig = div[i] + carry;
+ if (newdig < 0)
+ {
+ carry = -((-newdig-1)/NBASE) - 1;
+ newdig -= carry*NBASE;
+ }
+ else if (newdig >= NBASE)
+ {
+ carry = newdig/NBASE;
+ newdig -= carry*NBASE;
+ }
+ else
+ carry = 0;
+ div[i] = newdig;
}
+ newdig = div[qi] + carry;
+ div[qi] = newdig;
+ /*
+ * All the div[] digits except possibly div[qi] are now
+ * in the range 0..NBASE-1.
+ */
+ maxdiv = Abs(newdig) / (NBASE-1);
+ maxdiv = Max(maxdiv, 1);
+ /*
+ * Recompute the quotient digit since new info may have
+ * propagated into the top four dividend digits
+ */
+ fdividend = (double) div[qi];
+ for (i = 1; i < 4; i++)
+ {
+ fdividend *= NBASE;
+ if (qi+i <= div_ndigits)
+ fdividend += (double) div[qi+i];
+ }
+ /* Compute the (approximate) quotient digit */
+ fquotient = fdividend * fdivisorinverse;
+ qdigit = (fquotient >= 0.0) ? ((int) fquotient) :
+ (((int) fquotient) - 1); /* truncate towards -infinity */
+ maxdiv += Abs(qdigit);
}
- divisor[guess].weight = weight_tmp;
- divisor[guess].rscale = rscale_tmp;
-
- stat = cmp_abs(÷nd, &divisor[guess]);
- if (stat >= 0)
- break;
-
- guess--;
- }
+ /* Subtract off the appropriate multiple of the divisor */
+ if (qdigit != 0)
+ {
+ int istop = Min(var2ndigits, div_ndigits-qi+1);
- res_digits[ri + 1] = guess;
- if (stat == 0)
- {
- ri++;
- break;
+ for (i = 0; i < istop; i++)
+ div[qi+i] -= qdigit * var2digits[i];
+ }
}
- weight_tmp--;
- rscale_tmp++;
-
- if (guess == 0)
- continue;
-
- sub_abs(÷nd, &divisor[guess], ÷nd);
+ /*
+ * The dividend digit we are about to replace might still be nonzero.
+ * Fold it into the next digit position. We don't need to worry about
+ * overflow here since this should nearly cancel with the subtraction
+ * of the divisor.
+ */
+ div[qi+1] += div[qi] * NBASE;
- first_nextdigit = dividend.weight - weight_tmp;
- first_have = 0;
- if (first_nextdigit >= 0 && first_nextdigit < dividend.ndigits)
- first_have = dividend.digits[first_nextdigit];
- first_nextdigit++;
+ div[qi] = qdigit;
}
- result->ndigits = ri + 1;
- if (ri == res_ndigits + 1)
+ /*
+ * Approximate and store the last quotient digit (div[div_ndigits])
+ */
+ fdividend = (double) div[qi];
+ for (i = 1; i < 4; i++)
{
- int carry = (res_digits[ri] > 4) ? 1 : 0;
-
- result->ndigits = ri;
- res_digits[ri] = 0;
+ fdividend *= NBASE;
+ }
+ fquotient = fdividend * fdivisorinverse;
+ qdigit = (fquotient >= 0.0) ? ((int) fquotient) :
+ (((int) fquotient) - 1); /* truncate towards -infinity */
+ div[qi] = qdigit;
- while (carry && ri > 0)
+ /*
+ * Now we do a final carry propagation pass to normalize the result,
+ * which we combine with storing the result digits into the output.
+ * Note that this is still done at full precision w/guard digits.
+ */
+ alloc_var(result, div_ndigits+1);
+ res_digits = result->digits;
+ carry = 0;
+ for (i = div_ndigits; i >= 0; i--)
+ {
+ newdig = div[i] + carry;
+ if (newdig < 0)
{
- carry += res_digits[--ri];
- res_digits[ri] = carry % 10;
- carry /= 10;
+ carry = -((-newdig-1)/NBASE) - 1;
+ newdig -= carry*NBASE;
}
+ else if (newdig >= NBASE)
+ {
+ carry = newdig/NBASE;
+ newdig -= carry*NBASE;
+ }
+ else
+ carry = 0;
+ res_digits[i] = newdig;
}
+ Assert(carry == 0);
- while (result->ndigits > 0 && *(result->digits) == 0)
- {
- (result->digits)++;
- (result->weight)--;
- (result->ndigits)--;
- }
- while (result->ndigits > 0 && result->digits[result->ndigits - 1] == 0)
- (result->ndigits)--;
- if (result->ndigits == 0)
- result->sign = NUMERIC_POS;
+ pfree(div);
/*
- * Tidy up
+ * Finally, round the result to the requested precision.
*/
- digitbuf_free(dividend.buf);
- for (i = 1; i < 10; i++)
- digitbuf_free(divisor[i].buf);
+ result->weight = res_weight;
+ result->sign = res_sign;
+
+ /* Round to target rscale (and set result->dscale) */
+ round_var(result, rscale);
+
+ /* Strip leading and trailing zeroes */
+ strip_var(result);
}
/*
* Default scale selection for division
*
- * Returns the appropriate display scale for the division result,
- * and sets global_rscale to the result scale to use during div_var.
- *
- * Note that this must be called before div_var.
+ * Returns the appropriate result scale for the division result.
*/
static int
select_div_scale(NumericVar *var1, NumericVar *var2)
i;
NumericDigit firstdigit1,
firstdigit2;
- int res_dscale;
- int res_rscale;
+ int rscale;
/*
* The result scale of a division isn't specified in any SQL standard.
- * For PostgreSQL we select a display scale that will give at least
+ * For PostgreSQL we select a result scale that will give at least
* NUMERIC_MIN_SIG_DIGITS significant digits, so that numeric gives a
* result no less accurate than float8; but use a scale not less than
* either input's display scale.
- *
- * The result scale is NUMERIC_EXTRA_DIGITS more than the display scale,
- * to provide some guard digits in the calculation.
*/
/* Get the actual (normalized) weight and first digit of each input */
if (firstdigit1 <= firstdigit2)
qweight--;
- /* Select display scale */
- res_dscale = NUMERIC_MIN_SIG_DIGITS - qweight;
- res_dscale = Max(res_dscale, var1->dscale);
- res_dscale = Max(res_dscale, var2->dscale);
- res_dscale = Max(res_dscale, NUMERIC_MIN_DISPLAY_SCALE);
- res_dscale = Min(res_dscale, NUMERIC_MAX_DISPLAY_SCALE);
-
/* Select result scale */
- res_rscale = res_dscale + NUMERIC_EXTRA_DIGITS;
- global_rscale = res_rscale;
+ rscale = NUMERIC_MIN_SIG_DIGITS - qweight * DEC_DIGITS;
+ rscale = Max(rscale, var1->dscale);
+ rscale = Max(rscale, var2->dscale);
+ rscale = Max(rscale, NUMERIC_MIN_DISPLAY_SCALE);
+ rscale = Min(rscale, NUMERIC_MAX_DISPLAY_SCALE);
- return res_dscale;
+ return rscale;
}
-/* ----------
+/*
* mod_var() -
*
* Calculate the modulo of two numerics at variable level
- * ----------
*/
static void
mod_var(NumericVar *var1, NumericVar *var2, NumericVar *result)
{
NumericVar tmp;
- int save_global_rscale;
- int div_dscale;
+ int rscale;
init_var(&tmp);
/* ---------
* We do this using the equation
* mod(x,y) = x - trunc(x/y)*y
- * We set global_rscale the same way numeric_div and numeric_mul do
+ * We set rscale the same way numeric_div and numeric_mul do
* to get the right answer from the equation. The final result,
* however, need not be displayed to more precision than the inputs.
* ----------
*/
- save_global_rscale = global_rscale;
-
- div_dscale = select_div_scale(var1, var2);
-
- div_var(var1, var2, &tmp);
+ rscale = select_div_scale(var1, var2);
- tmp.dscale = div_dscale;
+ div_var(var1, var2, &tmp, rscale);
- /* do trunc() by forgetting digits to the right of the decimal point */
- tmp.ndigits = Max(0, Min(tmp.ndigits, tmp.weight + 1));
+ trunc_var(&tmp, 0);
- global_rscale = var2->rscale + tmp.rscale;
-
- mul_var(var2, &tmp, &tmp);
+ mul_var(var2, &tmp, &tmp, var2->dscale + tmp.dscale);
sub_var(var1, &tmp, result);
- result->dscale = Max(var1->dscale, var2->dscale);
+ round_var(result, Max(var1->dscale, var2->dscale));
- global_rscale = save_global_rscale;
free_var(&tmp);
}
-/* ----------
+/*
* ceil_var() -
*
* Return the smallest integer greater than or equal to the argument
* on variable level
- * ----------
*/
static void
ceil_var(NumericVar *var, NumericVar *result)
init_var(&tmp);
set_var_from_var(var, &tmp);
- tmp.rscale = 0;
- tmp.ndigits = Min(tmp.ndigits, Max(0, tmp.weight + 1));
- if (tmp.sign == NUMERIC_POS && cmp_var(var, &tmp) != 0)
+ trunc_var(&tmp, 0);
+
+ if (var->sign == NUMERIC_POS && cmp_var(var, &tmp) != 0)
add_var(&tmp, &const_one, &tmp);
set_var_from_var(&tmp, result);
}
-/* ----------
+/*
* floor_var() -
*
* Return the largest integer equal to or less than the argument
* on variable level
- * ----------
*/
static void
floor_var(NumericVar *var, NumericVar *result)
init_var(&tmp);
set_var_from_var(var, &tmp);
- tmp.rscale = 0;
- tmp.ndigits = Min(tmp.ndigits, Max(0, tmp.weight + 1));
- if (tmp.sign == NUMERIC_NEG && cmp_var(var, &tmp) != 0)
+ trunc_var(&tmp, 0);
+
+ if (var->sign == NUMERIC_NEG && cmp_var(var, &tmp) != 0)
sub_var(&tmp, &const_one, &tmp);
set_var_from_var(&tmp, result);
}
-/* ----------
+/*
* sqrt_var() -
*
* Compute the square root of x using Newton's algorithm
- * ----------
*/
static void
-sqrt_var(NumericVar *arg, NumericVar *result)
+sqrt_var(NumericVar *arg, NumericVar *result, int rscale)
{
NumericVar tmp_arg;
NumericVar tmp_val;
NumericVar last_val;
- int res_rscale;
- int save_global_rscale;
+ int local_rscale;
int stat;
- save_global_rscale = global_rscale;
- global_rscale += 8;
- res_rscale = global_rscale;
+ local_rscale = rscale + 8;
stat = cmp_var(arg, &const_zero);
if (stat == 0)
{
- set_var_from_var(&const_zero, result);
- result->rscale = res_rscale;
- result->sign = NUMERIC_POS;
- global_rscale = save_global_rscale;
+ zero_var(result);
+ result->dscale = rscale;
return;
}
/*
* Initialize the result to the first guess
*/
- digitbuf_free(result->buf);
- result->buf = digitbuf_alloc(1);
- result->digits = result->buf;
+ alloc_var(result, 1);
result->digits[0] = tmp_arg.digits[0] / 2;
if (result->digits[0] == 0)
result->digits[0] = 1;
- result->ndigits = 1;
result->weight = tmp_arg.weight / 2;
- result->rscale = res_rscale;
result->sign = NUMERIC_POS;
set_var_from_var(result, &last_val);
for (;;)
{
- div_var(&tmp_arg, result, &tmp_val);
+ div_var(&tmp_arg, result, &tmp_val, local_rscale);
add_var(result, &tmp_val, result);
- div_var(result, &const_two, result);
+ mul_var(result, &const_zero_point_five, result, local_rscale);
if (cmp_var(&last_val, result) == 0)
break;
free_var(&tmp_val);
free_var(&tmp_arg);
- global_rscale = save_global_rscale;
- div_var(result, &const_one, result);
+ /* Round to requested precision */
+ round_var(result, rscale);
}
-/* ----------
+/*
* exp_var() -
*
* Raise e to the power of x
- * ----------
*/
static void
-exp_var(NumericVar *arg, NumericVar *result)
+exp_var(NumericVar *arg, NumericVar *result, int rscale)
{
NumericVar x;
- NumericVar xpow;
- NumericVar ifac;
- NumericVar elem;
- NumericVar ni;
- int d;
- int i;
int xintval;
- int ndiv2 = 0;
bool xneg = FALSE;
- int save_global_rscale;
-
+ 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(&xpow);
- init_var(&ifac);
- init_var(&elem);
- init_var(&ni);
set_var_from_var(arg, &x);
x.sign = NUMERIC_POS;
}
- /* Select an appropriate scale for internal calculation */
+ /* Extract the integer part, remove it from x */
xintval = 0;
- for (i = x.weight, d = 0; i >= 0; i--, d++)
+ while (x.weight >= 0)
{
- xintval *= 10;
- if (d < x.ndigits)
- xintval += x.digits[d];
- if (xintval >= NUMERIC_MAX_RESULT_SCALE)
+ 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)
elog(ERROR, "argument for EXP() too big");
}
- save_global_rscale = global_rscale;
- global_rscale += xintval / 2 + 8;
+ /* Select an appropriate scale for internal calculation */
+ local_rscale = rscale + MUL_GUARD_DIGITS * 2;
+
+ /* Compute e^xfrac */
+ exp_var_internal(&x, result, local_rscale);
- /* Reduce input into range 0 <= x <= 0.1 */
- while (cmp_var(&x, &const_zero_point_one) > 0)
+ /* If there's an integer part, multiply by e^xint */
+ if (xintval > 0)
+ {
+ NumericVar e;
+
+ 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);
+ }
+
+ /* Compensate for input sign, and round to requested rscale */
+ if (xneg)
+ div_var(&const_one, result, result, rscale);
+ 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;
+
+ init_var(&x);
+ init_var(&xpow);
+ init_var(&ifac);
+ init_var(&elem);
+ init_var(&ni);
+
+ set_var_from_var(arg, &x);
+
+ Assert(x.sign == NUMERIC_POS);
+
+ local_rscale = rscale + 8;
+
+ /* Reduce input into range 0 <= x <= 0.01 */
+ while (cmp_var(&x, &const_zero_point_01) > 0)
{
ndiv2++;
- global_rscale++;
- div_var(&x, &const_two, &x);
+ local_rscale++;
+ mul_var(&x, &const_zero_point_five, &x, x.dscale+1);
}
/*
* exp(x) = 1 + x + x^2/2! + x^3/3! + ...
*
* Given the limited range of x, this should converge reasonably quickly.
- * We run the series until the terms fall below the global_rscale limit.
+ * 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);
for (;;)
{
add_var(&ni, &const_one, &ni);
- mul_var(&xpow, &x, &xpow);
- mul_var(&ifac, &ni, &ifac);
- div_var(&xpow, &ifac, &elem);
+ mul_var(&xpow, &x, &xpow, local_rscale);
+ mul_var(&ifac, &ni, &ifac, 0);
+ div_var(&xpow, &ifac, &elem, local_rscale);
if (elem.ndigits == 0)
break;
/* Compensate for argument range reduction */
while (ndiv2-- > 0)
- mul_var(result, result, result);
-
- /* Compensate for input sign, and round to caller's global_rscale */
- global_rscale = save_global_rscale;
-
- if (xneg)
- div_var(&const_one, result, result);
- else
- div_var(result, &const_one, result);
-
- result->sign = NUMERIC_POS;
+ mul_var(result, result, result, local_rscale);
free_var(&x);
free_var(&xpow);
}
-/* ----------
+/*
* ln_var() -
*
* Compute the natural log of x
- * ----------
*/
static void
-ln_var(NumericVar *arg, NumericVar *result)
+ln_var(NumericVar *arg, NumericVar *result, int rscale)
{
NumericVar x;
NumericVar xx;
NumericVar ni;
NumericVar elem;
NumericVar fact;
- int save_global_rscale;
+ int local_rscale;
if (cmp_var(arg, &const_zero) <= 0)
elog(ERROR, "math error on numeric - cannot compute LN of value <= zero");
- save_global_rscale = global_rscale;
- global_rscale += 8;
+ local_rscale = rscale + 8;
init_var(&x);
init_var(&xx);
init_var(&elem);
init_var(&fact);
- set_var_from_var(&const_two, &fact);
set_var_from_var(arg, &x);
+ set_var_from_var(&const_two, &fact);
/* Reduce input into range 0.9 < x < 1.1 */
while (cmp_var(&x, &const_zero_point_nine) <= 0)
{
- global_rscale++;
- sqrt_var(&x, &x);
- mul_var(&fact, &const_two, &fact);
+ local_rscale++;
+ sqrt_var(&x, &x, local_rscale);
+ mul_var(&fact, &const_two, &fact, 0);
}
while (cmp_var(&x, &const_one_point_one) >= 0)
{
- global_rscale++;
- sqrt_var(&x, &x);
- mul_var(&fact, &const_two, &fact);
+ local_rscale++;
+ sqrt_var(&x, &x, local_rscale);
+ mul_var(&fact, &const_two, &fact, 0);
}
/*
*/
sub_var(&x, &const_one, result);
add_var(&x, &const_one, &elem);
- div_var(result, &elem, result);
+ div_var(result, &elem, result, local_rscale);
set_var_from_var(result, &xx);
- mul_var(result, result, &x);
+ mul_var(result, result, &x, local_rscale);
set_var_from_var(&const_one, &ni);
for (;;)
{
add_var(&ni, &const_two, &ni);
- mul_var(&xx, &x, &xx);
- div_var(&xx, &ni, &elem);
+ mul_var(&xx, &x, &xx, local_rscale);
+ div_var(&xx, &ni, &elem, local_rscale);
if (elem.ndigits == 0)
break;
add_var(result, &elem, result);
- if (elem.weight < (result->weight - 2 * global_rscale))
+ if (elem.weight < (result->weight - local_rscale * 2/DEC_DIGITS))
break;
}
- /* Compensate for argument range reduction, round to caller's rscale */
- global_rscale = save_global_rscale;
-
- mul_var(result, &fact, result);
+ /* Compensate for argument range reduction, round to requested rscale */
+ mul_var(result, &fact, result, rscale);
free_var(&x);
free_var(&xx);
}
-/* ----------
+/*
* log_var() -
*
* Compute the logarithm of num in a given base.
*
- * Note: this routine chooses rscale and dscale of the result.
- * ----------
+ * Note: this routine chooses dscale of the result.
*/
static void
log_var(NumericVar *base, NumericVar *num, NumericVar *result)
{
NumericVar ln_base;
NumericVar ln_num;
- int save_global_rscale = global_rscale;
- int res_dscale;
+ int dec_digits;
+ int rscale;
+ int local_rscale;
init_var(&ln_base);
init_var(&ln_num);
- /* Set scale for ln() calculations */
- if (num->weight > 0)
- res_dscale = NUMERIC_MIN_SIG_DIGITS - (int) log10(num->weight);
- else if (num->weight < 0)
- res_dscale = NUMERIC_MIN_SIG_DIGITS - (int) log10(- num->weight);
+ /* 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
- res_dscale = NUMERIC_MIN_SIG_DIGITS;
+ rscale = NUMERIC_MIN_SIG_DIGITS;
- res_dscale = Max(res_dscale, base->dscale);
- res_dscale = Max(res_dscale, num->dscale);
- res_dscale = Max(res_dscale, NUMERIC_MIN_DISPLAY_SCALE);
- res_dscale = Min(res_dscale, NUMERIC_MAX_DISPLAY_SCALE);
+ rscale = Max(rscale, base->dscale);
+ rscale = Max(rscale, num->dscale);
+ rscale = Max(rscale, NUMERIC_MIN_DISPLAY_SCALE);
+ rscale = Min(rscale, NUMERIC_MAX_DISPLAY_SCALE);
- global_rscale = res_dscale + 8;
+ local_rscale = rscale + 8;
/* Form natural logarithms */
- ln_var(base, &ln_base);
- ln_var(num, &ln_num);
+ ln_var(base, &ln_base, local_rscale);
+ ln_var(num, &ln_num, local_rscale);
- ln_base.dscale = res_dscale;
- ln_num.dscale = res_dscale;
+ ln_base.dscale = rscale;
+ ln_num.dscale = rscale;
/* Select scale for division result */
- res_dscale = select_div_scale(&ln_num, &ln_base);
+ rscale = select_div_scale(&ln_num, &ln_base);
- div_var(&ln_num, &ln_base, result);
-
- result->dscale = res_dscale;
-
- global_rscale = save_global_rscale;
+ div_var(&ln_num, &ln_base, result, rscale);
free_var(&ln_num);
free_var(&ln_base);
}
-/* ----------
+/*
* power_var() -
*
* Raise base to the power of exp
*
- * Note: this routine chooses rscale and dscale of the result.
- * ----------
+ * Note: this routine chooses dscale of the result.
*/
static void
power_var(NumericVar *base, NumericVar *exp, NumericVar *result)
{
NumericVar ln_base;
NumericVar ln_num;
- int save_global_rscale = global_rscale;
- int res_dscale;
+ int dec_digits;
+ int rscale;
+ int local_rscale;
double val;
+ /* If exp can be represented as an integer, use power_var_int */
+ if (exp->ndigits == 0 || exp->ndigits <= exp->weight + 1)
+ {
+ /* exact integer, but does it fit in int? */
+ NumericVar x;
+ int64 expval64;
+
+ /* must copy because numericvar_to_int8() scribbles on input */
+ init_var(&x);
+ set_var_from_var(exp, &x);
+ if (numericvar_to_int8(&x, &expval64))
+ {
+ int expval = (int) expval64;
+
+ /* Test for overflow by reverse-conversion. */
+ if ((int64) expval == expval64)
+ {
+ /* 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, expval, result, rscale);
+
+ free_var(&x);
+ return;
+ }
+ }
+ free_var(&x);
+ }
+
init_var(&ln_base);
init_var(&ln_num);
/* Set scale for ln() calculation --- need extra accuracy here */
- if (base->weight > 0)
- res_dscale = NUMERIC_MIN_SIG_DIGITS*2 - (int) log10(base->weight);
- else if (base->weight < 0)
- res_dscale = NUMERIC_MIN_SIG_DIGITS*2 - (int) log10(- base->weight);
- else
- res_dscale = NUMERIC_MIN_SIG_DIGITS*2;
- res_dscale = Max(res_dscale, base->dscale * 2);
- res_dscale = Max(res_dscale, exp->dscale * 2);
- res_dscale = Max(res_dscale, NUMERIC_MIN_DISPLAY_SCALE);
- res_dscale = Min(res_dscale, NUMERIC_MAX_DISPLAY_SCALE);
+ /* Approx decimal digits before decimal point */
+ dec_digits = (base->weight + 1) * DEC_DIGITS;
- global_rscale = res_dscale + 8;
+ 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;
- ln_var(base, &ln_base);
+ 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);
- ln_base.dscale = res_dscale;
+ local_rscale = rscale + 8;
- mul_var(&ln_base, exp, &ln_num);
+ ln_var(base, &ln_base, local_rscale);
- ln_num.dscale = res_dscale;
+ mul_var(&ln_base, exp, &ln_num, local_rscale);
- /* Set scale for exp() */
+ /* Set scale for exp() -- compare numeric_exp() */
/* convert input to float8, ignoring overflow */
val = numericvar_to_double_no_overflow(&ln_num);
val = Max(val, -NUMERIC_MAX_RESULT_SCALE);
val = Min(val, NUMERIC_MAX_RESULT_SCALE);
- res_dscale = NUMERIC_MIN_SIG_DIGITS - (int) val;
- res_dscale = Max(res_dscale, base->dscale);
- res_dscale = Max(res_dscale, exp->dscale);
- res_dscale = Max(res_dscale, NUMERIC_MIN_DISPLAY_SCALE);
- res_dscale = Min(res_dscale, NUMERIC_MAX_DISPLAY_SCALE);
+ 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);
- global_rscale = res_dscale + 8;
+ exp_var(&ln_num, result, rscale);
- exp_var(&ln_num, result);
+ free_var(&ln_num);
+ free_var(&ln_base);
+}
- result->dscale = res_dscale;
+/*
+ * power_var_int() -
+ *
+ * Raise base to the power of exp, where exp is an integer.
+ */
+static void
+power_var_int(NumericVar *base, int exp, NumericVar *result, int rscale)
+{
+ bool neg;
+ NumericVar base_prod;
+ int local_rscale;
- global_rscale = save_global_rscale;
+ /* Detect some special cases, particularly 0^0. */
- free_var(&ln_num);
- free_var(&ln_base);
+ switch (exp)
+ {
+ case 0:
+ if (base->ndigits == 0)
+ elog(ERROR, "zero raised to zero is undefined");
+ set_var_from_var(&const_one, result);
+ result->dscale = rscale; /* no need to round */
+ return;
+ case 1:
+ set_var_from_var(base, result);
+ round_var(result, rscale);
+ return;
+ case -1:
+ div_var(&const_one, base, result, rscale);
+ return;
+ case 2:
+ mul_var(base, base, result, rscale);
+ return;
+ default:
+ break;
+ }
+
+ /*
+ * The general case repeatedly multiplies base according to the
+ * bit pattern of exp. We do the multiplications with some extra
+ * precision.
+ */
+ neg = (exp < 0);
+ exp = Abs(exp);
+
+ local_rscale = rscale + MUL_GUARD_DIGITS * 2;
+
+ init_var(&base_prod);
+ set_var_from_var(base, &base_prod);
+
+ if (exp & 1)
+ set_var_from_var(base, result);
+ else
+ set_var_from_var(&const_one, result);
+
+ while ((exp >>= 1) > 0)
+ {
+ mul_var(&base_prod, &base_prod, &base_prod, local_rscale);
+ if (exp & 1)
+ mul_var(&base_prod, result, result, local_rscale);
+ }
+
+ free_var(&base_prod);
+
+ /* Compensate for input sign, and round to requested rscale */
+ if (neg)
+ div_var(&const_one, result, result, rscale);
+ else
+ round_var(result, rscale);
}
static int
cmp_abs(NumericVar *var1, NumericVar *var2)
{
+ NumericDigit *var1digits = var1->digits;
+ NumericDigit *var2digits = var2->digits;
int i1 = 0;
int i2 = 0;
int w1 = var1->weight;
int w2 = var2->weight;
- int stat;
+
+ /* Check any digits before the first common digit */
while (w1 > w2 && i1 < var1->ndigits)
{
- if (var1->digits[i1++] != 0)
+ if (var1digits[i1++] != 0)
return 1;
w1--;
}
while (w2 > w1 && i2 < var2->ndigits)
{
- if (var2->digits[i2++] != 0)
+ if (var2digits[i2++] != 0)
return -1;
w2--;
}
+ /* At this point, either w1 == w2 or we've run out of digits */
+
if (w1 == w2)
{
while (i1 < var1->ndigits && i2 < var2->ndigits)
{
- stat = var1->digits[i1++] - var2->digits[i2++];
+ int stat = var1digits[i1++] - var2digits[i2++];
+
if (stat)
{
if (stat > 0)
}
}
+ /*
+ * At this point, we've run out of digits on one side or the other;
+ * so any remaining nonzero digits imply that side is larger
+ */
while (i1 < var1->ndigits)
{
- if (var1->digits[i1++] != 0)
+ if (var1digits[i1++] != 0)
return 1;
}
while (i2 < var2->ndigits)
{
- if (var2->digits[i2++] != 0)
+ if (var2digits[i2++] != 0)
return -1;
}
}
-/* ----------
+/*
* add_abs() -
*
* Add the absolute values of two variables into result.
* result might point to one of the operands without danger.
- * ----------
*/
static void
add_abs(NumericVar *var1, NumericVar *var2, NumericVar *result)
NumericDigit *res_digits;
int res_ndigits;
int res_weight;
- int res_rscale;
+ int res_rscale,
+ rscale1,
+ rscale2;
int res_dscale;
int i,
i1,
NumericDigit *var2digits = var2->digits;
res_weight = Max(var1->weight, var2->weight) + 1;
- res_rscale = Max(var1->rscale, var2->rscale);
+
res_dscale = Max(var1->dscale, var2->dscale);
+
+ /* Note: here we are figuring rscale in base-NBASE digits */
+ rscale1 = var1->ndigits - var1->weight - 1;
+ rscale2 = var2->ndigits - var2->weight - 1;
+ res_rscale = Max(rscale1, rscale2);
+
res_ndigits = res_rscale + res_weight + 1;
if (res_ndigits <= 0)
res_ndigits = 1;
- res_buf = digitbuf_alloc(res_ndigits);
- res_digits = res_buf;
+ res_buf = digitbuf_alloc(res_ndigits + 1);
+ res_buf[0] = 0; /* spare digit for later rounding */
+ res_digits = res_buf + 1;
i1 = res_rscale + var1->weight + 1;
i2 = res_rscale + var2->weight + 1;
if (i2 >= 0 && i2 < var2ndigits)
carry += var2digits[i2];
- if (carry >= 10)
+ if (carry >= NBASE)
{
- res_digits[i] = carry - 10;
+ res_digits[i] = carry - NBASE;
carry = 1;
}
else
Assert(carry == 0); /* else we failed to allow for carry out */
- while (res_ndigits > 0 && *res_digits == 0)
- {
- res_digits++;
- res_weight--;
- res_ndigits--;
- }
- while (res_ndigits > 0 && res_digits[res_ndigits - 1] == 0)
- res_ndigits--;
-
- if (res_ndigits == 0)
- res_weight = 0;
-
digitbuf_free(result->buf);
result->ndigits = res_ndigits;
result->buf = res_buf;
result->digits = res_digits;
result->weight = res_weight;
- result->rscale = res_rscale;
result->dscale = res_dscale;
+
+ /* Remove leading/trailing zeroes */
+ strip_var(result);
}
-/* ----------
- * sub_abs() -
+/*
+ * sub_abs()
*
* Subtract the absolute value of var2 from the absolute value of var1
* and store in result. result might point to one of the operands
* without danger.
*
* ABS(var1) MUST BE GREATER OR EQUAL ABS(var2) !!!
- * ----------
*/
static void
sub_abs(NumericVar *var1, NumericVar *var2, NumericVar *result)
NumericDigit *res_digits;
int res_ndigits;
int res_weight;
- int res_rscale;
+ int res_rscale,
+ rscale1,
+ rscale2;
int res_dscale;
int i,
i1,
NumericDigit *var2digits = var2->digits;
res_weight = var1->weight;
- res_rscale = Max(var1->rscale, var2->rscale);
+
res_dscale = Max(var1->dscale, var2->dscale);
+
+ /* Note: here we are figuring rscale in base-NBASE digits */
+ rscale1 = var1->ndigits - var1->weight - 1;
+ rscale2 = var2->ndigits - var2->weight - 1;
+ res_rscale = Max(rscale1, rscale2);
+
res_ndigits = res_rscale + res_weight + 1;
if (res_ndigits <= 0)
res_ndigits = 1;
- res_buf = digitbuf_alloc(res_ndigits);
- res_digits = res_buf;
+ res_buf = digitbuf_alloc(res_ndigits + 1);
+ res_buf[0] = 0; /* spare digit for later rounding */
+ res_digits = res_buf + 1;
i1 = res_rscale + var1->weight + 1;
i2 = res_rscale + var2->weight + 1;
if (borrow < 0)
{
- res_digits[i] = borrow + 10;
+ res_digits[i] = borrow + NBASE;
borrow = -1;
}
else
Assert(borrow == 0); /* else caller gave us var1 < var2 */
- while (res_ndigits > 0 && *res_digits == 0)
- {
- res_digits++;
- res_weight--;
- res_ndigits--;
- }
- while (res_ndigits > 0 && res_digits[res_ndigits - 1] == 0)
- res_ndigits--;
-
- if (res_ndigits == 0)
- res_weight = 0;
-
digitbuf_free(result->buf);
result->ndigits = res_ndigits;
result->buf = res_buf;
result->digits = res_digits;
result->weight = res_weight;
- result->rscale = res_rscale;
result->dscale = res_dscale;
+
+ /* Remove leading/trailing zeroes */
+ strip_var(result);
+}
+
+/*
+ * round_var
+ *
+ * Round the value of a variable to no more than rscale decimal digits
+ * after the decimal point. NOTE: we allow rscale < 0 here, implying
+ * rounding before the decimal point.
+ */
+static void
+round_var(NumericVar *var, int rscale)
+{
+ NumericDigit *digits = var->digits;
+ int di;
+ int ndigits;
+ int carry;
+
+ var->dscale = rscale;
+
+ /* decimal digits wanted */
+ di = (var->weight + 1) * DEC_DIGITS + rscale;
+
+ /*
+ * If di = 0, the value loses all digits, but could round up to 1
+ * if its first extra digit is >= 5. If di < 0 the result must be 0.
+ */
+ if (di < 0)
+ {
+ var->ndigits = 0;
+ var->weight = 0;
+ var->sign = NUMERIC_POS;
+ }
+ else
+ {
+ /* NBASE digits wanted */
+ ndigits = (di + DEC_DIGITS-1) / DEC_DIGITS;
+
+ /* 0, or number of decimal digits to keep in last NBASE digit */
+ di %= DEC_DIGITS;
+
+ if (ndigits < var->ndigits ||
+ (ndigits == var->ndigits && di > 0))
+ {
+ var->ndigits = ndigits;
+
+#if DEC_DIGITS == 1
+ /* di must be zero */
+ carry = (digits[ndigits] >= HALF_NBASE) ? 1 : 0;
+#else
+ if (di == 0)
+ {
+ carry = (digits[ndigits] >= HALF_NBASE) ? 1 : 0;
+ }
+ else
+ {
+ /* Must round within last NBASE digit */
+ int extra,
+ pow10;
+
+#if DEC_DIGITS == 4
+ pow10 = round_powers[di];
+#elif DEC_DIGITS == 2
+ pow10 = 10;
+#else
+#error unsupported NBASE
+#endif
+ extra = digits[--ndigits] % pow10;
+ digits[ndigits] -= extra;
+ carry = 0;
+ if (extra >= pow10/2)
+ {
+ pow10 += digits[ndigits];
+ if (pow10 >= NBASE)
+ {
+ pow10 -= NBASE;
+ carry = 1;
+ }
+ digits[ndigits] = pow10;
+ }
+ }
+#endif
+
+ /* Propagate carry if needed */
+ while (carry)
+ {
+ carry += digits[--ndigits];
+ if (carry >= NBASE)
+ {
+ digits[ndigits] = carry - NBASE;
+ carry = 1;
+ }
+ else
+ {
+ digits[ndigits] = carry;
+ carry = 0;
+ }
+ }
+
+ if (ndigits < 0)
+ {
+ Assert(ndigits == -1); /* better not have added > 1 digit */
+ Assert(var->digits > var->buf);
+ var->digits--;
+ var->ndigits++;
+ var->weight++;
+ }
+ }
+ }
+}
+
+/*
+ * trunc_var
+ *
+ * Truncate the value of a variable at rscale decimal digits after the
+ * decimal point. NOTE: we allow rscale < 0 here, implying
+ * truncation before the decimal point.
+ */
+static void
+trunc_var(NumericVar *var, int rscale)
+{
+ int di;
+ int ndigits;
+
+ var->dscale = rscale;
+
+ /* decimal digits wanted */
+ di = (var->weight + 1) * DEC_DIGITS + rscale;
+
+ /*
+ * If di <= 0, the value loses all digits.
+ */
+ if (di <= 0)
+ {
+ var->ndigits = 0;
+ var->weight = 0;
+ var->sign = NUMERIC_POS;
+ }
+ else
+ {
+ /* NBASE digits wanted */
+ ndigits = (di + DEC_DIGITS-1) / DEC_DIGITS;
+
+ if (ndigits <= var->ndigits)
+ {
+ var->ndigits = ndigits;
+
+#if DEC_DIGITS == 1
+ /* no within-digit stuff to worry about */
+#else
+ /* 0, or number of decimal digits to keep in last NBASE digit */
+ di %= DEC_DIGITS;
+
+ if (di > 0)
+ {
+ /* Must truncate within last NBASE digit */
+ NumericDigit *digits = var->digits;
+ int extra,
+ pow10;
+
+#if DEC_DIGITS == 4
+ pow10 = round_powers[di];
+#elif DEC_DIGITS == 2
+ pow10 = 10;
+#else
+#error unsupported NBASE
+#endif
+ extra = digits[--ndigits] % pow10;
+ digits[ndigits] -= extra;
+ }
+#endif
+ }
+ }
+}
+
+/*
+ * strip_var
+ *
+ * Strip any leading and trailing zeroes from a numeric variable
+ */
+static void
+strip_var(NumericVar *var)
+{
+ NumericDigit *digits = var->digits;
+ int ndigits = var->ndigits;
+
+ /* Strip leading zeroes */
+ while (ndigits > 0 && *digits == 0)
+ {
+ digits++;
+ var->weight--;
+ ndigits--;
+ }
+
+ /* Strip trailing zeroes */
+ while (ndigits > 0 && digits[ndigits - 1] == 0)
+ ndigits--;
+
+ /* If it's zero, normalize the sign and weight */
+ if (ndigits == 0)
+ {
+ var->sign = NUMERIC_POS;
+ var->weight = 0;
+ }
+
+ var->digits = digits;
+ var->ndigits = ndigits;
}
OSDN Git Service
