Re: Arbitrary precision modulo operation

["Dann Corbit" <DCorbit@connx.com>]

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

---------------------------(end of broadcast)---------------------------
TIP 4: Don't 'kill -9' the postmaster