bugGNU Octave - Bugs: bug #65495, nchoosek error "gcd: all...

 
 

bug #65495: nchoosek error "gcd: all values must be integers"

Submitter:  None
Submitted:  Wed 20 Mar 2024 03:32:14 PM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Regression
Status:  Ready For Test Assigned to:  None
Originator Name:  Originator Email:  -email is unavailable-
Open/Closed:  * Open Release:  * dev
Operating System:  * Any Fixed Release:  None
Planned Release:  10.1.0 (current default)
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Sun 21 Apr 2024 11:53:37 AM UTC, comment #18: 

I've added a commit message and a self-test (to avoid another regression for this example), adapted to Octave coding style, and slightly changed the diff from comment #13 to potentially reduce the number of calculations. (I hope I didn't mess up at the latter part.)
I pushed that change to the default branch with Hendrik as the author:
https://hg.savannah.gnu.org/hgweb/octave/rev/60772ed3e920

Thank you again for the contribution. And sorry for taking so long with the review.

Marking as ready for test.

Markus Mützel <mmuetzel>
Group administrator
Sat 13 Apr 2024 11:53:59 PM UTC, comment #17: 

Can someone review the proposed fix in comment #13?

Rik <rik5>
Group administrator
Fri 22 Mar 2024 10:27:33 AM UTC, comment #16: 

In Octave (and Matlab), large literal integer values are preserved exactly as integers if they appear as the argument to the int64 or uint64 functions.  In Octave, this behavior is implemented by storing a special integer value type that is either converted to a double value when it is first stored in a variable OR converted to an int64 or uint64 value by the int64/uint64 functions.

See bug #45945 and related patches.  This behavior probably does deserve some mention in the documentation.  Unfortunately, I don't see anything about it in a NEWS file either.  Oops, that's probably my fault.

See also https://hg.savannah.gnu.org/hgweb/octave/rev/96e7dc4c2214 for a change related to specifying the type of literal integer values.  Maybe there was one, but I haven't located an associated bug report.  That change in syntax and behavior was mentioned in the NEWS file for version 7, but may not be described in the manual.  It would probably be worth adding at least the info from the NEWS file to the section in the manual that covers literal values.

John W. Eaton <jwe>
Group administrator
Fri 22 Mar 2024 05:14:36 AM UTC, comment #15: 

Correct means of course analytically correct.

A numerical floating point solution is beyond flintmax not always analytically correct but is an approximate solution. Specifically in our case beyond flintmax one cannot assume, that i/(g1 * g2) is identical to one and therefore one needs to cater for these cases in the programming logic.


The analytical correct solution of the gcd between 879236959347161415 and 9 is actually 9.
That 879236959347161415 cannot exactly be represented as a double and that its "representation" is 879236959347161472 and that the gcd of its "representation" and 9 is 3 is true, but somewhat besides the point made below: one cannot assume that i/(g1 * g2) is equal to 1. I did not mean that the gcd does provide a numerically incorrect result beyond flintmax...


@Side note:
I am not aware of any major programming language which treats numerical literals/constants as doubles when they are used in the program to e.g. set integer variable values.  In fact, I would consider such a behaviour a major limitation/error.
So I am not sure that one needs to document normal/expected behaviour.


 


Hendrik K <koerhen>
Fri 22 Mar 2024 04:02:54 AM UTC, comment #14: 


> b)
> The other issue is that the gcd function does NOT always provide a correct result for floating point numbers beyond flintmax.
> Example:
> gcd (uint64(879236959347161415), 9)
> -> 9
> gcd (double(879236959347161415), 9)
> -> 3


`double(879236959347161415)` is the same as `879236959347161415` which has an actual value of 879236959347161472.  The GCD of this and 9 is 3.

Side note:  I thought numeric literals were always doubles, which would make `uint64(879236959347161415)` have an actual value of 879236959347161472, but I just checked in Octave 8.4.0 and it seems to have an actual value of 879236959347161415.  I don't recall seeing mention of this in the documentation.

Anonymous
Thu 21 Mar 2024 05:28:22 AM UTC, comment #13: 

Interesting case:

a)
One issue is that C may become Inf within the loop, so the gcd function (gcd(Inf)) complains and throws the observed error. This is easy to check and then break the loop within nchoosek accordingly.
Note: The gcd error message could be refined when receiving Inf compared to receiving e.g. a fraction, though.

b)
The other issue is that the gcd function does NOT always provide a correct result for floating point numbers beyond flintmax.
Example:
gcd (uint64(879236959347161415), 9)
-> 9
gcd (double(879236959347161415), 9)
-> 3

This impacts nchoosek: In theory, the final division by "i/(g1 * g2)" is NOT necessary, as it should be one by design. Unfortunately this is not always the case for floating point numbers as shown above (integer numbers have never these issues).

Once C times (v - k + i) reaches i_max and before the the (next) C iteration (may) reach Inf, we continue to work with the gcd for numerical precision reasons in each iteration. Or as the gcd function value for floats could be too small, we may temporarily reach Inf before dividing by "i / (g1 * g2)" which would not be identical to one in this case.

Below is the proposed changed code


    for i = 1:k
      if (C * (v - k + i) >= imax)
        ## Avoid overflow / precision loss by determining the smallest
        ## possible factor of (C * (n-k+i)) and i via the gcd.
        ## Note that by design in each iteration
        ##   1) C will always increase (factor is always > 1).
        ##   2) C will always be a whole number.
        ## Therefore, using the gcd will always provide the best possible
        ## solution until saturation / has the least precision loss.
        g1 = gcd (C, i);
        g2 = gcd (v - k + i, i/g1);
        C /= g1;

        ## In theory and (always for integers) i/(g1 * g2) is identical to 1 by
        ## design. Or for floats and beyond flintmax, the gcd may not be
        ## correctly derived by the gcd function and i/(g1 * g2) may not be 1.
        if (i/(g1 * g2) == 1) || !isinf(C * (v - k + i)/g2)
          C *= (v - k + i)/g2;
          C /= i/(g1 * g2);
        else
          C /= i/(g1 * g2);
          ## We have potential precision loss by dividing (too) early, but
          ## advantage is that we prevent possible interim overflows
          C *= (v - k + i)/g2;
        endif
        if (is_int && (C == imax)) || (! is_int && isinf(C))
          break;  # Stop here; saturation reached.
        endif
      else
        C *= (v - k + i);
        C /= i;
      endif
    endfor

 

@Markus-
you were right: When changing the order of multiplication and division, we are expected to get a precision loss in certain cases, so we only do this when it is absolutely necessary in the proposed code change.

 " || !isinf(C * (v - k + i)/g2))




(file #55873)

Hendrik K <koerhen>
Wed 20 Mar 2024 05:35:47 PM UTC, comment #12: 

Re-ordering the multiplications and division like this, leads to the same result for `nchoosek (1024, 512)` like in Octave 9:

diff -r a79e07f980ff scripts/specfun/nchoosek.m
--- a/scripts/specfun/nchoosek.m        Fri Mar 15 07:24:58 2024 -0400
+++ b/scripts/specfun/nchoosek.m        Wed Mar 20 18:34:14 2024 +0100
@@ -140,11 +140,11 @@
         g1 = gcd (C, i);
         g2 = gcd (v - k + i, i/g1);
         C /= g1;
+        C /= i/(g1 * g2);
         C *= (v - k + i)/g2;
         if (is_int && (C == imax))
           break;  # Stop here; saturation reached.
         endif
-        C /= i/(g1 * g2);
       else
         C *= (v - k + i);
         C /= i;


But that might un-do the precision improvements of Hendrik's change. (I haven't checked though.)

Markus Mützel <mmuetzel>
Group administrator
Wed 20 Mar 2024 05:20:03 PM UTC, comment #11: 

I wonder if we should interrupt the for loop when C reaches Inf.
Maybe something along the lines of the following?

diff -r a79e07f980ff scripts/specfun/nchoosek.m
--- a/scripts/specfun/nchoosek.m        Fri Mar 15 07:24:58 2024 -0400
+++ b/scripts/specfun/nchoosek.m        Wed Mar 20 18:19:07 2024 +0100
@@ -149,6 +149,9 @@
         C *= (v - k + i);
         C /= i;
       endif
+      if (isinf (C))
+        break;
+      endif
     endfor
     if (! is_int && C > imax)
       warning ("Octave:nchoosek:large-output-float", ...


CC'ing Hendrik in case they have a better idea.

Markus Mützel <mmuetzel>
Group administrator
Wed 20 Mar 2024 04:35:59 PM UTC, comment #10: 

I think the bug report is still valid (at least the error message should be improved), but the attribution is wrong.
"Release: " should be "dev" or 10.0.0.
Also "Operating System" should be "all"

Dmitri.
--

Dmitri A. Sergatskov <dasergatskov>
Wed 20 Mar 2024 04:31:36 PM UTC, comment #9: 

This is on linux, octave 10 (current dev snapshot)


octave:1> version
ans = 10.0.0
octave:2> version -hgid
ans = d5c5484322d1
octave:3> nchoosek (5000, 2000)
error: gcd: all values must be integers
error: called from
    nchoosek at line 140 column 12
octave:4> nchoosek (1024, 512)
error: gcd: all values must be integers
error: called from
    nchoosek at line 140 column 12


Dmitri.
--

Dmitri A. Sergatskov <dasergatskov>
Wed 20 Mar 2024 04:31:32 PM UTC, comment #8: 

Octave is installed by system admin on SecureCRT. They only said it was updated to 9.1.0 so I thought that is error but it was on wrong machine, not cloud machine but test machine.

Anonymous
Wed 20 Mar 2024 04:29:42 PM UTC, comment #7: 

Sorry sorry, was connected to wrong computer. Was supposed to be on HP cloud but was on app test station. Now version is 9.1.0 and nchoosek is correct, no gcd error.

Anonymous
Wed 20 Mar 2024 04:27:42 PM UTC, comment #6: 

You are running octave 10.0.0 (development version).
How did you install octave?

Dmitri.
--

Dmitri A. Sergatskov <dasergatskov>
Wed 20 Mar 2024 04:17:22 PM UTC, comment #5: 

Cannot copy paste from Octave window (SecureCRT) so typing. Exact error is "error: gcd: all values must be integers. Error: called from nchoosek line 140 column 12".

How to tell BLAS version? Typing "version" gives 10.0.0.

Anonymous
Wed 20 Mar 2024 04:07:13 PM UTC, comment #4: 

A search within nchoosek, gcd, and the full source code doesn't find a message "GCD should be integer" or variations around that.  Can you please copy and paste the exact output from octave around the error?

Nicholas Jankowski <nrjank>
Group Member
Wed 20 Mar 2024 04:00:50 PM UTC, comment #3: 

similarly for your second example:


>> version
ans = 8.4.0
>> nchoosek (1024, 512)
warning: nchoosek: possible loss of precision
warning: called from
    nchoosek at line 151 column 7

ans = 4.4813e+306



>> version
ans = 9.1.0
>> nchoosek (1024, 512)
warning: nchoosek: possible loss of precision
warning: called from
    nchoosek at line 151 column 7

ans = 4.4813e+306


I don't know if it's relevant but do you know what BLAS you are using? my results above are the same whether i use reference BLAS or OpenBLAS.


>> version -blas
ans = unknown or reference BLAS



Nicholas Jankowski <nrjank>
Group Member
Wed 20 Mar 2024 03:54:46 PM UTC, comment #2: 

I'm not seeing this.  Also on windows:


octave:1> version
ans = 8.4.0
octave:2> nchoosek (5000, 2000)
warning: nchoosek: possible loss of precision
warning: called from
    nchoosek at line 151 column 7

ans = Inf



ans = 9.1.0
>> nchoosek (5000, 2000)
warning: nchoosek: possible loss of precision
warning: called from
    nchoosek at line 151 column 7

ans = Inf


Nicholas Jankowski <nrjank>
Group Member
Wed 20 Mar 2024 03:40:35 PM UTC, comment #1: 

nchoosek (1024, 512) used to give 4.4813e306 but now gives error.

Anonymous
Wed 20 Mar 2024 03:32:14 PM UTC, original submission:  

nchoosek (5000, 2000)

On Octave 8.4.0 it gave Inf without error.

On Octave 9.1.0 it gives error "GCD should be integer".

Anonymous

 

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

Attach Files:
   
   
Comment:
   

Attached Files
file #55873:  patch_65495 added by koerhen (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 jwe (Posted a comment)
  • -email is unavailable- added by koerhen (Updated the item)
  • -email is unavailable- added by mmuetzel
  • -email is unavailable- added by mmuetzel (Updated the item)
  • -email is unavailable- added by dasergatskov (Posted a comment)
  • -email is unavailable- added by nrjank (Posted a comment)
  •  

    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 14 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2024-04-21 mmuetzel StatusPatch Submitted Ready For Test
        Planned Release9.2.0 (current stable) 10.1.0 (current default)
    2024-04-13 rik5 Item GroupUnexpected Error or Warning Regression
        StatusConfirmed Patch Submitted
        Planned ReleaseNone 9.2.0 (current stable)
    2024-03-21 koerhen Attached File- Added patch_65495, #55873
    2024-03-20 mmuetzel Carbon-Copy- Added koerhen
    2024-03-20 mmuetzel StatusNeed Info Confirmed
        Release9.1.0 dev
        Operating SystemMicrosoft Windows Any
        Summarynchoosek error &quot;GCD should be an integer&quot; nchoosek error "gcd: all values must be integers"
    2024-03-20 nrjank CategoryPerformance Octave Function
    2024-03-20 nrjank StatusNone Need Info
        SummaryGCD error nchoosek error "GCD should be an integer"

    Back to the top

    Powered by Savane 3.13-caa5.
    Corresponding source code