Re: loop problem with exit condition

["Dmitry A. Kazakov" <mailbox@dmitry-kazakov.de>]

On Tue, 26 Oct 2010 10:02:32 +0300, Niklas Holsti wrote:

> Yannick Duchêne (Hibou57) wrote:
>> Le Sun, 26 Sep 2010 11:10:24 +0200, Dmitry A. Kazakov 
>> <mailbox@dmitry-kazakov.de> a écrit:
>>> in your case the simplest way to remove bias would be to multiply the 
>>> step:
>>> […]
>>>          Inner := Inner_Start + Inner_Step * Long_Float (J - 1);
>> 
>> I did not knew that. Why is an iterative addition less accurate than a 
>> multiplication ? Does it have something to deal with adjustment of the 
>> exponent ? (the one of the float's internal representation)
> 
> Yes. If you add two floating-point numbers A+B, and the binary exponent 
> Ae of A is larger than the exponent Be of B, the mantissa of B will be 
> shifted right by Ae-Be bits before being added to the mantissa of A, so 
> the addition loses ("shifts out") about this many of the less 
> significant bits of the mantissa of B. You are effectively adding A+Bt, 
> where Bt is B with the lost bits set to zero (assuming positive numbers 
> -- I haven't figured out what happens with negative numbers).
> 
> If you compute A+N*B by iteratively adding B to A, N times, you are 
> effectively computing the approximation A+N*Bt (more or less; if the 
> binary exponent of the growing sum becomes larger than Ae, you will lose 
> even more precision). The error is N*(B-Bt), thus growing with N.
> 
> In the formula A+N*B, as suggested by Dmitry (and as indeed is good 
> practice), the multiplication N*B is as accurate as the floating-point 
> representation allows -- no mantissa bits are lost by exponent 
> adjustments. If the product is X, the addition A+X is again as accurate 
> as possible. Some bits of the smaller term (A or X) may again be lost in 
> the addition, but this loss is not cumulative.

For more complicated cases there exist summation algorithms reducing
precision loss, e.g. the Kahan algorithm:

   http://en.wikipedia.org/wiki/Kahan_summation_algorithm

(Numeric methods are exiting)

-- 
Regards,
Dmitry A. Kazakov
http://www.dmitry-kazakov.de