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

[anonymous]

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

Change this in the existing code:

  if (n == 1 && isnumeric (v))
    ## Improve precision at next step.
    k = min (k, v-k);
    C = round (prod ((v-k+1:v)./(1:k)));
    if (C*2*k*eps >= 0.5)
      warning ("nchoosek: possible loss of precision");
    endif
  elseif (k == 0)


to this:

  if (n == 1 && isnumeric (v))
    k = min (k, v-k);
    num = (v-k+1):v;
    denom = 2:k; # no need for 1

    # 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
  elseif (k == 0)


in order to achieve this difference:

>> nchoosek (58,29)   ### PATCHED
ans = 30067266499541040

>> nchoosek (58,29)   ### UNPATCHED
warning: nchoosek: possible loss of precision
warning: called from
    nchoosek at line 118 column 7

ans = 3.006726649954106e+16


If it is acceptable to cast into uint64, we can change the end of the patch as
follows:

    C = prod (num);
    if (C > flintmax)
      if (C < intmax("uint64")) # can cast up to get precision
        C = prod (uint64(num), "native");
      else # can't get precision from casting up
        warning ("nchoosek: possible loss of precision");
      end
    end


That will go all the way to 2^64-1 but will return a uint64 instead of a
double for values bigger than flintmax but less than intmax("uint64").

>> nchoosek (67,33)
ans = 14226520737620288370
>> class(nchoosek (67,33))
ans = uint64



    _______________________________________________________

Reply to this item at:

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

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