Boost :

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 )

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