[Octave-bug-tracker] [bug #61199] Allow 'char' input sets to nchoosek for Matlab compatibility

[anonymous]

Follow-up Comment #35, bug #61199 (project octave):

Rik, thank you for the profile information. I experimented with primes(k) and
factor() and different loop structures, which all turned out to be slower than
calling gcd() with the existing loops. This is likely because gcd() is
compiled while the others are interpreted. Replacing max() with find (..., 1)
was also not an improvement.

But it turns out that making num and denom have fewer elements by pairwise
multiplication does make it faster by some 15% to 20%. Even though the
individual elements are larger and each call to gcd() takes a little longer,
the number of elements is now about half what it was, and the number of
invocations of gcd() is reduced correspondingly. On my machine it falls from
280-285 microseconds to 230 microseconds. Still slower than the original code,
so it's accuracy at the expense of speed.

Patch is now as follows. Please test it and profile it.

  if (n == 1 && isnumeric (v))
    k = min (k, v-k);
    
    # make a list of integers for numerator and denominator,
    # then filter out common factors, and multiply what remains.
    # for a 15% to 20% performance boost, we multiply values pairwise so there
are fewer elements.
    if (rem(k,2)) # k is odd
      num = [(v-k+1:v-(k+1)/2) .* (v-(k-1)/2:v-1), v];
      denom = [(1:(k-1)/2) .* ((k+1)/2:k-1), k];
    else # k is even
      num = (v-k+1:v-k/2) .* (v-k/2+1:v);
      denom = (1:k/2) .* (k/2+1:k);
    end
    
    # remove common factors from num and denom
    while (~isempty(denom))
      for i = length(denom):-1:1
        tmp = gcd (denom(i), num);
        [m,f] = max(tmp);
        denom(i) /= m;
        num(f) /= m;
      endfor
      denom = denom(denom > 1);
      num = num(num > 1);
    endwhile

    C = prod (num);
    if (C > flintmax)
      warning ("nchoosek: possible loss of precision");
    end



    _______________________________________________________

Reply to this item at:

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

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