Arbitrary precision modulo operation

On Mon, Apr 26, 2004 at 10:18:52 -0400,
Chadwick Boggs <chadwickboggs@yahoo.com> wrote:

I need to perform modulo operations on extremely large numbers. The %
operator is giving me number out of range errors and the mod(x, y)
function simply seems to return the wrong results. Also, my numerator
is in the format of a quoted string, which the mod function can't take.

How large is extremely large?
You can cast the strings to a numeric type to solve the string problem.
'numeric' should work for numbers up to about 1000 digits.
For example:
area=> select '1234567890'::numeric % '123'::numeric;
?column?
----------
39
(1 row)

If you are getting wrong results you should post a specific example so
that the developers can figure out what is going wrong.

Example of wrong results from modulo operation of arbitrary precision
numbers:

# select '123456789012345678901234567890'::numeric % 123;
?column?
----------
-6
(1 row)

# select mod('123456789012345678901234567890'::numeric, 123);
mod
-----
-6
(1 row)

The correct result (at least according to another, unnamed, RDBMS):

select '123456789012345678901234567890' % 123;

+----------------------------------------+
| '123456789012345678901234567890' % 123 |
+----------------------------------------+
|                                     58 |
+----------------------------------------+
1 row in set (0.00 sec)

Bruno Wolff III wrote:

On Mon, Apr 26, 2004 at 10:18:52 -0400,
Chadwick Boggs <chadwickboggs@yahoo.com> wrote:

I need to perform modulo operations on extremely large numbers. The %
operator is giving me number out of range errors and the mod(x, y)
function simply seems to return the wrong results. Also, my numerator
is in the format of a quoted string, which the mod function can't take.

I tried some large values and seem to be getting reasonable results.
I did get some negative remainders, but I didn't find anything in the
mod documentation on whether or not these were allowed. For small values
I always got positive remainders so it isn't consistantly return the
remainder closest to zero. My guess is that it has to do rounding when
dividing numerics.

On Mon, Apr 26, 2004 at 13:30:17 -0400,
Chadwick Boggs <chadwickboggs@yahoo.com> wrote:

Example of wrong results from modulo operation of arbitrary precision
numbers:

# select '123456789012345678901234567890'::numeric % 123;
?column?
----------
-6
(1 row)

# select mod('123456789012345678901234567890'::numeric, 123);
mod
-----
-6
(1 row)

The correct result (at least according to another, unnamed, RDBMS):

select '123456789012345678901234567890' % 123;

+----------------------------------------+
| '123456789012345678901234567890' % 123 |
+----------------------------------------+
|                                     58 |
+----------------------------------------+
1 row in set (0.00 sec)

I checked this with bc and I got -6 (and 117) as being correct. I would
think the other database was wrong.

It wouldn't happen to be MYSQL would it?

Bruno, perhaps round is an issue. Thank you. Here is an example that
should involve no rounding and indeed it works:

Multiply the ten largest integer scale prime numbers:

# select 2147483477::numeric * 2147483489::numeric * 2147483497::numeric
* 2147483543::numeric * 2147483549::numeric * 2147483563::numeric *
2147483579::numeric * 2147483587::numeric * 2147483629::numeric *
2147483647::numeric;
?column?
------------------------------------------------------------------------------------------------
2085923946138988916149190605561960475118165298582929035878182900998428077414994652962618167119
(1 row)

Now, modulo this by any of the factors correctly returns 0:

# select
'2085923946138988916149190605561960475118165298582929035878182900998428077414994652962618167119'::numeric%
2147483563;
?column?
----------
0
(1 row)

Also, diviing out all of the factors correctly returns 1:

# select
2085923946138988916149190605561960475118165298582929035878182900998428077414994652962618167119
/ 2147483477 / 2147483489 / 2147483497 / 2147483543 / 2147483549 /
2147483563 / 2147483579 / 2147483587 / 2147483629 / 2147483647;
?column?
------------------------
1.00000000000000000000
(1 row)

This provides the solution to my problem: I merely needed to cast all
of the numbers to numberic for product and modulo operations.

Thank you,
Chadwick.

Bruno Wolff III wrote:

Maple output:
y := 123456789012345678901234567890 mod 123;
y := 117

CONNX output (which uses qfloat by S. Moshier):
select mod(123456789012345678901234567890 , 123) {nopassthrough}
117

PariGP output:
? Mod(123456789012345678901234567890, 123)
%4 = Mod(117, 123)
?

I see there are a few misconceptions about numeric and modulus on here:

(1) A modulus operation on a numeric type should NOT have rounding
errors. The whole point of numeric is that it is an arbitrary precision
BASE 10 representation of your number. The modulus returns the (whole
number) remainder as a result of a division.

(2) the modulus operator/function is, AFAIK, supposed to return the
modulus with the SAME SIGN as the divisor, so I think this is a bug.
That's what every other modulus operator that I have ever seen does.
Would you mind doing

foodb=> SELECT version();

(3) MySQL just rounds large numbers to the highest value that the type
will support, and apparently, no arbitrary precision types are listed on
this page:

http://dev.mysql.com/doc/mysql/en/Numberic_type_overview.html

You can't expect to get a modulus from an out-of-range number.

Paul Tillotson

Chadwick Boggs wrote:

-----Original Message-----
From: Paul Tillotson [mailto:pntil@shentel.net]
Sent: Monday, April 26, 2004 4:41 PM
To: pgsql-general@postgresql.org
Subject: Re: [GENERAL] Arbitrary precision modulo operation

I see there are a few misconceptions about numeric and
modulus on here:

(1) A modulus operation on a numeric type should NOT have
rounding errors. The whole point of numeric is that it is an
arbitrary precision BASE 10 representation of your number.

This is true

The modulus returns the (whole
number) remainder as a result of a division.

This is true if the numeric values are integers.

When the values are not integral, some non-integral results can be
returned.

2.50 % 2.50 is 0
But 13.89 modulo 3.50 is 3.39
If you work it out on paper, you will see that 3.39 is the correct
remainder.

(2) the modulus operator/function is, AFAIK, supposed to
return the modulus with the SAME SIGN as the divisor, so I
think this is a bug. That's what every other modulus operator
that I have ever seen does. Would you mind doing

I would agree that a positive modulus is preferable. However, the
negative result is also mathematically correct.

foodb=> SELECT version();

(3) MySQL just rounds large numbers to the highest value that
the type will support, and apparently, no arbitrary precision
types are listed on this page:

http://dev.mysql.com/doc/mysql/en/Numberic_type_overview.html

You can't expect to get a modulus from an out-of-range number.

Should it not (therefore) throw an error of some sort?

On Mon, Apr 26, 2004 at 12:48:45PM -0700, Dann Corbit wrote:

Maple output:
y := 123456789012345678901234567890 mod 123;
y := 117

PgSQL 7.3.6 gives the right answer (117), 7.4 gets it wrong (-6). Most
likely a bug was introduced when NUMERIC was rewritten. Strange it
hasn't been noticed before.

--
Alvaro Herrera (<alvherre[a]dcc.uchile.cl>)
"Vivir y dejar de vivir son soluciones imaginarias.
La existencia est� en otra parte" (Andre Breton)

In article <408D4729.303@yahoo.com>,
Chadwick Boggs <chadwickboggs@yahoo.com> writes:

Example of wrong results from modulo operation of arbitrary precision
numbers:

# select '123456789012345678901234567890'::numeric % 123;
?column?
----------
-6
(1 row)

# select mod('123456789012345678901234567890'::numeric, 123);
mod
-----
-6
(1 row)

The correct result (at least according to another, unnamed, RDBMS):

select '123456789012345678901234567890' % 123;

+----------------------------------------+
| '123456789012345678901234567890' % 123 |
+----------------------------------------+
|                                     58 |
+----------------------------------------+
1 row in set (0.00 sec)

Is the name of the other, unnamed RDBMS by chance starting with an
'M'? ;-)

Anyway, both are wrong. According to GNU MP, the correct answer is
117. Thus PostgreSQL is at least almost right (117 - 123 = -6).

Alvaro Herrera wrote:

On Mon, Apr 26, 2004 at 12:48:45PM -0700, Dann Corbit wrote:

Maple output:
y := 123456789012345678901234567890 mod 123;
y := 117

PgSQL 7.3.6 gives the right answer (117), 7.4 gets it wrong (-6). Most
likely a bug was introduced when NUMERIC was rewritten. Strange it
hasn't been noticed before.

mod(x, y) is computed as x - trunc(x / y) * y in the mod_var() function
(I think).

However, it appears that the division operator itself is rounding up,
such that the trunc() function (which ought to round down) does no good
as a round up has already occurred.

Thus, the value of (x / y) is 1 too large, and so x % y is actually
giving you (x % y) - y, a negative number. I tried looking at how the
division actually works, but it is over my head at least for the 30
minute perusal.

Regards,
Paul Tillotson

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

[paul@pjt4 paul]$ bc
bc 1.06
Copyright 1991-1994, 1997, 1998, 2000 Free Software Foundation, Inc.
This is free software with ABSOLUTELY NO WARRANTY.
For details type `warranty'.
111111111111111111 / 6
18518518518518518

[paul@pjt4 bin]$ ./psql -U postgres template1
Welcome to psql 7.4.2, the PostgreSQL interactive terminal.

Type: \copyright for distribution terms
\h for help with SQL commands
\? for help on internal slash commands
\g or terminate with semicolon to execute query
\q to quit

template1=# select 111111111111111111::numeric / 6;
?column?
-------------------
18518518518518519
(1 row)

template1=# select 111111111111111111 / 6;
?column?
-------------------
18518518518518518
(1 row)

template1=# select version();

version

---------------------------------------------------------------------------------------------------------
PostgreSQL 7.4.2 on i686-pc-linux-gnu, compiled by GCC gcc (GCC) 3.3.2
20031022 (Red Hat Linux 3.3.2-1)
(1 row)

Paul Tillotson <pntil@shentel.net> writes:

mod(x, y) is computed as x - trunc(x / y) * y in the mod_var() function

However, it appears that the division operator itself is rounding up,

This is because div_var rounds its output to the number of fractional
digits specified by select_div_scale, which saith

/*
* The result scale of a division isn't specified in any SQL standard.
* 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.
*/

You get the "correct" answer if you force a couple more digits to be
carried by increasing either input's dscale, eg

regression=# select 123456789012345678901234567890.00 % 123;
?column?
----------
117.00
(1 row)

It could be that select_div_scale needs to allow some more slop, or that
that routine is okay but mod_var shouldn't depend on it to select the
intermediate result scale for its internal division. Comments?

regards, tom lane

On Wed, Apr 28, 2004 at 02:01:00 -0400,
Tom Lane <tgl@sss.pgh.pa.us> wrote:

Paul Tillotson <pntil@shentel.net> writes:

mod(x, y) is computed as x - trunc(x / y) * y in the mod_var() function

It could be that select_div_scale needs to allow some more slop, or that
that routine is okay but mod_var shouldn't depend on it to select the
intermediate result scale for its internal division. Comments?

One option would be to define a separate division operator that always
returns an integral value and that is truncated toward 0 and use that
for the mod function. People might find this operator useful in itself.

Another option would be for mod to check if the remainder doesn't have the
same sign as the divisor and if so keeping adding the divisor to it until
it does.

-----Original Message-----
From: Bruno Wolff III [mailto:bruno@wolff.to]
Sent: Wednesday, April 28, 2004 12:05 AM
To: Tom Lane
Cc: Paul Tillotson; pgsql-general@postgresql.org
Subject: Re: [GENERAL] Arbitrary precision modulo operation

On Wed, Apr 28, 2004 at 02:01:00 -0400,
Tom Lane <tgl@sss.pgh.pa.us> wrote:

Paul Tillotson <pntil@shentel.net> writes:

mod(x, y) is computed as x - trunc(x / y) * y in the mod_var()
function

It could be that select_div_scale needs to allow some more slop, or
that that routine is okay but mod_var shouldn't depend on

it to select

the intermediate result scale for its internal division. Comments?

One option would be to define a separate division operator
that always returns an integral value and that is truncated
toward 0 and use that for the mod function. People might find
this operator useful in itself.

It will give some wrong results. The result of mod should be the
remainder after division, which is not always integral in the case of
numeric fixed point or floating point.
Consider the output of this program:

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <time.h>
int main(void)
{
int i;
double den,
other;
srand((unsigned) time(NULL));
for (i = 0; i < 10; i++) {
den = rand() + 1.0;
other = 1.0 + rand() / (rand() + 1.0);
printf("%f mod %f is %f\n", other, den, fmod(other, den));
printf("%f mod %f is %f\n", den, other, fmod(den, other));
}
return 0;
}

Another option would be for mod to check if the remainder
doesn't have the same sign as the divisor and if so keeping
adding the divisor to it until it does.

I would suggest computation in sufficient digits of accuracy to get a
correct answer. If any precision is lost in numeric calculations then
an error of PLOSS should be returned (unless somehow there is a total
loss of precision).

I do not know how division is performed in numeric quantities for
PostgreSQL.

However, if the algorithm is being revised, here is a possible outline:

1. Convert to double and perform the division to get an estimate
2. Use Newton's method for division, starting at 30 digits and doubling
the digits for each iteration.
3. Stop when you have reached 50% excess digits or so
4. Round to the correct value.

Since the computations use smaller significant digits until the last
iteration, there is typically a large savings compared to performing all
calculations in full precision.

So, if the target is 70 digits you will have:

~15 digits precision after the double precision floating point divide
30 digits after one iteration
60 digits after two iterations
120 digits after three iterations

Then round to the correct length.

See:
http://www.azillionmonkeys.com/qed/adiv.html

For a short discussion of Newton's method for division.

/*
** Hopefully, I am not annoying anyone with all of this.
** This is a brief demonstration of how Newton's division algorithm
works:
*/
#include <stdio.h>
double divide(double x, double y)
{
unsigned __int64 est;
double y1 = 1.0 / y;
double y2;
puts("\nForming reciprocal:");
est = y1 * 10000;
y1 = est / 10000.0; /* make an estimate good to 4 decimal
places... */
printf("x=%.7f, y=%.7f, y1=%.7f\n", x, y, y1);
y1 *= (2 - y * y1); /* should be 8 places */
printf("x=%.14f, y=%.14f, y1=%.14f\n", x, y, y1);
y1 *= (2 - y * y1); /* In theory ~16 digits */
printf("x=%.20f, y=%.20f, y1=%.20f\n", x, y, y1);
y1 *= (2 - y * y1); /* One more is really needed here... */
printf("x=%.20f, y=%.20f, y1=%.20f\n", x, y, y1);
puts("\ndividing:");
printf("result of division is:%f\n", x * y1);
return 0;
}

int main(void)
{
divide(1.0, 2.0);
divide(1.0, 17.0);
divide(3.0, 19.0);
divide(19.0, 3.0);
return 0;
}

On Wed, Apr 28, 2004 at 14:02:57 -0700,
Dann Corbit <DCorbit@connx.com> wrote:

-----Original Message-----
From: Bruno Wolff III [mailto:bruno@wolff.to]

One option would be to define a separate division operator
that always returns an integral value and that is truncated
toward 0 and use that for the mod function. People might find
this operator useful in itself.

It will give some wrong results. The result of mod should be the
remainder after division, which is not always integral in the case of
numeric fixed point or floating point.
Consider the output of this program:

The remainder may not be integral but the quotient should be. However the
idea of getting the quotient, multiplying by the divisor and subtracting
from the dividend is not very good from an efficiancy point of view.

On Wed, Apr 28, 2004 at 03:06:25PM -0700, Dann Corbit wrote:

/*
** Hopefully, I am not annoying anyone with all of this.
** This is a brief demonstration of how Newton's division algorithm
works:
*/

Well, have you looked at Postgres implementation? I don't know if it's
Newton or not, but it certainly is not simple (I didn't take the time to
understand it).

Have a look at src/backend/utils/adt/numeric.c function div_var()

Cc: list trimmed ...

--
Alvaro Herrera (<alvherre[a]dcc.uchile.cl>)
"�C�mo puedes confiar en algo que pagas y que no ves,
y no confiar en algo que te dan y te lo muestran?" (Germ�n Poo)

"Dann Corbit" <DCorbit@connx.com> writes:

I would suggest computation in sufficient digits of accuracy to get a
correct answer.

Fine. How many is that, exactly?

regards, tom lane