[Octave-bug-tracker] [bug #35611] Bug in bincoeff with large n and k

[Rik]

Follow-up Comment #2, bug #35611 (project octave):

Check the documentation for nchoosek and bincoeff.  They compute the same
mathematical value, but the underlying algorithms are very different.  In
particular, nchoosek is optimized for integer inputs and gives the correct
answer in both Octave and Matlab.  bincoeff is a much broader function and can
handle real inputs as well as negative inputs.

Quoting from the documentation for bincoeff:


In most cases, the `nchoosek' function is faster for small scalar
integer arguments.  It also warns about loss of precision for big
arguments.


Quoting from the documentation for nchoosek:


`nchoosek' works only for non-negative, integer arguments.  Use
`bincoeff' for non-integer and negative scalar arguments, or for
computing many binomial coefficients at once with vector inputs
for N or K.


nchoosek (100,10) = 1.7e13 which is representable exactly as an integer in
doubele precision arithmetic (you can check bitmax() which is 9e15 for
doubles).  Thus, nchoosek is really the correct function to use for this
example.

The line to take a look at is 98 of bincoeff.m which I reproduce below.


## Clean up rounding errors.
b(n_int) = round (b(n_int));


Before the rounding, the value of b is 17310309456440.8.  In this case, if we
used truncation, fix(), rather than round we would get the correct answer. 
However, I'm note sure this is a general result.  I believe round was chosen
so that 50% of the time we would guess low and 50% of the time we would guess
high and on average the algorithm would be more accurate.  It would be
valuable if you could prove any relation that might exist between the input
and the fractional part of the bincoeff algorithm.  For example, is it always
the case for integers less than bitmax that the result should use fix() rather
than round()?



    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/bugs/?35611>

_______________________________________________
  Message sent via/by Savannah
  http://savannah.gnu.org/