
Boost : 
From: Andras Erdei (aerdei_at_[hidden])
Date: 20041222 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_i2 / b_i1
p_i := p_i1 * a_i + p_i2
q_i := q_i1 * a_i + q_i2
if p_i or q_i overflows, return (p_i1,q_i1)
b_i := b_i2 % b_i1
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 MatulaKornerup: Foundations
of Finite Precision Rational Arithmetic, Computing Suppl. 2, 1980, pp 85111)
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)
 medianrounding: 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^311 (a 32 bit int) probability of rouding errors:
P( error > 1e18 ) < 0.43
P( error > 1e15 ) < 4.3e4
P( error > 1e11 ) < 4.3e8
P( error > 1e9.3 ) = 0
br,
andras
Boost list run by bdawes at acm.org, gregod at cs.rpi.edu, cpdaniel at pacbell.net, john at johnmaddock.co.uk