numeric.c

关系型数据库 Postgresql 6.5.2

 */
static void
div_var(NumericVar *var1, NumericVar *var2, NumericVar *result)
{
	NumericDigit *res_digits;
	int			res_ndigits;
	int			res_sign;
	int			res_weight;
	NumericVar	dividend;
	NumericVar	divisor[10];
	int			ndigits_tmp;
	int			weight_tmp;
	int			rscale_tmp;
	int			ri;
	int			i;
	long		guess;
	long		first_have;
	long		first_div;
	int			first_nextdigit;
	int			stat = 0;

	/* ----------
	 * First of all division by zero check
	 * ----------
	 */
	ndigits_tmp = var2->ndigits + 1;
	if (ndigits_tmp == 1)
	{
		free_allvars();
		elog(ERROR, "division by zero on numeric");
	}

	/* ----------
	 * 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;

	/* ----------
	 * Now result zero check
	 * ----------
	 */
	if (var1->ndigits == 0)
	{
		digitbuf_free(result->buf);
		result->buf = digitbuf_alloc(0);
		result->digits = ((NumericDigit *) (result->buf)) + sizeof(NumericDigitBuf);
		result->ndigits = 0;
		result->weight = 0;
		result->rscale = global_rscale;
		result->sign = NUMERIC_POS;
		return;
	}

	/* ----------
	 * Initialize local variables
	 * ----------
	 */
	init_var(&dividend);
	for (i = 1; i < 10; i++)
		init_var(&divisor[i]);


	/* ----------
	 * Make a copy of the divisor which has one leading zero digit
	 * ----------
	 */
	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 = ((NumericDigit *) (divisor[1].buf)) +
		sizeof(NumericDigitBuf);
	divisor[1].digits[0] = 0;
	memcpy(&(divisor[1].digits[1]), var2->digits, ndigits_tmp - 1);

	/* ----------
	 * Make a copy of the dividend
	 * ----------
	 */
	dividend.ndigits = var1->ndigits;
	dividend.weight = 0;
	dividend.rscale = var1->ndigits;
	dividend.sign = NUMERIC_POS;
	dividend.buf = digitbuf_alloc(var1->ndigits);
	dividend.digits = ((NumericDigit *) (dividend.buf)) + sizeof(NumericDigitBuf);
	memcpy(dividend.digits, var1->digits, var1->ndigits);

	/* ----------
	 * Setup the result
	 * ----------
	 */
	digitbuf_free(result->buf);
	result->buf = digitbuf_alloc(res_ndigits + 2);
	res_digits = ((NumericDigit *) (result->buf)) + sizeof(NumericDigitBuf);
	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];

	first_have = 0;
	first_nextdigit = 0;

	weight_tmp = 1;
	rscale_tmp = divisor[1].rscale;

	for (ri = 0; ri < res_ndigits + 1; ri++)
	{
		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;

		while (guess > 0)
		{
			if (divisor[guess].buf == NULL)
			{
				int			i;
				long		sum = 0;

				memcpy(&divisor[guess], &divisor[1], sizeof(NumericVar));
				divisor[guess].buf = digitbuf_alloc(divisor[guess].ndigits);
				divisor[guess].digits = ((NumericDigit *) (divisor[guess].buf) +
										 sizeof(NumericDigitBuf));
				for (i = divisor[1].ndigits - 1; i >= 0; i--)
				{
					sum += divisor[1].digits[i] * guess;
					divisor[guess].digits[i] = sum % 10;
					sum /= 10;
				}
			}

			divisor[guess].weight = weight_tmp;
			divisor[guess].rscale = rscale_tmp;

			stat = cmp_abs(&dividend, &divisor[guess]);
			if (stat >= 0)
				break;

			guess--;
		}

		res_digits[ri + 1] = guess;
		if (stat == 0)
		{
			ri++;
			break;
		}

		weight_tmp--;
		rscale_tmp++;

		if (guess == 0)
			continue;

		sub_abs(&dividend, &divisor[guess], &dividend);

		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++;
	}

	result->ndigits = ri + 1;
	if (ri == res_ndigits + 1)
	{
		long		carry = (res_digits[ri] > 4) ? 1 : 0;

		result->ndigits = ri;
		res_digits[ri] = 0;

		while (carry && ri > 0)
		{
			carry += res_digits[--ri];
			res_digits[ri] = carry % 10;
			carry /= 10;
		}
	}

	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;

	/*
	 * Tidy up
	 *
	 */
	digitbuf_free(dividend.buf);
	for (i = 1; i < 10; i++)
		digitbuf_free(divisor[i].buf);
}


/* ----------
 * 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;

	init_var(&tmp);

	/* ----------
	 * We do it by fiddling around with global_rscale and truncating
	 * the result of the division.
	 * ----------
	 */
	save_global_rscale = global_rscale;
	global_rscale = var2->rscale + 2;

	div_var(var1, var2, &tmp);
	tmp.rscale = var2->rscale;
	tmp.ndigits = MAX(0, MIN(tmp.ndigits, tmp.weight + tmp.rscale + 1));

	global_rscale = var2->rscale;
	mul_var(var2, &tmp, &tmp);

	sub_var(var1, &tmp, result);

	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)
{
	NumericVar	tmp;

	init_var(&tmp);
	set_var_from_var(var, &tmp);

	tmp.rscale = 0;
	tmp.ndigits = MAX(0, tmp.weight + 1);
	if (tmp.sign == NUMERIC_POS && cmp_var(var, &tmp) != 0)
		add_var(&tmp, &const_one, &tmp);

	set_var_from_var(&tmp, result);
	free_var(&tmp);
}


/* ----------
 * floor_var() -
 *
 *	Return the largest integer equal to or less than the argument
 *	on variable level
 * ----------
 */
static void
floor_var(NumericVar *var, NumericVar *result)
{
	NumericVar	tmp;

	init_var(&tmp);
	set_var_from_var(var, &tmp);

	tmp.rscale = 0;
	tmp.ndigits = MAX(0, tmp.weight + 1);
	if (tmp.sign == NUMERIC_NEG && cmp_var(var, &tmp) != 0)
		sub_var(&tmp, &const_one, &tmp);

	set_var_from_var(&tmp, result);
	free_var(&tmp);
}


/* ----------
 * sqrt_var() -
 *
 *	Compute the square root of x using Newtons algorithm
 * ----------
 */
static void
sqrt_var(NumericVar *arg, NumericVar *result)
{
	NumericVar	tmp_arg;
	NumericVar	tmp_val;
	NumericVar	last_val;
	int			res_rscale;
	int			save_global_rscale;
	int			stat;

	save_global_rscale = global_rscale;
	global_rscale += 8;
	res_rscale = global_rscale;

	stat = cmp_var(arg, &const_zero);
	if (stat == 0)
	{
		set_var_from_var(&const_zero, result);
		result->rscale = res_rscale;
		result->sign = NUMERIC_POS;
		return;
	}

	if (stat < 0)
	{
		free_allvars();
		elog(ERROR, "math error on numeric - cannot compute SQRT of negative value");
	}

	init_var(&tmp_arg);
	init_var(&tmp_val);
	init_var(&last_val);

	set_var_from_var(arg, &tmp_arg);
	set_var_from_var(result, &last_val);

	/* ----------
	 * Initialize the result to the first guess
	 * ----------
	 */
	digitbuf_free(result->buf);
	result->buf = digitbuf_alloc(1);
	result->digits = ((NumericDigit *) (result->buf)) + sizeof(NumericDigitBuf);
	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;

	for (;;)
	{
		div_var(&tmp_arg, result, &tmp_val);

		add_var(result, &tmp_val, result);
		div_var(result, &const_two, result);

		if (cmp_var(&last_val, result) == 0)
			break;
		set_var_from_var(result, &last_val);
	}

	free_var(&last_val);
	free_var(&tmp_val);
	free_var(&tmp_arg);

	global_rscale = save_global_rscale;
	div_var(result, &const_one, result);
}


/* ----------
 * exp_var() -
 *
 *	Raise e to the power of x
 * ----------
 */
static void
exp_var(NumericVar *arg, NumericVar *result)
{
	NumericVar	x;
	NumericVar	xpow;
	NumericVar	ifac;
	NumericVar	elem;
	NumericVar	ni;
	int			d;
	int			i;
	int			ndiv2 = 0;
	bool		xneg = FALSE;
	int			save_global_rscale;

	init_var(&x);
	init_var(&xpow);
	init_var(&ifac);
	init_var(&elem);
	init_var(&ni);

	set_var_from_var(arg, &x);

	if (x.sign == NUMERIC_NEG)
	{
		xneg = TRUE;
		x.sign = NUMERIC_POS;
	}

	save_global_rscale = global_rscale;
	global_rscale = 0;
	for (i = x.weight, d = 0; i >= 0; i--, d++)
	{
		global_rscale *= 10;
		if (d < x.ndigits)
			global_rscale += x.digits[d];
		if (global_rscale >= 1000)
		{
			free_allvars();
			elog(ERROR, "argument for EXP() too big");
		}
	}

	global_rscale = global_rscale / 2 + save_global_rscale + 8;

	while (cmp_var(&x, &const_one) > 0)
	{
		ndiv2++;
		global_rscale++;
		div_var(&x, &const_two, &x);
	}

	add_var(&const_one, &x, result);
	set_var_from_var(&x, &xpow);
	set_var_from_var(&const_one, &ifac);
	set_var_from_var(&const_one, &ni);

	for (i = 2; TRUE; i++)
	{
		add_var(&ni, &const_one, &ni);
		mul_var(&xpow, &x, &xpow);
		mul_var(&ifac, &ni, &ifac);
		div_var(&xpow, &ifac, &elem);

		if (elem.ndigits == 0)
			break;

		add_var(result, &elem, result);
	}

	while (ndiv2-- > 0)
		mul_var(result, result, result);

	global_rscale = save_global_rscale;
	if (xneg)
		div_var(&const_one, result, result);
	else
		div_var(result, &const_one, result);

	result->sign = NUMERIC_POS;

	free_var(&x);
	free_var(&xpow);
	free_var(&ifac);
	free_var(&elem);
	free_var(&ni);
}


/* ----------
 * ln_var() -
 *
 *	Compute the natural log of x
 * ----------
 */
static void
ln_var(NumericVar *arg, NumericVar *result)
{
	NumericVar	x;
	NumericVar	xx;
	NumericVar	ni;
	NumericVar	elem;
	NumericVar	fact;
	int			i;
	int			save_global_rscale;

	if (cmp_var(arg, &const_zero) <= 0)
	{
		free_allvars();
		elog(ERROR, "math error on numeric - cannot compute LN of value <= zero");
	}

	save_global_rscale = global_rscale;
	global_rscale += 8;

	init_var(&x);
	init_var(&xx);
	init_var(&ni);
	init_var(&elem);
	init_var(&fact);

	set_var_from_var(&const_two, &fact);
	set_var_from_var(arg, &x);

	while (cmp_var(&x, &const_two) >= 0)
	{
		sqrt_var(&x, &x);
		mul_var(&fact, &const_two, &fact);
	}
	set_var_from_str("0.5", &elem);
	while (cmp_var(&x, &elem) <= 0)
	{
		sqrt_var(&x, &x);
		mul_var(&fact, &const_two, &fact);
	}

	sub_var(&x, &const_one, result);
	add_var(&x, &const_one, &elem);
	div_var(result, &elem, result);
	set_var_from_var(result, &xx);
	mul_var(result, result, &x);

	set_var_from_var(&const_one, &ni);

	for (i = 2; TRUE; i++)
	{
		add_var(&ni, &const_two, &ni);
		mul_var(&xx, &x, &xx);
		div_var(&xx, &ni, &elem);

		if (cmp_var(&elem, &const_zero) == 0)
			break;

		add_var(result, &elem, result);
	}

	global_rscale = save_global_rscale;
	mul_var(result, &fact, result);

	free_var(&x);
	free_var(&xx);
	free_var(&ni);
	free_var(&elem);
	free_var(&fact);
}


/* ----------
 * log_var() -
 *
 *	Compute the logarithm of x in a given base
 * ----------
 */
static void
log_var(NumericVar *base, NumericVar *num, NumericVar *result)
{
	NumericVar	ln_base;
	NumericVar	ln_num;

	global_rscale += 8;

	init_var(&ln_base);
	init_var(&ln_num);

	ln_var(base, &ln_base);
	ln_var(num, &ln_num);

	global_rscale -= 8;

	div_var(&ln_num, &ln_base, result);

	free_var(&ln_num);
	free_var(&ln_base);
}


/* ----------
 * power_var() -
 *
 *	Raise base to the power of exp
 * ----------
 */
static void
power_var(NumericVar *base, NumericVar *exp, NumericVar *result)
{
	NumericVar	ln_base;
	NumericVar	ln_num;
	int			save_global_rscale;

	save_global_rscale = global_rscale;
	global_rscale += global_rscale / 3 + 8;

	init_var(&ln_base);
	init_var(&ln_num);

	ln_var(base, &ln_base);
	mul_var(&ln_base, exp, &ln_num);

	global_rscale = save_global_rscale;

	exp_var(&ln_num, result);

	free_var(&ln_num);
	free_var(&ln_base);

}


/* ----------------------------------------------------------------------
 *
 * Following are the lowest level functions that operate unsigned
 * on the variable level
 *
 * ----------------------------------------------------------------------
 */


/* ----------
 * cmp_abs() -
 *
 *	Compare the absolute values of var1 and var2
 *	Returns:	-1 for ABS(var1) < ABS(var2)
 *				0  for ABS(var1) == ABS(var2)
 *				1  for ABS(var1) > ABS(var2)
 * ----------
 */
static int
cmp_abs(NumericVar *var1, NumericVar *var2)
{
	int			i1 = 0;
	int			i2 = 0;
	int			w1 = var1->weight;
	int			w2 = var2->weight;
	int			stat;

	while (w1 > w2 && i1 < var1->ndigits)
	{
		if (var1->digits[i1++] != 0)
			return 1;
		w1--;
	}
	while (w2 > w1 && i2 < var2->ndigits)
	{
		if (var2->digits[i2++] != 0)
			return -1;
		w2--;
	}

	if (w1 == w2)
	{
		while (i1 < var1->ndigits && i2 < var2->ndigits)
		{
			stat = var1->digits[i1++] - var2->digits[i2++];
			if (stat)
			{
				if (stat > 0)
					return 1;
				return -1;
			}
		}
	}

	while (i1 < var1->ndigits)
	{
		if (var1->digits[i1++] != 0)
			return 1;
	}
	while (i2 < var2->ndigits)
	{
		if (var2->digits[i2++] != 0)
			return -1;
	}

	return 0;
}


/* ----------
 * 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)
{
	NumericDigitBuf *res_buf;
	NumericDigit *res_digits;
	int			res_ndigits;
	int			res_weight;
	int			res_rscale;
	int			res_dscale;
	int			i,
				i1,
				i2;
	int			carry = 0;

	res_weight = MAX(var1->weight, var2->weight) + 1;
	res_rsca