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

 
 

bug #32924: lcm does not warn when losing precision

Submitter:  Cornelius <cornim>
Submitted:  Sun 27 Mar 2011 11:53:43 AM UTC
   
 
Category:  Libraries Severity:  1 - Wish
Priority:  3 - Low Item Group:  Feature Request
Status:  Fixed Assigned to:  None
Originator Name:  Cornelius Open/Closed:  * Closed
Release:  * dev Operating System:  * GNU/Linux
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Sat 17 Sep 2022 08:26:12 AM UTC, comment #28: 

Thanks, Rik. Pushed here: https://hg.savannah.gnu.org/hgweb/octave/rev/6646f2b5a3d1

Closing as fixed.

Arun Giridhar <arungiridhar>
Group Member
Sat 17 Sep 2022 01:11:39 AM UTC, comment #27: 

The patch looks good to me.  I would change the last "end" to "endif" and put a newline after it to separate it from "endfunction".

Rik <rik5>
Group administrator
Thu 15 Sep 2022 09:00:48 PM UTC, comment #26: 

Added some BISTs. Updated diff / patch attached. Passes BISTs.


(file #53691)

Arun Giridhar <arungiridhar>
Group Member
Thu 15 Sep 2022 08:12:43 PM UTC, comment #25: 

Going through old bug reports and this looks like a simple change analogous to what we already do in nchoosek.m, so I copied and adapted these warnings from nchoosek:


diff -r 0dec459a4064 scripts/specfun/lcm.m
--- a/scripts/specfun/lcm.m     Wed Sep 14 09:59:31 2022 -0400
+++ b/scripts/specfun/lcm.m     Thu Sep 15 16:00:16 2022 -0400
@@ -53,6 +53,13 @@ function l = lcm (varargin)
     l(msk) = 0;
   endfor

+  if (isfloat (l) && l > flintmax (l))
+    warning ("Octave:lcm:large-output-float", ...
+             "lcm: possible loss of precision");
+  elseif (isinteger (l) && l == intmax (l))
+    warning ("Octave:lcm:large-output-integer", ...
+             "lcm: result may have saturated at intmax");
+  end
 endfunction


For floating point inputs:

octave:15> lcm (num2cell(1:40){:})
ans = 5342931457063200

octave:16> lcm (num2cell(1:41){:})
warning: lcm: possible loss of precision
warning: called from
    lcm at line 57 column 5

ans = 2.190601897395912e+17


For integer inputs:

octave:27> lcm (num2cell(uint64(1:46)){:})
ans = 9419588158802421600

octave:28> lcm (num2cell(uint64(1:47)){:})
warning: lcm: result may have saturated at intmax
warning: called from
    lcm at line 60 column 5

ans = 18446744073709551615


If there are no objection to this change in the next few days, I'll push it.





@nrjank: If you're trying to recreate in Matlab, this sequence should work:

l = 26; for i = 27:52, l = lcm(l, i), end             % forward, double

l = 26; for i = 52:-1:27, l = lcm(l, i), end          % reverse, double

l = uint64(26); for i = 27:52, l = lcm(l, i), end     % forward, integer

l = uint64(26); for i = 52:-1:27, l = lcm(l, i), end  % reverse, integer


Arun Giridhar <arungiridhar>
Group Member
Sun 24 Mar 2019 03:47:56 AM UTC, comment #24: 

adding that (a) this still appears to be valid in 5.1.0, as no warnings are generated by lcm

(b) lcm appears to no longer(?) accept input as given in the example. I couldn't find an example of matlab lcm accepting input that way either, except as a shadowed lcm in the symbolic toolbox.

Matlab 2018b:

>> lcm([26:52])
Not enough input arguments.

>> lcm(num2cell(26:52){:})
Error: ()-indexing must appear last in an index expression.

>> a = num2cell(26:52);
>> lcm(a{:})
Error using lcm
Too many input arguments.


Octave 5.1.0

>> lcm([26:52])
error: Invalid call to lcm.  Correct usage is:
 -- lcm (X, Y)
 -- lcm (X, Y, ...)

>> lcm(num2cell([26:52]){:})
ans =    2.3971e+25

>> lcm(num2cell([52:-1:26]){:})
ans =    1.1213e+30

>> lcm(num2cell(25:52){:})
ans =    1.1866e+28


so, it's still here, but I can't recreate it in matlab to know if that exact situation still produces an error/warning.

Nicholas Jankowski <nrjank>
Group Member
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>
Thu 12 Nov 2015 07:07:35 AM UTC, comment #22: 


Is this still a valid bug report?

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

I have added bug #41894 for MOD.

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

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

Rik <rik5>
Group 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>
Group 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>
Group 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>
Group Member
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>
Group 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>
Group Member
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>
Group Member
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>
Group Member
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>
Group Member
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.

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.

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.
And this is neither du to incorrect input values nor giving integers which are larger than a double as input values (+verbatim+ lcm([26:52]) -verbatim-).
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:
+verbatim+ Inputs contain values larger than the largest consecutive flint. Result may be inaccurate. -verbatim-
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: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.)

Whereas in my example the resulting value seems like a perfectly fine and possibly correct result even though it is by orders of magnitude of. This error is neither due to incorrect input values nor is any of the input values larger than what can be represented as a double (+verbatim+ lcm([26:52]) -verbatim-).
The issue only arises in the for loop inside the function and therefore it is the functions responsiblity to warn the user.

FYI: In the exact same case Matlab issues the following warning:
+verbatim+ Inputs contain values larger than the largest consecutive flint. Result may be inaccurate. -verbatim-
Admittedly the warning stems from the gcd function which it should in octave too considering the way the lcm function is implemented.

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>
Group Member
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>
Group 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>
Group Member
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


lcm([26:52]) == lcm([52:-1:26]) <= lcm([25:52])


but in reality


>>> lcm([26:52])
ans =  8.8544e+19
>>> lcm([52:-1:26])
ans =  5.1651e+20
>>> lcm([25:52])
ans =  4.0955e+17


Cornelius <cornim>

 

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

Attach Files:
   
   
Comment:
   

Attached Files
file #53691:  lcm.patch added by arungiridhar (1KiB - text/x-patch)
file #23077:  diff.txt added by cornim (2KiB - text/plain)

 

Depends on the following items: None found

Digest:
   bug dependencies.

 

Carbon-Copy List
  • -email is unavailable- added by arungiridhar (Posted a comment)
  • -email is unavailable- added by nrjank (Posted a comment)
  • -email is unavailable- added by lachlan (Updated the item)
  • -email is unavailable- added by avinoam (Posted a comment)
  • -email is unavailable- added by juanpi (Posted a comment)
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by jordigh (Posted a comment)
  • -email is unavailable- added by cornim (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 20 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2022-09-17 arungiridhar StatusPatch Reviewed Fixed
        Open/ClosedOpen Closed
    2022-09-17 rik5 StatusPatch Submitted Patch Reviewed
    2022-09-15 arungiridhar Attached File- Added lcm.patch, #53691
    2022-09-15 arungiridhar StatusIn Progress Patch Submitted
    2016-07-02 lachlan StatusPatch Submitted In Progress
    2016-02-08 lachlan Assigned tojordigh None
    2012-03-01 jordigh Dependencies- bugs #35680 is dependent
    2011-04-15 jordigh Assigned toNone jordigh
    2011-04-05 jordigh StatusConfirmed Patch Submitted
    2011-04-02 cornim Attached File- Added diff.txt, #23077
    2011-03-29 jordigh Severity3 - 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
    2011-03-27 jordigh StatusNone Invalid
        Open/ClosedOpen Closed

    Back to the top

    Powered by Savane 3.13-cf05.
    Corresponding source code