Re: [boost] rational with finite precision

From: Andras Erdei (aerdei_at_[hidden])
Date: 2004-12-22 04:56:54

On Mon, 20 Dec 2004 19:52:30 -0500, Daryle Walker <darylew_at_[hidden]> wrote:

>I would like to see precise directions of the algorithms that give the
>closest representable result (that have to try avoiding internal overflow).

i'm not sure what you mean by "have to try avoiding internal overflow":
the algorithm is the Euclidean algorithm for computing the GCD, and it
computes the best rational approximation when there is an overflow, but
it does not avoid the overflow

euclid(n,m)

 b_-2 = n ; p_-2 = 0 ; q_-2 = 1
 b_-1 = m ; p_-1 = 1 ; q_-1 = 0

 loop

   a_i := b_i-2 / b_i-1

   p_i := p_i-1 * a_i + p_i-2
   q_i := q_i-1 * a_i + q_i-2

   if p_i or q_i overflows, return (p_i-1,q_i-1)

   b_i := b_i-2 % b_i-1

   if 0 == b_i, return (p_i,q_i)

with our previous example (191/181 + 197/193 = 72520/34933)

 i  a   b      p      q      p/q

-2      72520      0      1
-1      34933      1      0

 0   2   2654      2      1  2.00...
 1  13    431     27     13  2.076...
 2   6     68    164     79  2.07594...
 3   6     23   1011    487  2.075975...
 4   2     22   2186   1053  2.075973...
 5   1      1   3197   1540  2.07597402...
 6  22      0  72520  34933  2.075974007... overflow, return (3197,1540)

* * * *

the math background (no proofs, simplified, from Matula-Kornerup: Foundations
of Finite Precision Rational Arithmetic, Computing Suppl. 2, 1980, pp 85-111)

the convergents p_i/q_i and algorithm have the properties:

- best rational approximation:

        if r/s != p/q, s <= q_i then abs( r/s - p/q ) > abs( p_i/q_i - p/q )

- quadratic convergence:

        1/( q_i * ( q_i+1 + q_i ) ) < abs( p_i/q_i = p/q ) <= 1/( q_i * q_i+1 )

- alternating convergence:

        p_0/q_0 < p_2/q_2 < ... <= p/q <= ... < p_3/q_3 < p_1/q_1

- if a convergent cannot be represented (overflow) then no
  subsequent convergent can be represented

- fixed point: iff p/q can be represented, then euclid(p/q)==(p,q)

- monotonic: if x < y then euclid(x) <= euclid(y)

- median-rounding: if a value falls between two representable
  numbers, the median of those is the splitting point for the
  rounding; the bounds are rounded towards the simpler rational

- with rational<n> the rouding error at most 1/n, but usually it's
  around 1/n^2; the wide rounding intervals belong to simpler
  fractions (where q is small compared to n)

with n = 2^31-1 (a 32 bit int) probability of rouding errors:

  P( error > 1e-18 ) < 0.43
  P( error > 1e-15 ) < 4.3e-4
  P( error > 1e-11 ) < 4.3e-8
  P( error > 1e-9.3 ) = 0

br,
andras