Re: a-b+b != a

[stn021]

Well, then it is easy to solve. Simply round by some power of 2:

mul = 2^20 ;
r1 = round( mul*randn(5) ) / mul ;
r2 = round( mul*randn(5) ) / mul ;
r3 = r1 - r2 + r2 ;
not_zero = r1 - r3

not_zero =

   0   0   0   0   0
   0   0   0   0   0
   0   0   0   0   0
   0   0   0   0   0
   0   0   0   0   0



2017-09-04 17:58 GMT+02:00 James Sherman Jr. <address@hidden>:
> On Mon, Sep 4, 2017 at 10:53 AM, stn021 <address@hidden> wrote:
>> Hi,
>>
>> this is weird:
>>
>> r1 = round( 1e3*randn(5) ) / 1e3 ;
>> r2 = round( 1e3*randn(5) ) / 1e3 ;
>> r3 = r1 - r2 + r2 ;
>> not_zero = r1 - r3
>>
>> not_zero =
>>
>>    0.0000e+00   0.0000e+00   5.5511e-17   0.0000e+00   1.1102e-16
>>    5.5511e-17   0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00
>>    0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00  -2.0817e-16
>>   -5.5511e-17   0.0000e+00   0.0000e+00  -1.1102e-16   0.0000e+00
>>    2.7756e-17   0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00
>>
>> So a-b+b != a
>>
>> The difference to zero is small, around 1e-17. But iterations can
>> cause this error to increase.
>>
>> I use leasqr() and in each iteration the last line is
>>   retval = retval - someval + someval
>>
>> With that additional line I get quite different results compared to
>> the same program without this line even though they should be
>> identical.
>>
>> Is there a way to avoid this phenomenon ?
>>
>> THX
>> Stefan
>>
>> _______________________________________________
>> Help-octave mailing list
>> address@hidden
>> https://lists.gnu.org/mailman/listinfo/help-octave
>
> This is fundamental issue caused by the use of floating point numbers.
> I'd suggest reading the article to get a better background in floating
> point, but the key section is
> https://en.wikipedia.org/wiki/Floating-point_arithmetic#Accuracy_problems
>
> One way to avoid this issue to never test equality (or not equality)
> when dealing with floating point numbers.  So instead of
>>if ( a == b )
>>  do stuff
>>end
> instead use something like:
>>if ( abs(a-b) < eps(a) )
>> do stuff
>>end
>
> Hope this helps,
> James Sherman