[Octave-bug-tracker] [bug #65495] nchoosek error "gcd: all values must be integers"

[Hendrik K]

Follow-up Comment #13, bug #65495 (group octave):

Interesting case:

a)
One issue is that C may become Inf within the loop, so the gcd function
(gcd(Inf)) complains and throws the observed error. This is easy to check and
then break the loop within nchoosek accordingly.
Note: The gcd error message could be refined when receiving Inf compared to
receiving e.g. a fraction, though.

b)
The other issue is that the gcd function does NOT always provide a correct
result for floating point numbers beyond flintmax.
Example:
gcd (uint64(879236959347161415), 9)
-> 9
gcd (double(879236959347161415), 9)
-> 3

This impacts nchoosek: In theory, the final division by "i/(g1 * g2)" is NOT
necessary, as it should be one by design. Unfortunately this is not always the
case for floating point numbers as shown above (integer numbers have never
these issues).

Once C times (v - k + i) reaches i_max and before the the (next) C iteration
(may) reach Inf, we continue to work with the gcd for numerical precision
reasons in each iteration. Or as the gcd function value for floats could be
too small, we may temporarily reach Inf before dividing by "i / (g1 * g2)"
which would not be identical to one in this case.

Below is the proposed changed code


    for i = 1:k
      if (C * (v - k + i) >= imax)
        ## Avoid overflow / precision loss by determining the smallest
        ## possible factor of (C * (n-k+i)) and i via the gcd.
        ## Note that by design in each iteration
        ##   1) C will always increase (factor is always > 1).
        ##   2) C will always be a whole number.
        ## Therefore, using the gcd will always provide the best possible
        ## solution until saturation / has the least precision loss.
        g1 = gcd (C, i);
        g2 = gcd (v - k + i, i/g1);
        C /= g1;

        ## In theory and (always for integers) i/(g1 * g2) is identical to 1
by
        ## design. Or for floats and beyond flintmax, the gcd may not be
        ## correctly derived by the gcd function and i/(g1 * g2) may not be
1.
        if (i/(g1 * g2) == 1) || !isinf(C * (v - k + i)/g2)
          C *= (v - k + i)/g2;
          C /= i/(g1 * g2);
        else
          C /= i/(g1 * g2);
          ## We have potential precision loss by dividing (too) early, but
          ## advantage is that we prevent possible interim overflows
          C *= (v - k + i)/g2;
        endif
        if (is_int && (C == imax)) || (! is_int && isinf(C))
          break;  # Stop here; saturation reached.
        endif
      else
        C *= (v - k + i);
        C /= i;
      endif
    endfor
 

  

@Markus-
you were right: When changing the order of multiplication and division, we are
expected to get a precision loss in certain cases, so we only do this when it
is absolutely necessary in the proposed code change.

 " || !isinf(C * (v - k + i)/g2))




(file #55873)

    _______________________________________________________

Additional Item Attachment:

File name: patch_65495                    Size:1 KB
    <https://file.savannah.gnu.org/file/patch_65495?file_id=55873>


    AGPL NOTICE

These attachments are served by Savane. You can download the corresponding
source code of Savane at
https://git.savannah.nongnu.org/cgit/administration/savane.git/snapshot/savane-cf05868ce827c9612cb6083f86556077b7ab06b3.tar.gz


    _______________________________________________________

Reply to this item at:

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

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