bugGNU Octave - Bugs: bug #62319, Another nchoosek bug for integer...


bug #62319: Another nchoosek bug for integer inputs

Submitter:  Arun Giridhar <arungiridhar>
Submitted:  Fri 15 Apr 2022 03:20:21 PM UTC
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Regression
Status:  Fixed Assigned to:  None
Originator Name:  Open/Closed:  * Closed
Release:  * 7.1.0 Operating System:  * Any
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup

Jump to the original submission

Sat 16 Apr 2022 10:11:06 PM UTC, comment #6: 

@Rik in comment #5: sounds like a very good tradeoff. I had tried a similar if-check but instead of isinteger like your fix I was comparing the pairwise products to intmax and falling back to the slower code if necessary. The isinteger use is definitely faster, and avoids double-computation of numer and denom.

Arun Giridhar <arungiridhar>
Group Member
Sat 16 Apr 2022 07:45:49 PM UTC, comment #5: 

Library functions in Octave like nchoosek need to be accurate first, but performance is the next most important consideration.  Unlike user code that will get executed 10-1000 times, we can expect that library routines will be executed millions of times and hence even small improvements should be considered even at the cost of code complexity.

I wrote a simple benchmark script for nchoosek which I show below and is also attached.

N = 1e4;

for i = 1:N
  y = nchoosek (63,19);

I then execute with

tic; bm_nchoosek; bm = toc

Results are:

New Code

Old Code

The original comment in the code was correct that there is about a 25% savings from ~4 seconds to ~3 seconds.

Given that the majority of the time the input to nchoosek is going to be a floating point double (even though it won't have a fractional part), I don't think we should optimize for the integer input case.

Octave just needs to check the inputs to make sure that neither one is an integer before using the fast path code.  I checked in this change on the stable branch

diff -r 670eb988dd6a scripts/specfun/nchoosek.m
--- a/scripts/specfun/nchoosek.m        Fri Apr 15 14:46:23 2022 -0400
+++ b/scripts/specfun/nchoosek.m        Sat Apr 16 12:40:09 2022 -0700
@@ -115,8 +115,23 @@ function C = nchoosek (v, k)
     ## Steps: 1) Make a list of integers for numerator and denominator,
     ## 2) filter out common factors, 3) multiply what remains.
     k = min (k, v-k);
-    numer = (v-k+1):v;
-    denom = (1:k);
+    if (isinteger (v) || isinteger (k))
+      numer = (v-k+1):v;
+      denom = (1:k);
+    else
+      ## For a ~25% performance boost, multiply values pairwise so there
+      ## are fewer elements in do/until loop which is the slow part.
+      ## Since Odd*Even is guaranteed to be Even, also take out a factor
+      ## of 2 from numerator and denominator.
+      if (rem (k, 2))  # k is odd
+        numer = [((v-k+1:v-(k+1)/2) .* (v-1:-1:v-(k-1)/2)) / 2, v];
+        denom = [((1:(k-1)/2) .* (k-1:-1:(k+1)/2)) / 2, k];
+      else             # k is even
+        numer = ((v-k+1:v-k/2) .* (v:-1:v-k/2+1)) / 2;
+        denom = ((1:k/2) .* (k:-1:k/2+1)) / 2;
+      endif
+    endif

Performance is very slightly degraded by the extra isinteger() tests (~3%), but overall this restores the performance of the original code and yet it still passes the new BIST tests for saturating inputs.

(file #53098)

Rik <rik5>
Group administrator
Sat 16 Apr 2022 03:58:35 PM UTC, comment #4: 

It is slower by some 20%, yes, but it has only gone from about 0.1 milliseconds to 0.12 milliseconds per call when measured over some 5e4 calls. That's a 20-microsecond penalty per call to avoid the multiplication problems in calling nchoosek (int8 (32), 5) using the previous method.

Arun Giridhar <arungiridhar>
Group Member
Sat 16 Apr 2022 02:52:15 PM UTC, comment #3: 

Haven't looked into this in any detail. But I notice that the section of the code that was removed contained a comment that the prior implementation would save about 25% execution time.
Is that performance trick lost with this change? In case it is, is there an alternative implementation that maintains the previous performance but still avoids this issue?

Markus M├╝tzel <mmuetzel>
Group administrator
Fri 15 Apr 2022 06:58:12 PM UTC, comment #2: 

I pushed this fix to stable here: https://hg.savannah.gnu.org/hgweb/octave/rev/670eb988dd6a

Marking as Ready for test.

Arun Giridhar <arungiridhar>
Group Member
Fri 15 Apr 2022 03:24:03 PM UTC, comment #1: 

Here is the patch with the commit message.

I have added some BISTs that would cause a freeze without this code patch.

Marking as patch submitted.

(file #53092)

Arun Giridhar <arungiridhar>
Group Member
Fri 15 Apr 2022 03:20:21 PM UTC, original submission:  

This is to get a bug ID for an infinite loop with nchoosek for certain integer inputs.

Test: nchoosek (int8 (32), 5)

Numerically correct answer: 201376

Acceptable answer: saturate at intmax == 255 and give a warning so the user can switch to a different data type

Actual behavior: it freezes and needs control-C.

The problem is to do with the way the gcd calculation is set up, which causes an infinite loop for even modestly sized inputs when the integer bitwidth is small.

A fix will be attached with this bug ID.

Arun Giridhar <arungiridhar>
Group Member


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

Attach Files:

Attached Files
file #53098:  bm_nchoosek.m added by rik5 (54B - text/x-matlab)
file #53092:  nchoosek.patch added by arungiridhar (2KiB - text/x-patch)


Depends on the following items: None found

Items that depend on this one: None found


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

    Date Changed by Updated Field Previous Value => Replaced by
    2022-04-17 rik5 StatusReady For Test Fixed
        Open/ClosedOpen Closed
    2022-04-16 rik5 Attached File- Added bm_nchoosek.m, #53098
    2022-04-15 arungiridhar StatusPatch Submitted Ready For Test
    2022-04-15 arungiridhar Attached File- Added nchoosek.patch, #53092
        StatusNone Patch Submitted

    Back to the top

    Powered by Savane 3.13-7c56.
    Corresponding source code