bugGNU Octave - Bugs: bug #61565, nchoosek broken for integer inputs

 
 

bug #61565: nchoosek broken for integer inputs

Submitter:  None
Submitted:  Fri 26 Nov 2021 04:02:53 PM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Regression
Status:  Fixed Assigned to:  None
Originator Name:  Originator Email:  -email is unavailable-
Open/Closed:  * Closed Release:  * 7.0.90
Release:  Operating System:  * Any
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Mon 06 Dec 2021 12:35:24 PM UTC, comment #21: 

Seems to be working now.

Closing as fixed.

Markus Mützel <mmuetzel>
Group administrator
Thu 02 Dec 2021 02:13:18 PM UTC, comment #20: 

Oops. Forgot that this still depends on bug #61132. While that isn't fixed, the new BIST will be failing.

Markus Mützel <mmuetzel>
Group administrator
Thu 02 Dec 2021 01:51:02 PM UTC, comment #19: 

I pushed a modified version of your patch here:
https://hg.savannah.gnu.org/hgweb/octave/rev/391b35ef8b24

Not sure if performance is an issue here, but I moved the check for (Octave's default) floating point numbers to the front. Additionally, I tried to make the warning messages more concise. Lastly, I added a BIST for the new warning.

Marking as ready for test.

Markus Mützel <mmuetzel>
Group administrator
Thu 02 Dec 2021 01:08:32 PM UTC, comment #18: 

I've attached a warning patch. I see that your edits to nchoosek are on stable and my diff is against default. If it doesn't apply cleanly, please merge stable to default and I will reupload the warning patch.

(file #52406)

Arun Giridhar <arungiridhar>
Group Member
Thu 02 Dec 2021 10:22:34 AM UTC, comment #17: 

Looks good to me. Maybe try to break long lines (at 80 characters) and add a warning id to both warnings so user could turn them off in case they are bothered by them.

Could you please prepare a patch with these changes?

Markus Mützel <mmuetzel>
Group administrator
Thu 02 Dec 2021 09:59:32 AM UTC, comment #16: 

Sounds fair. This warning code works for me. Please edit as needed:


C = prod (numer, "native");
if (isinteger (v) && C >= intmax (class (v)))
  warning ("nchoosek: result is saturated at intmax. Cast inputs to double if you need more range.");
elseif (isfloat (v) && C >= flintmax (class (v)))
  warning ("nchoosek: result exceeds flintmax and may not be exact.")
end


Arun Giridhar <arungiridhar>
Group Member
Thu 02 Dec 2021 09:50:27 AM UTC, comment #15: 

Well, I'd find it more surprising if the output type "suddenly" changed depending on the magnitude...
If a user is interested in a floating point approximation, they could explicitly cast their input to a floating point type. Otherwise, we should assume that the user knows what they are doing. But a warning would definitely help.

Markus Mützel <mmuetzel>
Group administrator
Thu 02 Dec 2021 09:25:21 AM UTC, comment #14: 

Oh wow. I was not expecting them to saturate at intmax without a warning.

For integer inputs, should we also saturate at intmax of the same class as the inputs, or should we return a doublle that is at least in the correct order of magnitude? I think the second is more useful and much less surprising to the user, but would like feedback on whether it would break Matlab compatibility.

Thinking of something like this: calculate C = prod (numer, "native"), and if it saturates at the corresponding intmax or exceeds flintmax("single"), depending on the class of input, then we can recalculate it as a double and return it with a warning that it is beyond intmax or the single-precision flintmax. If the inputs are double already, we retain the current check for the double-precision flintmax. Thoughts on this?

Arun Giridhar <arungiridhar>
Group Member
Thu 02 Dec 2021 09:03:05 AM UTC, comment #13: 

Forgot the second part.
Also in Matlab R2021a:

>> nchoosek (uint64 (80), uint64 (40))

ans =

  uint64

   18446744073709551615


No warning either. We could probably do a better job and emit a warning if the output saturated to `intmax`.

Markus Mützel <mmuetzel>
Group administrator
Thu 02 Dec 2021 09:00:53 AM UTC, comment #12: 

In Matlab R2021a:

>> nchoosek (uint64 (60), uint64 (30))

ans =

  uint64

   118264581564861424


No warning (which makes sense).

We should probably tweak that condition. And also use the correct `flintmax` depending on whether the output `C` is double or single precision.

Markus Mützel <mmuetzel>
Group administrator
Thu 02 Dec 2021 08:15:30 AM UTC, comment #11: 

Does Matlab still give that warning for nchoosek (uint64 (60), uint64 (30)) and for nchoosek (uint64 (80), uint64 (40))? If so, then no further changes are needed for Octave's nchoosek, otherwise the threshold may need to be tweaked.

Arun Giridhar <arungiridhar>
Group Member
Thu 02 Dec 2021 01:05:06 AM UTC, comment #10: 

Matlab 2021a:


>> nchoosek (60, 30)
Warning: Result may not be exact. Coefficient is greater than
9.007199e+15 and is only accurate to 15 digits
> In nchoosek (line 92)
ans =
     1.182645815648614e+17
>> nchoosek (80, 40)
Warning: Result may not be exact. Coefficient is greater than
9.007199e+15 and is only accurate to 15 digits
> In nchoosek (line 92)
ans =
     1.075072087333362e+23


Nicholas Jankowski <nrjank>
Group Member
Wed 01 Dec 2021 10:28:30 PM UTC, comment #9: 

That is even better if Matlab returns an integer. I was basing my statement on what I read in this bug report about nchoosek: https://savannah.gnu.org/bugs/index.php?61199#comment25

Could you please check Matlab whether it warns of precision for values more than flintmax? If not, that can be removed or changed from nchoosek too. Good tests are: nchoosek (60, 30) gives a value above flintmax but less than 2^64, so passing both inputs as uint64 will help. Another is nchoosek (80, 40) which is well above 2^64, so any precision warning should be triggered in Matlab.

Arun Giridhar <arungiridhar>
Group Member
Wed 01 Dec 2021 09:27:18 PM UTC, comment #8: 

I checked in Matlab R2021a:

>> nchoosek(uint8(10), uint8(5))

ans =

  uint8

   252


Which version of Matlab did you use that returns a double value for that input?

Markus Mützel <mmuetzel>
Group administrator
Wed 01 Dec 2021 09:21:11 PM UTC, comment #7: 

No, the colon operator has a lower precedence than the arithmetic ones, so the parentheses don't make a difference.

It turned out the first nchoosek-specific bug was in this line, where k was being rounded differently when it was an integer (rounded up for odd integers) compared to a float (always rounded down):

    denom = [(1:k/2) .* (k-1:-1:(k+1)/2) / 2, k];


That line needed to be changed to:

    denom = [(1:(k-1)/2) .* (k-1:-1:(k+1)/2) / 2, k];


The reason this bug was not exposed until now was that the previous line in nchoosek would trigger an error about ranges, so we didn't get to see this behavior until the range bugs were fixed first. With the change above, nchoosek (uint8(10), uint8(5)) evaluates correctly, and it is no longer a regression from Octave 6.

But that change is not enough, because nchoosek still complains about mismatched types when given different integer types:

>> nchoosek (uint8(11), int8(5))
error: binary operator '-' not implemented for 'uint8 scalar' by 'int8 scalar' operations
error: called from
    nchoosek at line 117 column 7


This problem is internal to the integer type system and not because of the definition of the binomial coefficient. The simplest solution is to cast both v and k to doubles. Losing precision for v and k above flintmax is not a problem because nchoosek is required (for Matlab compatibility) to return a double and any result that exceeds flintmax will trigger a warning anyway. (Casting v and k to double does not introduce new loss of precision.)

So the two changes to be made are to convert k/2 to (k-1)/2 in the line above, and to convert v and k to double. Casting to double makes the first change optional.

Arun Giridhar <arungiridhar>
Group Member
Wed 01 Dec 2021 09:20:32 PM UTC, comment #6: 

I pushed a change here that allows using this function with integer type input:
https://hg.savannah.gnu.org/hgweb/octave/rev/181c44079f6c

This still depends on the changes for bug #61132.

Afaict, this should be Matlab compatible. They also fail for mixed integer type input.

Marking as ready for test.

Markus Mützel <mmuetzel>
Group administrator
Wed 01 Dec 2021 07:33:26 AM UTC, comment #5: 

Would adding bracket in these calculations help to make sure that all intermediate results stay integer?

Markus Mützel <mmuetzel>
Group administrator
Tue 30 Nov 2021 10:32:40 PM UTC, comment #4: 

Note: Depending on the parity of k, and whether n and k are signed or unsigned, the errors can either be about descending ranges (being fixed in bug #61132) or a length mismatch which is caused by rounding integers differently from floating point when divided by 2. This was exposed by recent patches to bug #61300.


    121     ## Since Odd*Even is guaranteed to be Even, also take out a factor
    122     ## of 2 from numerator and denominator.
    123     if (rem (k, 2))  # k is odd
    124       numer = [(v-k+1:v-(k+1)/2) .* (v-1:-1:v-(k-1)/2) / 2, v];
    125       denom = [(1:k/2) .* (k-1:-1:(k+1)/2) / 2, k];
    126     else             # k is even
    127       numer = (v-k+1:v-k/2) .* (v:-1:v-k/2+1) / 2;
    128       denom = (1:k/2) .* (k:-1:k/2+1) / 2;
    129     endif


Which means that independent of fixes to bug #61132 or bug #61300, we do need to convert n and k to double inside nchoosek before calculating the binomial coefficient otherwise they will round differently for float and int and cause a length mismatch. Patch given in file #52354 does that and adds a BIST.

Arun Giridhar <arungiridhar>
Group Member
Fri 26 Nov 2021 05:15:22 PM UTC, comment #3: 

@OP, thank you for reporting this.

@jwe: Independent of the integer ranges in bug #61132, nchoosek can use the freedom to mix and match types like this:


octave:1> nchoosek (uint8 (10), uint8 (5))
ans = 252
octave:2> nchoosek (uint8 (10), uint16 (5))
ans = 252
octave:3> nchoosek (uint8 (10), uint64 (5))
ans = 252
octave:4> nchoosek (uint8 (10), int64 (5))
ans = 252
octave:5> nchoosek (int8 (10), int64 (5))
ans = 252


This is not currently possible, with errors like this, which will persist even if the integer range behavior is changed:


octave:1> nchoosek (int8 (10), int64 (5))
error: binary operator '-' not implemented for 'int8 scalar' by 'int64 scalar' operations
error: called from
    nchoosek at line 119 column 7


I have attached a small patch for nchoosek, with a BIST. It converts both v and k to double for the binomial calculation. This allows for calling it with all combinations of integers and floating points like the above, also solving OP's problem. Please review and edit as required.

(file #52354)

Arun Giridhar <arungiridhar>
Group Member
Fri 26 Nov 2021 04:24:03 PM UTC, comment #2: 

Adding "v = double (v);" just before the line "k = min (k, v-k);" fixes the problem.

Anonymous
Fri 26 Nov 2021 04:23:58 PM UTC, comment #1: 

Thanks for the report.  The root cause of this problem has been reported in bug #61132

and also discussed in this thread on the discourse forum:

https://octave.discourse.group/t/improving-correctness-and-reducing-maintenance-burden/1605/4

I've marked bug #61132 as a blocker for the version 7 release.

John W. Eaton <jwe>
Group administrator
Fri 26 Nov 2021 04:02:53 PM UTC, original submission:  

Hello Octave,

On Octave 7, nchoosek fails for integer inputs:


octave:1> nchoosek (uint8 (10), uint8 (5))
error: colon operator increment invalid (not an integer or out of range for given integer type)
error: called from
    nchoosek at line 124 column 13


On Octave 6 at https://octave-online.net/, the result is correct:


octave:1> nchoosek (uint8 (10), uint8 (5))
ans = 252


Anonymous

 

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

Attach Files:
   
   
Comment:
   

Attached Files
file #52354:  nchoosek.patch added by arungiridhar (853B - text/x-patch)

 

Digest:
   bug dependencies.

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by nrjank (Posted a comment)
  • -email is unavailable- added by mmuetzel (Posted a comment)
  • -email is unavailable- added by arungiridhar (Updated the item)
  • -email is unavailable- added by jwe (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 8 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2021-12-06 mmuetzel StatusReady For Test Fixed
        Open/ClosedOpen Closed
    2021-12-02 mmuetzel StatusNone Ready For Test
    2021-12-02 arungiridhar Attached File- Added nchoosek_warning.patch, #52406
    2021-12-02 mmuetzel StatusReady For Test None
    2021-12-01 mmuetzel StatusNone Ready For Test
        Dependencies- Depends on bugs #61132
    2021-11-26 arungiridhar Attached File- Added nchoosek.patch, #52354

    Back to the top

    Powered by Savane 3.15-e6e5.
    Corresponding source code