Re: Arbitrary precision modulo operation - Mailing list pgsql-general

> -----Original Message-----
> From: Tom Lane [mailto:tgl@sss.pgh.pa.us]
> Sent: Wednesday, April 28, 2004 9:14 PM
> To: Dann Corbit
> Cc: Bruno Wolff III; Paul Tillotson; pgsql-general@postgresql.org
> Subject: Re: [GENERAL] Arbitrary precision modulo operation
>
>
> "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?

Here is what I would suggest:

Using the outline I proposed before (starting with a floating point
divide of DBL_DIG digits of precision), keep doubling the precision
until the precision is 5 digits larger than either operand.  If the last
doubling makes the precision larger (quite likely) simply reduce it to
the smaller margin.

Something like this (pseudocode):

numeric          divide(numeric x, numeric y)
{
    /* double starting estimate of quotient*/
    numeric y1 = 1.0 / y::double;
    numeric two = 2:numeric;
    /* Need to collect the maximal precision of either operand */
    total_precision = get_max_precision(x,y);

    y1 *= (two - y * y1); /* use numeric 30 math or total precision+5
whichever is less */
    if (total_precision <= 25) return x*y; /* correctly rounded -->
banker's rounding? */
    y1 *= (two - y * y1); /* use numeric 60 math or total precision+5
whichever is less */
    if (total_precision <= 50) return x*y; /* correctly rounded -->
banker's rounding? */
    y1 *= (two - y * y1); /* use numeric 120 math or total precision+5
whichever is less */
    if (total_precision <= 110) return x*y; /* correctly rounded -->
banker's rounding? */
    y1 *= (two - y * y1); /* use numeric 240 math or total precision+5
whichever is less */
    if (total_precision <= 230) return x*y; /* correctly rounded -->
banker's rounding? */
    y1 *= (two - y * y1); /* use numeric 480 math or total precision+5
whichever is less */
    if (total_precision <= 470) return x*y; /* correctly rounded -->
banker's rounding? */
    y1 *= (two - y * y1); /* use numeric 960 math or total precision+5
whichever is less */
    if (total_precision <= 950) return x*y; /* correctly rounded -->
banker's rounding? */
    y1 *= (two - y * y1); /* use maximum precision math or total
precision+5 whichever is less */
    return x*y;
}