Re: Arbitrary precision modulo operation

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