bugGNU Octave - Bugs: bug #46859, floating point error in range...

 
 

bug #46859: floating point error in range results in non-integer value

Submitter:  None
Submitted:  Sat 09 Jan 2016 02:21:21 PM UTC
   
 
Category:  Interpreter Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Incorrect Result
Status:  Fixed Assigned to:  None
Originator Name:  Francesco Menoncin Originator Email:  -email is unavailable-
Open/Closed:  * Closed Release:  * dev
Operating System:  * Any Fixed Release:  None
Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Fri 17 Aug 2018 12:40:04 AM UTC, comment #12: 

I checked in a changeset that forces the internal rng_limit to be an integer if the base and increment are integers.  This is executed just once per construction of a range so there is no performance hit.  I also added BIST tests to check the behavior.  See https://hg.savannah.gnu.org/hgweb/octave/rev/dc47c9e48801.

Marking bug as fixed and closing report.

Rik <rik5>
Group administrator
Thu 16 Aug 2018 08:31:23 PM UTC, comment #11: 

The remaining issue concerns how to handle the endpoint of a Range when it is not exactly an integer, but is within eps of an integer.  In Matlab, and more generally in most users expectations, it appears that the endpoint should be rounded to the integer value.  Marking the bug as Confirmed and changing the category to Matlab Compatibility.

Rik <rik5>
Group administrator
Thu 16 Aug 2018 08:28:51 PM UTC, comment #10: 

I think this issue is important and should be re-visited.  The error message that results from indexing with a double value that is close to an integer obscures what is going on.  Today's behavior:


octave:1> x = [1:5];
octave:2> x(1+100*eps)
error: x(1): subscripts must be either integers 1 to (2^63)-1 or logicals


The patch from Lachlan to print out more detail about the invalid index no longer applies.  I updated it to target the correct file and to use C++ rather than C library functions.  I checked that in here: https://hg.savannah.gnu.org/hgweb/octave/rev/0360ed7c39a8.

The new behavior is


octave:2> x(1+100*eps)
error: x(1+2.22045e-14): subscripts must be either integers 1 to (2^63)-1 or logicals



Rik <rik5>
Group administrator
Thu 14 Jan 2016 04:00:24 AM UTC, comment #9: 

Oops.  The previous patch printed the fractional part too often (|| instead of &&).  Version 2 should work.

(file #36022)

Lachlan Andrew <lachlan>
Mon 11 Jan 2016 10:07:12 AM UTC, comment #8: 

Lachlan's patch for the error message seems good
to me.

Also, the "and only return values up to 1000 as Matlab does,"
is likely a good choice. But, there may be other ideas??

Michael

Michael Godfrey <godfrey>
Group Member
Mon 11 Jan 2016 06:52:11 AM UTC, comment #7: 

Here is a patch to show small fractional parts of indices in the error message, as Michael suggested.  It shows


>> a = 1:(1001/250)/(1/250);
>> b(a(end)) = 1
error: b(1001-1.13687e-13): subscripts must be either integers 1 to (2^31)-1 or logicals


It looks a little unwieldy, but I think it is clearer than hex, and certainly less misleading than omitting the fractional part.

However, the main problem remains: given  sum(ones(1001,1))==1001, why does the range contain a non-integer final value?  I assume the problem is in


Range::max (void) const
{
  double retval = 0.0;
  if (rng_numel > 0)
    {
      if (rng_inc > 0)
        {
          retval = rng_base + (rng_numel - 1) * rng_inc;

          // On some machines (x86 with extended precision floating point
          // arithmetic, for example) it is possible that we can overshoot the
          // limit by approximately the machine precision even though we were
          // very careful in our calculation of the number of elements.
          // Therefore, we clip the result to the limit if it overshoots.
          // The test also includes equality (>= rng_limit) to have expressions
          // such as -5:1:-0 result in a -0 endpoint.
          if (retval >= rng_limit)
            retval = rng_limit;


We could check for integrality of the base and step and not return a fractional rng_limit, but that would slow things down.

My vote would be to treat the upper limit literally, and only return values up to 1000 as Matlab does, even if that confuses some people...

(file #35984)

Lachlan Andrew <lachlan>
Mon 11 Jan 2016 05:03:45 AM UTC, comment #6: 

Michale, I like your idea of showing that the answer isn't exactly an integer.

Many of the comments below are about the fact that the end-point is not well defined if floating point is involved.  However, the real puzzle is that adding integers much less than the mantissa size gives a non-integer result.

Representing 1000 only takes 10 bits, whereas a double should be able to do integer arithmetic for numbers up to 52 bits without error, shouldn't it?

Lachlan Andrew <lachlan>
Sat 09 Jan 2016 10:04:24 PM UTC, comment #5: 

I did a little playing around in Matlab and it seems that they consider float values to be equal if they are within +/- eps(x)/2.

E.g.
 
 % Matlab
 x = 500
 x == (x+eps(x)/2)
 x == (x+eps(x)/1.99999999)

 500
 1
 0

Likewise,

 numel(1:(x-eps(x)/2))
 numel(1:(x-eps(x)/1.99999999))

 500
 499

Doing the same thing in Octave gives

 % Octave
 500
 1
 0

and

 500
 500

So it looks like Octave's range function is using different logic to determine whether to include the endpoint. It seems to allow a wider tolerance up to about 6-7 times eps(x):

 numel(1:(x-eps(x)*6))
 500

 numel(1:(x-eps(x)*7))
 499

Ceral Paquet <octavebugs>
Sat 09 Jan 2016 05:29:01 PM UTC, comment #4: 

John: Right, finding the "best" solution is not easy.
Matlab compatibility is a consideration, but does not
dictate the right way to do it.

A part of this that may confuse users is:


>>zz =T/dt
zz =  1001.0

>>xx = 1001
xx = 1001

but, of course:
format hex
>> zz = T/dt
zz = 408f47ffffffffff

>> xx
xx = 408f480000000000


Maybe enhancing the message:
error: x(1000): subscripts must be either integers 1 to (2^31)-1 or logicals

by displaying the failing subscript in hex or at least somehow
indicating explicitly that it does not have an exact integer
value.

Michael Godfrey <godfrey>
Group Member
Sat 09 Jan 2016 04:35:38 PM UTC, comment #3: 

I agree that there is a bug of some sort, and certainly an incompatibility with Matlab.  But I don't know what the best solution is.

In any case, my experience has been that when people come across things I like this they typically expect something different.  the usual expectation is that the result of T/dt should be exactly 1001 and that the range 2:T/dt should have 1000 elements, not 999.  And that's what Octave tries to do here, but then there is still a possible bug because it sets the end point of the range to the result of the floating point calculation T/dt instead of 1001 even though the beginning of the range and the increment are both integers, and (again, if we were doing exact arithmetic) the endpoint would be as well.

These sorts of things have given us all kinds of trouble, and trying to be "smart" about them by going to great lengths to compute the "proper" number of elements for ranges can NEVER give the results that please everyone all the time.

John W. Eaton <jwe>
Group administrator
Sat 09 Jan 2016 03:43:37 PM UTC, comment #2: 

Matlab 2015a says:


>> T=1001/250;
>> dt=1/250;
>> x(1)=10;
>> for i=2:T/dt
x(i)=x(i-1)+randn();
end
>> whos
  Name      Size              Bytes  Class     Attributes

  T         1x1                   8  double
  dt        1x1                   8  double
  i         1x1                   8  double
  x         1x1000             8000  double

>>


So, it is a bug.

Michael Godfrey <godfrey>
Group Member
Sat 09 Jan 2016 03:42:23 PM UTC, comment #1: 

Also Matlab 2013a runs this script without complaint. The final values of the variables are:


>> T


T =

    4.0040

>> dt


dt =

    0.0040

>> i


i =

        1000

>>

Ceral Paquet <octavebugs>
Sat 09 Jan 2016 02:21:21 PM UTC, original submission:  

T=1001/250;
dt=1/250;
x(1)=10;
for i=2:T/dt
x(i)=x(i-1)+randn();
end

These commands give the following error message (since T/dt is not exactly equal to 1001):

error: subscript indices must be either positive integers less than 2^31 or logicals

When T=1000/250, for instance, the error is not produced and the result is correct.

I have checked with other software (Scilab and R) and this mistake is not produced.

Anonymous

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

Attached Files
file #36022:  bug_46859_gripe.v2.cset added by lachlan (1KiB - application/octet-stream)
file #35984:  bug_46859_gripe.cset added by lachlan (1KiB - application/octet-stream)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by jordigh (Updated the item)
  • -email is unavailable- added by mtmiller (Updated the item)
  • -email is unavailable- added by lachlan (Posted a comment)
  • -email is unavailable- added by jwe (Posted a comment)
  • -email is unavailable- added by godfrey (Posted a comment)
  • -email is unavailable- added by octavebugs (Posted a comment)
  • -email is unavailable- added by None (Submitted the item)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 13 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2018-08-17 rik5 Item GroupMatlab Compatibility Incorrect Result
        StatusConfirmed Fixed
        Open/ClosedOpen Closed
    2018-08-16 rik5 Item GroupUnexpected Error or Warning Matlab Compatibility
        StatusPatch Submitted Confirmed
        Release4.0.0 dev
    2018-08-16 rik5 Item GroupIncorrect Result Unexpected Error or Warning
    2016-06-18 jordigh SummaryError in &quot;for&quot; cycle floating point error in range results in non-integer value
    2016-02-21 mtmiller CategoryNone Interpreter
        StatusNone Patch Submitted
        Operating SystemGNU/Linux Any
    2016-01-14 lachlan Attached File- Added bug_46859_gripe.v2.cset, #36022
    2016-01-11 lachlan Attached File- Added bug_46859_gripe.cset, #35984

    Back to the top

    Powered by Savane 3.13-4448.
    Corresponding source code