bugGNU Octave - Bugs: bug #32924, lcm does not warn when losing...

 
 

bug #32924: lcm does not warn when losing precision

Submitted by:  Cornelius <cornim>
Submitted on:  Sun 27 Mar 2011 11:53:43 AM UTC  
 
Category: LibrariesSeverity: 1 - Wish
Priority: 3 - LowItem Group: Feature Request
Status: In ProgressAssigned to: None
Originator Name: CorneliusOpen/Closed: Open
Release: devOperating System: GNU/Linux

Add a New Comment(Rich Markup)
   

You are not logged in

Please log in, so followups can be emailed to you.

 

(Jump to the original submission Jump to the original submission)

Sat 02 Jul 2016 06:11:46 AM UTC, comment #23:

Retagging as "In progress". This seems not ready for 4.2.0.

Lachlan Andrew <lachlan>
Project Member
Thu 12 Nov 2015 07:07:35 AM UTC, comment #22:

Is this still a valid bug report?

Avinoam Kalma <avinoam>
Project Member
Tue 18 Mar 2014 10:42:40 AM UTC, comment #21:

I have added bug #41894 for MOD.

Juan Pablo Carbajal <juanpi>
Project Member
Fri 28 Feb 2014 04:01:23 PM UTC, comment #20:

I think the changes to 'mod' represent another bug report.

Rik <rik5>
Project Administrator
Fri 28 Feb 2014 09:21:04 AM UTC, comment #19:

In the mailing list we have a report about mod. The matlab implementation indeed corrects the result when roundoff generates an answer different than expected. Example follows

Gamma = 1.62e7;
duration = 10/Gamma;
dt = 0.0025/Gamma;
t = 0:dt:duration;
y = mod (t, 0.2/Gamma);
find (y==0,3,'first')

octave
1 241 401

Matlab r2008b
1 81 161

Reading the help of mod in matlab it says
MOD(x,y) is x - n.*y where n = floor(x./y) if y ~= 0. If y is not an
integer and the quotient x./y is within roundoff error of an integer,
then n is that integer.

So indeed matlab is giving a result considering roundoff error, I assume they do something like

function m = mod_ml(x,y)
if fix(y) != y
err = abs (x./y - round(x./y)) < sqrt (eps);
m = mod (x,y);
m(err) = 0;
endif
endfunction

We could issue a warning as well. Shall I open a new bug(matlab incompatibility) for mod?

Juan Pablo Carbajal <juanpi>
Project Member
Mon 21 Oct 2013 04:14:45 PM UTC, comment #18:

Ping! This patch was almost accepted 2 years ago.

Jordi, you were going to generalize it a bit. Do you think that is possible for the 3.8 release or should we wait for the first bugfix release?

Rik <rik5>
Project Administrator
Sun 02 Oct 2011 12:53:10 AM UTC, comment #17:

Thanks for the nudge, Rik. I applied Cornelius's patch, but then I realised it wasn't exactly right. Turns out it's a bit more work to do this correctly in order to account for all possible Octave types (floats, integrals, complexes, sparse). I spend some time during the past two days trying to get this to work with templates, but I couldn't do it, and I think I will just macro it instead. Sadface.

Please bear with me. I'll try to get the right patch in during the next couple of days. I'm trying to make it general enough so that other integer functions can warn when floats are used outside the integer range. After working on this patch, I realised how important it was to make this work in general.

Jordi Gutiérrez Hermoso <jordigh>
Project Administrator
Fri 30 Sep 2011 05:42:37 AM UTC, comment #16:

Jordi, Do you want to apply the patch from Cornelius and close this report? It seems it was only waiting on an e-mail address from Cornelius for crediting his work.

Rik <rik5>
Project Administrator
Mon 01 Aug 2011 07:47:29 PM UTC, comment #15:

Ha, makes sense.
I sent you a pm with my email.

Cornelius <cornim>
Mon 01 Aug 2011 07:43:04 PM UTC, comment #14:

I mean, do you have a name and an email? If not, I guess I'll say it was anonymous.

Jordi Gutiérrez Hermoso <jordigh>
Project Administrator
Mon 01 Aug 2011 06:39:32 PM UTC, comment #13:

No worries.
I have no idea how someone is usually credited for a minor patch like this. Just do what ever is common.
And thanks for taking the time.

Cornelius <cornim>
Wed 27 Jul 2011 03:23:12 AM UTC, comment #12:

Sorry for the long delay. I think I may apply your patch now. How would you like to be credited?

Jordi Gutiérrez Hermoso <jordigh>
Project Administrator
Tue 05 Apr 2011 02:23:30 PM UTC, comment #11:

Marking as patch submitted so that I don't forget later to apply your patch. It looks good at a first glance, but I need to do some basic testing to make sure.

Jordi Gutiérrez Hermoso <jordigh>
Project Administrator
Sat 02 Apr 2011 11:36:43 AM UTC, comment #10:

Hi

I didn't get a chance to look at the octave source till now but now I did and I wrote a little patch that you will hopefully be able to use.
Unfortunately because of the way octave does libraries I was not able to integrate the patch in such a way that it was testable and since I spend more than 2 hours now trying to figure out how where my source will be complied to I'm giving up now, in the hope that you or somebody else will be able to do the same work in a matter of minutes.
What I did figure out so far is that my code (which is in lo-mappers.cc) gets complied into liboctave.la (which is unfortunately not a library I can link against). I could not find where liboctave.so gets assembled and I also wasn't able to figure out way liboctave.la is not loaded when gcd.oct is loaded. Therefore I get an undefined symbol error in that case.
Well like I said, hopefully somebody with more experience will be able to fix this in a matter of minutes.
Have a look at my patch pls and let me know what you think.

(file #23077)

Cornelius <cornim>
Tue 29 Mar 2011 01:40:22 AM UTC, comment #9:

I don't see a way to remove comments in Savannah, only to flag them as spam. Sorry.

I'm going to be a little stubborn about this. But I'll reopen this bug and make it wishlist for a warning, plus retitling this bug. I don't have any particular inclination myself to write a warning about this. It seems a bit like adding a warning for the other examples I gave. Floats have limited precision. That's all there is to it. I could give examples of things being too large and losing precision instead of being too small and also losing precision. Note that there are many other locations in Octave where doubles are used for what should be integers (bincoeff and some permutation functions come to mind), so if you really want to warn for all of these, perhaps you should look into functions like xisinteger in liboctave/lo-mappers.h and put the warnings in there, or create a new function that checks if a double is too large to be an integer and warn accordingly.

If you write a patch to add a warning to src/DLD-FUNCTIONS/gcd.cc, I'll apply it. Make sure your check doesn't slow down the function considerably. Also note that the types of arguments to gcd and lcm have changed in the dev sources and vector arguments have are now interpreted for elementwise gcd and lcm, and you need at least two arguments. So you would have to rewrite your test case as lcm(25,26,27, ... , 52), or patch the functions again to recover that kind of passing the arguments.

This really seems like unnecessary effort to me. Octave is a not a number theoretical environment. These number theoretical functions should be kept simple and used sparingly, and if you need heavy number theory lifting, use Pari or Sage. If you decide to proceed with this effort, however, I will not dismiss it and apply your patch.

Jordi Gutiérrez Hermoso <jordigh>
Project Administrator
Mon 28 Mar 2011 06:31:10 PM UTC, comment #8:

The examples you bring up are
a) either wrong by about the numerical inaccuracy of the data representation or
b) give as a result a value which is not a number even though it might not be the correct one. (Would depend on the limit you are taking.)
Whereas in this case the output looks like a perfectly valid result even though it is of by orders of magnitude.
And this is neither du to incorrect input values nor giving integers which are larger than a double as input values (lcm([26:52])).
The issues arises inside the function during the for loop and therefore it is the responsibility of the function to warn the user.

FYI: In the exact same case Matlab issues the following warning:

Inputs contain values larger than the largest consecutive flint. Result may be inaccurate.

Admittedly the warning is issued inside the gcd function which should also be the case in octave considering how the lcm function is implemented.

P.S.: Please delete my previous 3 comments since they got badly mangled.

Cornelius <cornim>
Mon 28 Mar 2011 06:28:17 PM UTC, comment #7:

don't now why my answers get mangled.
but how do they say. Third time is the charm:
The examples you bring up are
a) either wrong by about the numerical inaccuracy of the data representation or
b) give as a result a value which is not a number even though it might not be the correct one. (Would depend on the limit you are taking.)
Whereas in this case the output looks like a perfectly valid result even though it is of by orders of magnitude.
The issues arises inside the function during the for loop and therefore it is the responsibility of the function to warn the user.

FYI: In the exact same case Matlab issues the following warning:
Admittedly the warning is issued inside the gcd function which should also be the case in octave considering how the lcm function is implemented.

Cornelius <cornim>
Mon 28 Mar 2011 06:27:04 PM UTC, comment #6:

The examples you bring up are
a) either wrong by about the numerical inaccuracy of the data representation or
b) give as a result a value which is not a number even though it might not be the correct one. (Would depend on the limit you are taking.)
Whereas in this case the output looks like a perfectly valid result even though it is of by orders of magnitude.

Cornelius <cornim>
Mon 28 Mar 2011 06:19:29 PM UTC, comment #5:

The examples you bring up are
a) either wrong in by about the numerical inaccuracy of the data representation or
b) give as a result a value which is not a number even though it might not be the correct one. (Would depend on the limit you are taking.)

Cornelius <cornim>
Mon 28 Mar 2011 02:06:26 PM UTC, comment #4:

Octave is a numerical computing environment. Hardware floats are the native type, and the user is thus expected to know about the limitations of hardware floats. The reason why it's using hardware floats is for speed, and the reason why many operations are unchecked is also for speed.

Octave has many other functions that silently return incorrect results if the user inputs incorrect values. For example, try tan(pi/2), sin(pi), or 0.3 - 0.2 - 0.1. All of these silently give incorrect results due to the limitations of hardware floats. Or try gamma(-1), which really should give NaN, not inf, as the limit around that simple pole doesn't exist.

Like I said before, lcm is meant as a simple convenience, not a powerful number theory function, because Octave is not a symbolic computation nor a number theoretical environment. Building in checks into lcm would slow the function down, because the principal aim here is speed, not correctness. It is up to the user to not give the function incorrect values and not use integers larger than what can be represented with a double.

Jordi Gutiérrez Hermoso <jordigh>
Project Administrator
Mon 28 Mar 2011 08:14:15 AM UTC, comment #3:

Yes, that seems to be the reason. I do have to disagree with the decision to mark this bug as invalid though. I think that no function in octave should ever return a value which is wrong/undefined. Admittedly it is a bit of an issue here but I looked at lcm.m for a bit and since calculating the lcm is an iterative process my suggestion would be that the lcm function should issue at the very least a warning once the intermediate result is so large that all significant bits do not fit into the mantissa anymore. I personally would even prefer that at this point the function should terminate with no result and an explanation for the termination. Any reason why this is not possible/unreasonable?

Cornelius <cornim>
Mon 28 Mar 2011 12:40:03 AM UTC, comment #2:

To add a bit more to Jordi's analysis, Octave is using the double class to represent the integers for the lcm analysis. The largest integer that can be represented is found with the 'bitmax' function. For standard IEEE 64-bit doubles the result is 9.0072e+15. The fact that you are getting results of order 1e20 means that you have gone well beyond the ability of the double class to represent the numbers involved.

Even using the 'uint64' class will only get you numbers up to ~1.8e19. If you need something more you will need to look into a package with arbitrary precision.

Rik <rik5>
Project Administrator
Sun 27 Mar 2011 04:05:39 PM UTC, comment #1:

The lcm function is meant as a simple convenience, but it doesn't replace a number theory package with arbitrary precision like Pari or Sage. The errors you are seeing are due to the size of your inputs and that the data type is a finite-precision double float, not an arbitrary precision integer. Octave doesn't have this data type natively, although the symbolic package in Octave-Forge defines the vpa function.

The problem with the order of the arguments is a simple fact of numerical analysis with finite precision; the canonical example usually is computing exp(x) from its power series expansion starting with the largest terms first vs starting with the smallest terms first.

I'm therefore marking this bug as invalid. I hope you agree.

Jordi Gutiérrez Hermoso <jordigh>
Project Administrator
Sun 27 Mar 2011 11:53:43 AM UTC, original submission:

Quite obviously the result of lcm function should be independent of the order of the input and should upon adding a number to the array always yield a result which is >= the previous result.
There I would expect the following

but in reality

Cornelius <cornim>

 

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

Attach File(s):
   
   
Comment:
   

Attached Files
file #23077:  diff.txt added by cornim (2KiB - text/plain)

 

Depends on the following items: None found

Digest:
   bug dependencies.

 

Carbon-Copy List
  • -unavailable- added by lachlan (Updated the item)
  • -unavailable- added by avinoam (Posted a comment)
  • -unavailable- added by juanpi (Posted a comment)
  • -unavailable- added by rik5 (Posted a comment)
  • -unavailable- added by jordigh (Posted a comment)
  • -unavailable- added by cornim (Submitted the item)
  •  

    Do you think this task is very important?
    If so, you can click here to add your encouragement to it.
    This task has 0 encouragements so far.

    Only project members can vote.

     

    Please enter the title of George Orwell's famous dystopian book (it's a date):

     

     

    Follow 15 latest changes.

    Date Changed By Updated Field Previous Value => Replaced By
    Sat 02 Jul 2016 06:11:46 AM UTClachlanStatusPatch Submitted=>In Progress
    Mon 08 Feb 2016 03:58:49 AM UTClachlanAssigned tojordigh=>None
    Thu 01 Mar 2012 04:32:14 PM UTCjordighDependencies-=>bugs #35680 is dependent
    Fri 15 Apr 2011 02:09:36 AM UTCjordighAssigned toNone=>jordigh
    Tue 05 Apr 2011 02:23:29 PM UTCjordighStatusConfirmed=>Patch Submitted
    Sat 02 Apr 2011 11:36:43 AM UTCcornimAttached File-=>Added diff.txt, #23077
    Tue 29 Mar 2011 01:40:22 AM UTCjordighSeverity3 - Normal=>1 - Wish
      Priority5 - Normal=>3 - Low
      Item GroupIncorrect Result=>Feature Request
      StatusInvalid=>Confirmed
      Open/ClosedClosed=>Open
      Release3.2.4=>dev
      Summarylcm yields wrong results=>lcm does not warn when losing precision
    Sun 27 Mar 2011 04:05:39 PM UTCjordighStatusNone=>Invalid
      Open/ClosedOpen=>Closed

    Back to the top


    Powered by Savane 3.1-cleanup1