Re-implement division for numeric values using the traditional "schoolbook"

algorithm.  This is a good deal slower than our old roundoff-error-prone
code for long inputs, so we keep the old code for use in the transcendental
functions, where everything is approximate anyway.  Also create a
user-accessible function div(numeric, numeric) to provide access to the
exact result of trunc(x/y) --- since the regular numeric / operator will
round off its result, simply computing that expression in SQL doesn't
reliably give the desired answer.  This fixes bug #3387 and various related
corner cases, and improves the usefulness of PG for high-precision integer
arithmetic.

6 files changed, 457 insertions, 31 deletions

diff --git a/src/backend/utils/adt/numeric.c b/src/backend/utils/adt/numeric.c
index 2ebfee52e45..86765d5d532 100644
--- a/src/backend/utils/adt/numeric.c
+++ b/src/backend/utils/adt/numeric.c
@@ -14,7 +14,7 @@
  * Copyright (c) 1998-2008, PostgreSQL Global Development Group
  *
  * IDENTIFICATION
- *	  $PostgreSQL: pgsql/src/backend/utils/adt/numeric.c,v 1.108 2008/01/01 19:45:52 momjian Exp $
+ *	  $PostgreSQL: pgsql/src/backend/utils/adt/numeric.c,v 1.109 2008/04/04 18:45:36 tgl Exp $
  *
  *-------------------------------------------------------------------------
  */
@@ -53,7 +53,7 @@
  * 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
+ * sqrt(INT_MAX), so that there is "headroom" for mul_var and div_var_fast to
  * postpone processing carries.
  * ----------
  */
@@ -90,6 +90,10 @@ typedef int16 NumericDigit;
 
 
 /* ----------
+ * NumericVar is the format we use for arithmetic.  The digit-array part
+ * is the same as the NumericData storage format, but the header is more
+ * complex.
+ *
  * The value represented by a NumericVar is determined by the sign, weight,
  * ndigits, and digits[] array.
  * Note: the first digit of a NumericVar's value is assumed to be multiplied
@@ -100,7 +104,7 @@ typedef int16 NumericDigit;
  * NumericVar.	digits points at the first digit in actual use (the one
  * 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
+ * a carry out of the top digit without reallocating space.  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.)
@@ -265,6 +269,8 @@ static void mul_var(NumericVar *var1, NumericVar *var2, NumericVar *result,
 		int rscale);
 static void div_var(NumericVar *var1, NumericVar *var2, NumericVar *result,
 		int rscale, bool round);
+static void div_var_fast(NumericVar *var1, NumericVar *var2, NumericVar *result,
+		int rscale, bool round);
 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);
@@ -1420,6 +1426,52 @@ numeric_div(PG_FUNCTION_ARGS)
 
 
 /*
+ * numeric_div_trunc() -
+ *
+ *	Divide one numeric into another, truncating the result to an integer
+ */
+Datum
+numeric_div_trunc(PG_FUNCTION_ARGS)
+{
+	Numeric		num1 = PG_GETARG_NUMERIC(0);
+	Numeric		num2 = PG_GETARG_NUMERIC(1);
+	NumericVar	arg1;
+	NumericVar	arg2;
+	NumericVar	result;
+	Numeric		res;
+
+	/*
+	 * Handle NaN
+	 */
+	if (NUMERIC_IS_NAN(num1) || NUMERIC_IS_NAN(num2))
+		PG_RETURN_NUMERIC(make_result(&const_nan));
+
+	/*
+	 * Unpack the arguments
+	 */
+	init_var(&arg1);
+	init_var(&arg2);
+	init_var(&result);
+
+	set_var_from_num(num1, &arg1);
+	set_var_from_num(num2, &arg2);
+
+	/*
+	 * Do the divide and return the result
+	 */
+	div_var(&arg1, &arg2, &result, 0, false);
+
+	res = make_result(&result);
+
+	free_var(&arg1);
+	free_var(&arg2);
+	free_var(&result);
+
+	PG_RETURN_NUMERIC(res);
+}
+
+
+/*
  * numeric_mod() -
  *
  *	Calculate the modulo of two numerics
@@ -4036,14 +4088,293 @@ mul_var(NumericVar *var1, NumericVar *var2, NumericVar *result,
 /*
  * div_var() -
  *
- *	Division on variable level. Quotient of var1 / var2 is stored
- *	in result.	Result is rounded to no more than rscale fractional digits.
+ *	Division on variable level. Quotient of var1 / var2 is stored in result.
+ *	The quotient is figured to exactly rscale fractional digits.
+ *	If round is true, it is rounded at the rscale'th digit; if false, it
+ *	is truncated (towards zero) at that digit.
  */
 static void
 div_var(NumericVar *var1, NumericVar *var2, NumericVar *result,
 		int rscale, bool round)
 {
 	int			div_ndigits;
+	int			res_ndigits;
+	int			res_sign;
+	int			res_weight;
+	int			carry;
+	int			borrow;
+	int			divisor1;
+	int			divisor2;
+	NumericDigit *dividend;
+	NumericDigit *divisor;
+	NumericDigit *res_digits;
+	int			i;
+	int			j;
+
+	/* copy these values into local vars for speed in inner loop */
+	int			var1ndigits = var1->ndigits;
+	int			var2ndigits = var2->ndigits;
+
+	/*
+	 * First of all division by zero check; we must not be handed an
+	 * unnormalized divisor.
+	 */
+	if (var2ndigits == 0 || var2->digits[0] == 0)
+		ereport(ERROR,
+				(errcode(ERRCODE_DIVISION_BY_ZERO),
+				 errmsg("division by zero")));
+
+	/*
+	 * Now result zero check
+	 */
+	if (var1ndigits == 0)
+	{
+		zero_var(result);
+		result->dscale = rscale;
+		return;
+	}
+
+	/*
+	 * Determine the result sign, weight and number of digits to calculate.
+	 * The weight figured here is correct if the emitted quotient has no
+	 * leading zero digits; otherwise strip_var() will fix things up.
+	 */
+	if (var1->sign == var2->sign)
+		res_sign = NUMERIC_POS;
+	else
+		res_sign = NUMERIC_NEG;
+	res_weight = var1->weight - var2->weight;
+	/* The number of accurate result digits we need to produce: */
+	res_ndigits = res_weight + 1 + (rscale + DEC_DIGITS - 1) / DEC_DIGITS;
+	/* ... but always at least 1 */
+	res_ndigits = Max(res_ndigits, 1);
+	/* If rounding needed, figure one more digit to ensure correct result */
+	if (round)
+		res_ndigits++;
+	/*
+	 * The working dividend normally requires res_ndigits + var2ndigits
+	 * digits, but make it at least var1ndigits so we can load all of var1
+	 * into it.  (There will be an additional digit dividend[0] in the
+	 * dividend space, but for consistency with Knuth's notation we don't
+	 * count that in div_ndigits.)
+	 */
+	div_ndigits = res_ndigits + var2ndigits;
+	div_ndigits = Max(div_ndigits, var1ndigits);
+
+	/*
+	 * We need a workspace with room for the working dividend (div_ndigits+1
+	 * digits) plus room for the possibly-normalized divisor (var2ndigits
+	 * digits).  It is convenient also to have a zero at divisor[0] with
+	 * the actual divisor data in divisor[1 .. var2ndigits].  Transferring the
+	 * digits into the workspace also allows us to realloc the result (which
+	 * might be the same as either input var) before we begin the main loop.
+	 * Note that we use palloc0 to ensure that divisor[0], dividend[0], and
+	 * any additional dividend positions beyond var1ndigits, start out 0.
+	 */
+	dividend = (NumericDigit *)
+		palloc0((div_ndigits + var2ndigits + 2) * sizeof(NumericDigit));
+	divisor = dividend + (div_ndigits + 1);
+	memcpy(dividend + 1, var1->digits, var1ndigits * sizeof(NumericDigit));
+	memcpy(divisor + 1, var2->digits, var2ndigits * sizeof(NumericDigit));
+
+	/*
+	 * Now we can realloc the result to hold the generated quotient digits.
+	 */
+	alloc_var(result, res_ndigits);
+	res_digits = result->digits;
+
+	if (var2ndigits == 1)
+	{
+		/*
+		 * If there's only a single divisor digit, we can use a fast path
+		 * (cf. Knuth section 4.3.1 exercise 16).
+		 */
+		divisor1 = divisor[1];
+		carry = 0;
+		for (i = 0; i < res_ndigits; i++)
+		{
+			carry = carry * NBASE + dividend[i + 1];
+			res_digits[i] = carry / divisor1;
+			carry = carry % divisor1;
+		}
+	}
+	else
+	{
+		/*
+		 * The full multiple-place algorithm is taken from Knuth volume 2,
+		 * Algorithm 4.3.1D.
+		 *
+		 * We need the first divisor digit to be >= NBASE/2.  If it isn't,
+		 * make it so by scaling up both the divisor and dividend by the
+		 * factor "d".  (The reason for allocating dividend[0] above is to
+		 * leave room for possible carry here.)
+		 */
+		if (divisor[1] < HALF_NBASE)
+		{
+			int		d = NBASE / (divisor[1] + 1);
+
+			carry = 0;
+			for (i = var2ndigits; i > 0; i--)
+			{
+				carry += divisor[i] * d;
+				divisor[i] = carry % NBASE;
+				carry = carry / NBASE;
+			}
+			Assert(carry == 0);
+			carry = 0;
+			/* at this point only var1ndigits of dividend can be nonzero */
+			for (i = var1ndigits; i >= 0; i--)
+			{
+				carry += dividend[i] * d;
+				dividend[i] = carry % NBASE;
+				carry = carry / NBASE;
+			}
+			Assert(carry == 0);
+			Assert(divisor[1] >= HALF_NBASE);
+		}
+		/* First 2 divisor digits are used repeatedly in main loop */
+		divisor1 = divisor[1];
+		divisor2 = divisor[2];
+
+		/*
+		 * Begin the main loop.  Each iteration of this loop produces the
+		 * j'th quotient digit by dividing dividend[j .. j + var2ndigits]
+		 * by the divisor; this is essentially the same as the common manual
+		 * procedure for long division.
+		 */
+		for (j = 0; j < res_ndigits; j++)
+		{
+			/* Estimate quotient digit from the first two dividend digits */
+			int		next2digits = dividend[j] * NBASE + dividend[j+1];
+			int		qhat;
+
+			/*
+			 * If next2digits are 0, then quotient digit must be 0 and there's
+			 * no need to adjust the working dividend.  It's worth testing
+			 * here to fall out ASAP when processing trailing zeroes in
+			 * a dividend.
+			 */
+			if (next2digits == 0)
+			{
+				res_digits[j] = 0;
+				continue;
+			}
+
+			if (dividend[j] == divisor1)
+				qhat = NBASE - 1;
+			else
+				qhat = next2digits / divisor1;
+			/*
+			 * Adjust quotient digit if it's too large.  Knuth proves that
+			 * after this step, the quotient digit will be either correct
+			 * or just one too large.  (Note: it's OK to use dividend[j+2]
+			 * here because we know the divisor length is at least 2.)
+			 */
+			while (divisor2 * qhat >
+				   (next2digits - qhat * divisor1) * NBASE + dividend[j+2])
+				qhat--;
+
+			/* As above, need do nothing more when quotient digit is 0 */
+			if (qhat > 0)
+			{
+				/*
+				 * Multiply the divisor by qhat, and subtract that from the
+				 * working dividend.  "carry" tracks the multiplication,
+				 * "borrow" the subtraction (could we fold these together?)
+				 */
+				carry = 0;
+				borrow = 0;
+				for (i = var2ndigits; i >= 0; i--)
+				{
+					carry += divisor[i] * qhat;
+					borrow -= carry % NBASE;
+					carry = carry / NBASE;
+					borrow += dividend[j + i];
+					if (borrow < 0)
+					{
+						dividend[j + i] = borrow + NBASE;
+						borrow = -1;
+					}
+					else
+					{
+						dividend[j + i] = borrow;
+						borrow = 0;
+					}
+				}
+				Assert(carry == 0);
+
+				/*
+				 * If we got a borrow out of the top dividend digit, then
+				 * indeed qhat was one too large.  Fix it, and add back the
+				 * divisor to correct the working dividend.  (Knuth proves
+				 * that this will occur only about 3/NBASE of the time; hence,
+				 * it's a good idea to test this code with small NBASE to be
+				 * sure this section gets exercised.)
+				 */
+				if (borrow)
+				{
+					qhat--;
+					carry = 0;
+					for (i = var2ndigits; i >= 0; i--)
+					{
+						carry += dividend[j + i] + divisor[i];
+						if (carry >= NBASE)
+						{
+							dividend[j + i] = carry - NBASE;
+							carry = 1;
+						}
+						else
+						{
+							dividend[j + i] = carry;
+							carry = 0;
+						}
+					}
+					/* A carry should occur here to cancel the borrow above */
+					Assert(carry == 1);
+				}
+			}
+
+			/* And we're done with this quotient digit */
+			res_digits[j] = qhat;
+		}
+	}
+
+	pfree(dividend);
+
+	/*
+	 * Finally, round or truncate the result to the requested precision.
+	 */
+	result->weight = res_weight;
+	result->sign = res_sign;
+
+	/* Round or truncate to target rscale (and set result->dscale) */
+	if (round)
+		round_var(result, rscale);
+	else
+		trunc_var(result, rscale);
+
+	/* Strip leading and trailing zeroes */
+	strip_var(result);
+}
+
+
+/*
+ * div_var_fast() -
+ *
+ *	This has the same API as div_var, but is implemented using the division
+ *	algorithm from the "FM" library, rather than Knuth's schoolbook-division
+ *	approach.  This is significantly faster but can produce inaccurate
+ *	results, because it sometimes has to propagate rounding to the left,
+ *	and so we can never be entirely sure that we know the requested digits
+ *	exactly.  We compute DIV_GUARD_DIGITS extra digits, but there is
+ *	no certainty that that's enough.  We use this only in the transcendental
+ *	function calculation routines, where everything is approximate anyway.
+ */
+static void
+div_var_fast(NumericVar *var1, NumericVar *var2, NumericVar *result,
+			 int rscale, bool round)
+{
+	int			div_ndigits;
 	int			res_sign;
 	int			res_weight;
 	int		   *div;
@@ -4367,30 +4698,21 @@ static void
 mod_var(NumericVar *var1, NumericVar *var2, NumericVar *result)
 {
 	NumericVar	tmp;
-	int			rscale;
 
 	init_var(&tmp);
 
 	/* ---------
 	 * We do this using the equation
 	 *		mod(x,y) = x - trunc(x/y)*y
-	 * 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.
+	 * div_var can be persuaded to give us trunc(x/y) directly.
 	 * ----------
 	 */
-	rscale = select_div_scale(var1, var2);
+	div_var(var1, var2, &tmp, 0, false);
 
-	div_var(var1, var2, &tmp, rscale, false);
-
-	trunc_var(&tmp, 0);
-
-	mul_var(var2, &tmp, &tmp, var2->dscale + tmp.dscale);
+	mul_var(var2, &tmp, &tmp, var2->dscale);
 
 	sub_var(var1, &tmp, result);
 
-	round_var(result, Max(var1->dscale, var2->dscale));
-
 	free_var(&tmp);
 }
 
@@ -4497,7 +4819,7 @@ sqrt_var(NumericVar *arg, NumericVar *result, int rscale)
 
 	for (;;)
 	{
-		div_var(&tmp_arg, result, &tmp_val, local_rscale, true);
+		div_var_fast(&tmp_arg, result, &tmp_val, local_rscale, true);
 
 		add_var(result, &tmp_val, result);
 		mul_var(result, &const_zero_point_five, result, local_rscale);
@@ -4587,7 +4909,7 @@ exp_var(NumericVar *arg, NumericVar *result, int rscale)
 
 	/* Compensate for input sign, and round to requested rscale */
 	if (xneg)
-		div_var(&const_one, result, result, rscale, true);
+		div_var_fast(&const_one, result, result, rscale, true);
 	else
 		round_var(result, rscale);
 
@@ -4652,7 +4974,7 @@ exp_var_internal(NumericVar *arg, NumericVar *result, int rscale)
 		add_var(&ni, &const_one, &ni);
 		mul_var(&xpow, &x, &xpow, local_rscale);
 		mul_var(&ifac, &ni, &ifac, 0);
-		div_var(&xpow, &ifac, &elem, local_rscale, true);
+		div_var_fast(&xpow, &ifac, &elem, local_rscale, true);
 
 		if (elem.ndigits == 0)
 			break;
@@ -4736,7 +5058,7 @@ ln_var(NumericVar *arg, NumericVar *result, int rscale)
 	 */
 	sub_var(&x, &const_one, result);
 	add_var(&x, &const_one, &elem);
-	div_var(result, &elem, result, local_rscale, true);
+	div_var_fast(result, &elem, result, local_rscale, true);
 	set_var_from_var(result, &xx);
 	mul_var(result, result, &x, local_rscale);
 
@@ -4746,7 +5068,7 @@ ln_var(NumericVar *arg, NumericVar *result, int rscale)
 	{
 		add_var(&ni, &const_two, &ni);
 		mul_var(&xx, &x, &xx, local_rscale);
-		div_var(&xx, &ni, &elem, local_rscale, true);
+		div_var_fast(&xx, &ni, &elem, local_rscale, true);
 
 		if (elem.ndigits == 0)
 			break;
@@ -4816,7 +5138,7 @@ log_var(NumericVar *base, NumericVar *num, NumericVar *result)
 	/* Select scale for division result */
 	rscale = select_div_scale(&ln_num, &ln_base);
 
-	div_var(&ln_num, &ln_base, result, rscale, true);
+	div_var_fast(&ln_num, &ln_base, result, rscale, true);
 
 	free_var(&ln_num);
 	free_var(&ln_base);
@@ -4990,7 +5312,7 @@ power_var_int(NumericVar *base, int exp, NumericVar *result, int rscale)
 
 	/* Compensate for input sign, and round to requested rscale */
 	if (neg)
-		div_var(&const_one, result, result, rscale, true);
+		div_var_fast(&const_one, result, result, rscale, true);
 	else
 		round_var(result, rscale);
 }
@@ -5361,8 +5683,8 @@ round_var(NumericVar *var, int rscale)
 /*
  * trunc_var
  *
- * Truncate the value of a variable at rscale decimal digits after the
- * decimal point.  NOTE: we allow rscale < 0 here, implying
+ * Truncate (towards zero) 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
diff --git a/src/include/catalog/catversion.h b/src/include/catalog/catversion.h
index 009cf8abf66..f4821984e91 100644
--- a/src/include/catalog/catversion.h
+++ b/src/include/catalog/catversion.h
@@ -37,7 +37,7 @@
  * Portions Copyright (c) 1996-2008, PostgreSQL Global Development Group
  * Portions Copyright (c) 1994, Regents of the University of California
  *
- * $PostgreSQL: pgsql/src/include/catalog/catversion.h,v 1.444 2008/03/23 00:24:19 tgl Exp $
+ * $PostgreSQL: pgsql/src/include/catalog/catversion.h,v 1.445 2008/04/04 18:45:36 tgl Exp $
  *
  *-------------------------------------------------------------------------
  */
@@ -53,6 +53,6 @@
  */
 
 /*							yyyymmddN */
-#define CATALOG_VERSION_NO	200803222
+#define CATALOG_VERSION_NO	200804041
 
 #endif
diff --git a/src/include/catalog/pg_proc.h b/src/include/catalog/pg_proc.h
index fba36a018d3..2abaeea211e 100644
--- a/src/include/catalog/pg_proc.h
+++ b/src/include/catalog/pg_proc.h
@@ -7,7 +7,7 @@
  * Portions Copyright (c) 1996-2008, PostgreSQL Global Development Group
  * Portions Copyright (c) 1994, Regents of the University of California
  *
- * $PostgreSQL: pgsql/src/include/catalog/pg_proc.h,v 1.486 2008/04/04 16:57:21 momjian Exp $
+ * $PostgreSQL: pgsql/src/include/catalog/pg_proc.h,v 1.487 2008/04/04 18:45:36 tgl Exp $
  *
  * NOTES
  *	  The script catalog/genbki.sh reads this file and generates .bki
@@ -1115,7 +1115,7 @@ DESCR("does not match LIKE expression");
 DATA(insert OID =  860 (  bpchar		   PGNSP PGUID 12 1 0 f f t f i 1 1042 "18" _null_ _null_ _null_	char_bpchar - _null_ _null_ ));
 DESCR("convert char to char()");
 
-DATA(insert OID = 861 ( current_database	   PGNSP PGUID 12 1 0 f f t f i 0 19 "" _null_ _null_ _null_ current_database - _null_ _null_ ));
+DATA(insert OID = 861 ( current_database	   PGNSP PGUID 12 1 0 f f t f s 0 19 "" _null_ _null_ _null_ current_database - _null_ _null_ ));
 DESCR("returns the current database");
 DATA(insert OID = 817 (  current_query        PGNSP PGUID 12 1 0 f f f f v 0 25  "" _null_ _null_ _null_  current_query - _null_ _null_ ));
 DESCR("returns the currently executing query");
@@ -2573,6 +2573,10 @@ DATA(insert OID = 1745 ( float4					PGNSP PGUID 12 1 0 f f t f i 1 700 "1700" _n
 DESCR("(internal)");
 DATA(insert OID = 1746 ( float8					PGNSP PGUID 12 1 0 f f t f i 1 701 "1700" _null_ _null_ _null_	numeric_float8 - _null_ _null_ ));
 DESCR("(internal)");
+DATA(insert OID = 1973 ( div					PGNSP PGUID 12 1 0 f f t f i 2 1700 "1700 1700" _null_ _null_ _null_	numeric_div_trunc - _null_ _null_ ));
+DESCR("trunc(x/y)");
+DATA(insert OID = 1980 ( numeric_div_trunc		PGNSP PGUID 12 1 0 f f t f i 2 1700 "1700 1700" _null_ _null_ _null_	numeric_div_trunc - _null_ _null_ ));
+DESCR("trunc(x/y)");
 DATA(insert OID = 2170 ( width_bucket			PGNSP PGUID 12 1 0 f f t f i 4 23 "1700 1700 1700 23" _null_ _null_ _null_	width_bucket_numeric - _null_ _null_ ));
 DESCR("bucket number of operand in equidepth histogram");
 
diff --git a/src/include/utils/builtins.h b/src/include/utils/builtins.h
index f80802f9bff..c14fcab72e8 100644
--- a/src/include/utils/builtins.h
+++ b/src/include/utils/builtins.h
@@ -7,7 +7,7 @@
  * Portions Copyright (c) 1996-2008, PostgreSQL Global Development Group
  * Portions Copyright (c) 1994, Regents of the University of California
  *
- * $PostgreSQL: pgsql/src/include/utils/builtins.h,v 1.311 2008/04/04 16:57:21 momjian Exp $
+ * $PostgreSQL: pgsql/src/include/utils/builtins.h,v 1.312 2008/04/04 18:45:36 tgl Exp $
  *
  *-------------------------------------------------------------------------
  */
@@ -845,6 +845,7 @@ extern Datum numeric_add(PG_FUNCTION_ARGS);
 extern Datum numeric_sub(PG_FUNCTION_ARGS);
 extern Datum numeric_mul(PG_FUNCTION_ARGS);
 extern Datum numeric_div(PG_FUNCTION_ARGS);
+extern Datum numeric_div_trunc(PG_FUNCTION_ARGS);
 extern Datum numeric_mod(PG_FUNCTION_ARGS);
 extern Datum numeric_inc(PG_FUNCTION_ARGS);
 extern Datum numeric_smaller(PG_FUNCTION_ARGS);
diff --git a/src/test/regress/expected/numeric.out b/src/test/regress/expected/numeric.out
index 7840990bf7d..1c5047455e8 100644
--- a/src/test/regress/expected/numeric.out
+++ b/src/test/regress/expected/numeric.out
@@ -1260,3 +1260,84 @@ SELECT * FROM num_input_test;
  -555.50
 (5 rows)
 
+--
+-- Test some corner cases for division
+--
+select 999999999999999999999::numeric/1000000000000000000000;
+        ?column?        
+------------------------
+ 1.00000000000000000000
+(1 row)
+
+select div(999999999999999999999::numeric,1000000000000000000000);
+ div 
+-----
+   0
+(1 row)
+
+select mod(999999999999999999999::numeric,1000000000000000000000);
+          mod          
+-----------------------
+ 999999999999999999999
+(1 row)
+
+select div(-9999999999999999999999::numeric,1000000000000000000000);
+ div 
+-----
+  -9
+(1 row)
+
+select mod(-9999999999999999999999::numeric,1000000000000000000000);
+          mod           
+------------------------
+ -999999999999999999999
+(1 row)
+
+select div(-9999999999999999999999::numeric,1000000000000000000000)*1000000000000000000000 + mod(-9999999999999999999999::numeric,1000000000000000000000);
+        ?column?         
+-------------------------
+ -9999999999999999999999
+(1 row)
+
+select mod (70.0,70) ;
+ mod 
+-----
+ 0.0
+(1 row)
+
+select div (70.0,70) ;
+ div 
+-----
+   1
+(1 row)
+
+select 70.0 / 70 ;
+        ?column?        
+------------------------
+ 1.00000000000000000000
+(1 row)
+
+select 12345678901234567890 % 123;
+ ?column? 
+----------
+       78
+(1 row)
+
+select 12345678901234567890 / 123;
+      ?column?      
+--------------------
+ 100371373180768845
+(1 row)
+
+select div(12345678901234567890, 123);
+        div         
+--------------------
+ 100371373180768844
+(1 row)
+
+select div(12345678901234567890, 123) * 123 + 12345678901234567890 % 123;
+       ?column?       
+----------------------
+ 12345678901234567890
+(1 row)
+
diff --git a/src/test/regress/sql/numeric.sql b/src/test/regress/sql/numeric.sql
index dc1452f9ef0..9fd6bba31ee 100644
--- a/src/test/regress/sql/numeric.sql
+++ b/src/test/regress/sql/numeric.sql
@@ -805,3 +805,21 @@ INSERT INTO num_input_test(n1) VALUES ('');
 INSERT INTO num_input_test(n1) VALUES (' N aN ');
 
 SELECT * FROM num_input_test;
+
+--
+-- Test some corner cases for division
+--
+
+select 999999999999999999999::numeric/1000000000000000000000;
+select div(999999999999999999999::numeric,1000000000000000000000);
+select mod(999999999999999999999::numeric,1000000000000000000000);
+select div(-9999999999999999999999::numeric,1000000000000000000000);
+select mod(-9999999999999999999999::numeric,1000000000000000000000);
+select div(-9999999999999999999999::numeric,1000000000000000000000)*1000000000000000000000 + mod(-9999999999999999999999::numeric,1000000000000000000000);
+select mod (70.0,70) ;
+select div (70.0,70) ;
+select 70.0 / 70 ;
+select 12345678901234567890 % 123;
+select 12345678901234567890 / 123;
+select div(12345678901234567890, 123);
+select div(12345678901234567890, 123) * 123 + 12345678901234567890 % 123;