On Mon, 20 Dec 2004 19:52:30 -0500, Daryle Walker <darylew@hotmail.com> 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