[Octave-bug-tracker] [bug #62319] Another nchoosek bug for integer inputs

[Rik]

Follow-up Comment #5, bug #62319 (project octave):

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);
endfor


I then execute with


tic; bm_nchoosek; bm = toc


Results are:


New Code
3.911489963531494
3.924179077148438
3.934047937393188

Old Code
3.055461168289185
2.918849945068359
2.936843872070312


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)

    _______________________________________________________

Additional Item Attachment:

File name: bm_nchoosek.m                  Size:0 KB
    <https://file.savannah.gnu.org/file/bm_nchoosek.m?file_id=53098>



    _______________________________________________________

Reply to this item at:

  <https://savannah.gnu.org/bugs/?62319>

_______________________________________________
  Message sent via Savannah
  https://savannah.gnu.org/