BUG #3387: mod on non-integer returns bad result

"Filip Krska" <filip.krska@comstar.cz> writes:

select mod (70.0,70) from dual;
returns
70.0
instead of
0.0

PG 8.0 gets this right. I think this demonstrates that Bruce's 8.1 patch
http://archives.postgresql.org/pgsql-committers/2005-06/msg00045.php
didn't actually fix anything, merely move the failure cases around.

regards, tom lane

Gregory Stark <stark@enterprisedb.com> writes:

The source of the problem is the floating point arithmetic which is used to do
the individual steps in the long division.

I don't think so. The ultimate source of the problem is that div_var
can only report a finite number of digits.

regards, tom lane

"Tom Lane" <tgl@sss.pgh.pa.us> writes:

"Filip Krska" <filip.krska@comstar.cz> writes:

select mod (70.0,70) from dual;
returns
70.0
instead of
0.0

PG 8.0 gets this right. I think this demonstrates that Bruce's 8.1 patch
http://archives.postgresql.org/pgsql-committers/2005-06/msg00045.php
didn't actually fix anything, merely move the failure cases around.

Well I think that patch is right in itself.

The source of the problem is the floating point arithmetic which is used to do
the individual steps in the long division. It does seem odd to me that it
isn't using integer arithmetic for that step. It says it has to do avoid
overflow though?

I would start with this though. 1/70 isn't exactly representable so 70.0 *
1.0/70 doesn't work properly whereas 70.0 / 70.0 has at least a chance of
working. But I do think switching this one way or another to integer math
would be the real solution.

Index: numeric.c
===================================================================
RCS file: /home/stark/src/REPOSITORY/pgsql/src/backend/utils/adt/numeric.c,v
retrieving revision 1.104
diff -u -r1.104 numeric.c
--- numeric.c	9 Jun 2007 15:52:30 -0000	1.104
+++ numeric.c	15 Jun 2007 16:22:11 -0000
@@ -4052,7 +4052,6 @@
 	NumericDigit *res_digits;
 	double		fdividend,
 				fdivisor,
-				fdivisorinverse,
 				fquotient;
 	int			qi;
 	int			i;
@@ -4128,7 +4127,6 @@
 		if (i < var2ndigits)
 			fdivisor += (double) var2digits[i];
 	}
-	fdivisorinverse = 1.0 / fdivisor;

 	/*
 	 * maxdiv tracks the maximum possible absolute value of any div[] entry;
@@ -4152,7 +4150,7 @@
 				fdividend += (double) div[qi + i];
 		}
 		/* Compute the (approximate) quotient digit */
-		fquotient = fdividend * fdivisorinverse;
+		fquotient = fdividend / fdivisor;
 		qdigit = (fquotient >= 0.0) ? ((int) fquotient) :
 			(((int) fquotient) - 1);	/* truncate towards -infinity */

@@ -4203,7 +4201,7 @@
 						fdividend += (double) div[qi + i];
 				}
 				/* Compute the (approximate) quotient digit */
-				fquotient = fdividend * fdivisorinverse;
+				fquotient = fdividend / fdivisor;
 				qdigit = (fquotient >= 0.0) ? ((int) fquotient) :
 					(((int) fquotient) - 1);	/* truncate towards -infinity */
 				maxdiv += Abs(qdigit);
@@ -4236,7 +4234,7 @@
 	fdividend = (double) div[qi];
 	for (i = 1; i < 4; i++)
 		fdividend *= NBASE;
-	fquotient = fdividend * fdivisorinverse;
+	fquotient = fdividend / fdivisor;
 	qdigit = (fquotient >= 0.0) ? ((int) fquotient) :
 		(((int) fquotient) - 1);	/* truncate towards -infinity */
 	div[qi] = qdigit;

--
Gregory Stark
EnterpriseDB http://www.enterprisedb.com

"Tom Lane" <tgl@sss.pgh.pa.us> writes:

Gregory Stark <stark@enterprisedb.com> writes:

The source of the problem is the floating point arithmetic which is used to do
the individual steps in the long division.

I don't think so. The ultimate source of the problem is that div_var
can only report a finite number of digits.

In the case reported div_var was getting 70/70 = 0.99999. Which is really just
wrong. The only reason was because 1.0/70 isn't representable so (1.0/70) * 70
is slightly more than 1 which ... div_var really ought not have any trouble
representing an integer even in its finite number of digits.

--
Gregory Stark
EnterpriseDB http://www.enterprisedb.com

Gregory Stark <stark@enterprisedb.com> writes:

In the case reported div_var was getting 70/70 = 0.99999. Which is
really just wrong.

Agreed, but I think your proposed patch is just a band-aid. The real
problem with div_var is that it generates inaccurate output digits at
all --- it's assuming that computing a few guard digits before rounding
will suffice to ensure that all the delivered digits are correct, but
I think there will always be corner cases where it fails.

I just blew the dust off my old copy of Knuth vol 2, and see that his
algorithm for multi-precision division generates output digits that are
correct to start with (or at least he never needs to revisit a digit
after moving on to the next). ISTM we should go over to an approach
like that. For the truncated-output case the result will be exact,
and for the rounded-output case you generate exactly one guard digit
to see if it's >= base/2.

regards, tom lane

This has been saved for the 8.4 release:

http://momjian.postgresql.org/cgi-bin/pgpatches_hold

---------------------------------------------------------------------------

Tom Lane wrote:

Gregory Stark <stark@enterprisedb.com> writes:

In the case reported div_var was getting 70/70 = 0.99999. Which is
really just wrong.

Agreed, but I think your proposed patch is just a band-aid. The real
problem with div_var is that it generates inaccurate output digits at
all --- it's assuming that computing a few guard digits before rounding
will suffice to ensure that all the delivered digits are correct, but
I think there will always be corner cases where it fails.

I just blew the dust off my old copy of Knuth vol 2, and see that his
algorithm for multi-precision division generates output digits that are
correct to start with (or at least he never needs to revisit a digit
after moving on to the next). ISTM we should go over to an approach
like that. For the truncated-output case the result will be exact,
and for the rounded-output case you generate exactly one guard digit
to see if it's >= base/2.

regards, tom lane

--
Bruce Momjian <bruce@momjian.us> http://momjian.us
EnterpriseDB http://www.enterprisedb.com

+ If your life is a hard drive, Christ can be your backup. +