Re: A problem with range objects and floating point numbers

[s.jo]

Daniel Sebald wrote
> On 09/23/2014 01:12 AM, Oliver Heimlich wrote:
>> Am 23.09.2014 03:27, schrieb s.jo:
>>> John W. Eaton wrote
>>>> On 06/20/2014 04:28 PM, Julien Bect wrote:
>>>>> Hello all,
>>>>>
>>>>> I would like to bring the following bug report to your attention :
>>>>>
>>>>> http://savannah.gnu.org/bugs/?42589
>>>>>
>>>>> A simple example of the regression described in this report is the
>>>>> following :
>>>>>
>>>>>       octave:1>  (9:10) * .1 - 1
>>>>>       ans =
>>>>>
>>>>>          -1.0000e-01    2.7756e-17
>>>>>
>>>>> The second element should be 0.
>>>>>
>>>>> Rik thinks that this bug might be hard to solve... which is why I'm
>>>>> bringing this up on this list :)
>>>>
>>>> Please see my comments on the bug tracker.
>>>>
>>>> jwe
>>>
>>> I have the similar problem with octave 3.8.1 on cygwin.
>>> See below:
>>> octave:50>  t=[-2:0.1:0]
>>> t =
>>>
>>>   Columns 1 through 5:
>>>
>>>                      -2                  -1.9                  -1.8
>>> -1.7                  -1.6
>>>
>>>   Columns 6 through 10:
>>>
>>>                    -1.5                  -1.4                  -1.3
>>> -1.2                  -1.1
>>>
>>>   Columns 11 through 15:
>>>
>>>                      -1                  -0.9                  -0.8
>>> -0.7                  -0.6
>>>
>>>   Columns 16 through 20:
>>>
>>>                    -0.5                  -0.4                  -0.3
>>> -0.2   -0.0999999999999999
>>>
>>>   Column 21:
>>>
>>>                       0
>>>
>>> Range command is not good with floating-point step only starting from
>>> negative integer as above.
>>>
>>> I am trying to upgrade octave 3.8.2 to see if this problem is fixed.
>>>
>>> I post it just in case.
>>
>> The results are perfectly okay. This is how floating point works, see
>> IEEE 754. Above mentioned bug report also gives some explanation in the
>> comments.
>>
>> The “problem” occurs because .1 is not a binary floating point number.
>> and the approximation
>> 0.1000000000000000055511151231257827021181583404541015625 is used
>> instead. 10 * 0.1 suffers (or in this case benefits) from rounding.
>>
>> In the last example with t=[-2:0.1:0] only the numbers -2, -1.5, -1,
>> -0.5, and 0 have the exact value that is displayed. Column 20
>> additionally yields a propagation of errors. You can check this with the
>> bit format:
>>
>> octave:3>  format bit
>> octave:4>  -.1
>> ans = 1011111110111001100110011001100110011001100110011001100110011010
>> (this is -0.1000000000000000055511151231257827021181583404541015625)
>>
>> octave:5>  -2 + 19 * .1
>> ans = 1011111110111001100110011001100110011001100110011001100110010000
>> (this is -0.09999999999999997779553950749686919152736663818359375)
>>
>> octave:6>  t(20)
>> ans = 1011111110111001100110011001100110011001100110011001100110010000
> 
> Thank you for the explanation Oliver.
> 
> Would there be any utility to attempting to "integerize" the range, if 
> possible, and thereby eliminate the accumulation of errors?  For 
> example, [-2:0.1:0] is equivalent to [-20:1:0]*.1, so if the limits turn 
> out to be factorable, then internally the range could be represented 
> slightly different than what the user types.
> 
>    scale = 1.0;
> ...
>    if (range.limits_are_factorable ())
>      {
>        start /= step;
>        scale = step;
>      }
> ...
>    ridx.xelem (i) = (start + i*step)*scale;
> 
> To be factorable doesn't necessarily require that the floating point 
> step size be an exact number system equivalent, i.e., that .1 equal 
> exactly 1/10.  It just requires the machine arithmetic to be as expected 
> when producing an operation result, i.e.,
> 
> octave-cli:1> rem(-2,.1) == 0
> ans =  1
> octave-cli:2> rem(0,.1) == 0
> ans =  1
> octave-cli:3> [-20:1:0]*.1
> ans =
> 
>   Columns 1 through 5:
> 
>    -2.00000  -1.90000  -1.80000  -1.70000  -1.60000
> 
>   Columns 6 through 10:
> 
>    -1.50000  -1.40000  -1.30000  -1.20000  -1.10000
> 
>   Columns 11 through 15:
> 
>    -1.00000  -0.90000  -0.80000  -0.70000  -0.60000
> 
>   Columns 16 through 20:
> 
>    -0.50000  -0.40000  -0.30000  -0.20000  -0.10000
> 
>   Column 21:
> 
>     0.00000
> 
> 
> Jo, what do you get for the above scripts?
> 
> Dan

Dan and others,
Sorry for late reply.

Your integerized stepping seems to be a decimation process.
You used step-size 1. Any bigger integers such as 10 or 10^17 are allowed in
the decimation process.

Comparing the results of colon operators with floating-point numbers,
Matlab seems to use such decimation process, but Octave does not.
I agree that this decimation feature of colon operator is far from any IEEE
standard.

Also I point out that the results of linspace function and colon operator in
Matlab are different.
This may imply decimation processes are different. 
By inspection, the linspace result is much close to decimal number.
I guess that the decimation process of Matlab's colon operator is simpler
for speed.

I rephase my question: 
Does Octave's colon operator with floating-point numbers need to use
decimation as Matlab does?

At the early conversation, I think that the compatibility with Matlab is
most important for usual Octave users. Now I understand that this may
mislead other unknown features.

So far, I learned that integer range can be generated 'safely' with colon
operator.
Regarding to floating-point number range, we may want to use the function
"linspace" for the best decimation.

I don't know who is responsible for the colon operator source code.
But I think that if there is a document for comparing colon operator with
linspace function,
that will be helpful.

There is another question: I found that (-20:1:1)*.1 and [-20:1:1]*0.1 have
different results.


The range with [...] is better decimated than the case of (...).
More interestingly, if I transpose the range array  such as (-20:1:1)'*.1,
then I have the same result in [-20:1:1]'*.1.
Is this a bug? Did someone say that there is a patch for this?
I am useing Ocatve 3.8.2 on cygwin.

Thank you guys for many comments.

Best regards,

jo

ps. 
I re-opened this topic for my curiosity. 
It is fun to learn how Octave works and to compare it with Matlab.
You has better email me to see any more comparison results. 
Sharing raw data seems on this conversation looks like spam.





--
View this message in context: 
http://octave.1599824.n4.nabble.com/A-problem-with-range-objects-and-floating-point-numbers-tp4664900p4666710.html
Sent from the Octave - Maintainers mailing list archive at Nabble.com.