Re: Is there a bug in the inversion of floats

[JuanPi]

On Wed, Apr 27, 2016 at 6:26 PM, Oliver Heimlich <address@hidden> wrote:
> Am 27. April 2016 14:12:23 MESZ, schrieb JuanPi <address@hidden>:
>>I get the following results and I am not sure this is a bug or there
>>is something I do not understand in the way the floats are
>>represented/treated
>>
>>format native-hex
>>
>>v = complex (0,1)
>>inv(exp(v)),exp(-v), 1./exp(v)
>>ans = 8c06b50f284ae13f  ee0c098f54edeabfi
>>ans = 8c06b50f284ae13f  ee0c098f54edeabfi
>>ans = 8c06b50f284ae13f  ee0c098f54edeabfi
>>
>>v = complex (1,0)
>>inv(exp(v)),exp(-v), 1./exp(v)
>>ans = 38ef2c36568bd73f
>>ans = 38ef2c36568bd73f
>>ans = 38ef2c36568bd73f
>>
>>so far so good, but
>>
>>inv(exp(v)),exp(-v), 1./exp(v)
>>ans = 21d8befb2a71c93f  545f8539d4cfd3bfi
>>ans = 23d8befb2a71c93f  555f8539d4cfd3bfi
>>ans = 23d8befb2a71c93f  555f8539d4cfd3bfi
>>
>>(inv gives different results) and thought: ok but who uses inv to
>>invert numbers...
>>so I did
>>
>>inv(exp(v/99)),exp(-v/99), 1./exp(v/99)
>>ans = f30d005d41adef3f  52ef0272867a84bfi
>>ans = f30d005d41adef3f  50ef0272867a84bfi
>>ans = f30d005d41adef3f  52ef0272867a84bfi
>>
>>(exp(-a) gives different results) and started worrying
>>
>>Essentially I hoped to get a definite false here, but I did not
>>
>>t = linspace (0, 1, 100);
>>v = complex(0,1);
>>V = v * ( t - t.');
>>M = exp (-V) != 1 ./ exp (V);
>>any(M(:))
>>
>>In this case the result is also true for any v
>>
>>Should one desire that the result of exp(-v) is exactly the same as
>>1/exp(v) for consistency?
>
> Regarding your last question,  exp (-v) computes the value without 
> intermediate rounding errors whereas 1 ./ exp (v) has to invert an 
> intermediate result that has rounding errors and inversion isn't error free.
>
> So I would expect exp (-v) to computer a better value in general.
>
> Oliver
I guess it depends on how the interpreter parses 1 ./ exp (v), but if
indeed you do the the two operations they propagate more error.
The underlying question is what does matlab returns here...

t = linspace (0, 1, 100);
V = bsxfun (@minus,  t , t.');
M = exp (-V) != 1 ./ exp (V);
any(M(:))


-- 
JuanPi Carbajal
Public GnuPG key: 9C5B72BF
-----
The end of funding: "Many researchers were caught up in a web of
increasing exaggeration."
- Hans Moravec