failed nchoosek test in mercurial tip

[Francesco Potortì]

>>>>>>> processing /home/thorsten/hg/octave/scripts/specfun/nchoosek.m
>>  ***** assert (nchoosek(100,45), bincoeff(100,45))
>> !!!!! test failed
>> assert (nchoosek (100, 45),bincoeff (100, 45)) expected
>> 6.1448e+28
>> but got
>> 6.1448e+28
>> values do not match

On December, 1st, I posted a new changeset which corrects this.

>> Trying directly in octave, I get:
>> octave:2> (nchoosek (100, 45)-bincoeff (100, 45))/nchoosek(100,45)
>> ans =  2.8343e-14
>>
>> The difference certainly is not much, but still significantly larger  
>> than eps. So, the question is:
>> is it relevant?

Well, that depends on the purpose.  After the latest changeset I posted,
nchoosek should produce exact results, or else emit a warning.

>> regards
>>
>> Thorsten
>
>I see this as well. I've cc'd Francesco as he recently added this test.
>
>	http://hg.savannah.gnu.org/hgweb/octave/rev/cf620941af1a

Yes, that's the previous changeset.  The last one corrects this, but
maybe it has not been applied yet:
 <https://www-old.cae.wisc.edu/pipermail/octave-maintainers/2008-December/009672.html>

>Running Matlab 2008b
>
> >> fprintf('nchoosek=%15.f\n',nchoosek(100,45))
>Warning: Result may not be exact. Coefficient is greater than 1.000000e 
>+15
>  and is only accurate to 15 digits.
> > In nchoosek at 66
>nchoosek=61448471214136181560759549952

With the latest changeset applied I see:

octave> fprintf('nchoosek=%15.f\n',nchoosek(100,45))
warning: nchoosek: possible loss of precision
nchoosek=61448471214138292623084879872

Note that all those digits simply do not make sense, as only about 15
digits are significant:

octave> 1/eps
ans =  4.5036e+15

>If we trust the competitor, it appears that bincoeff is not as  
>accurate at nchoosek.

Yes, this is generally true, but not really important for this
particular case, as we are by far exceeding the available accuracy of
the machine (for 8-byte doubles, that is).

>Perhaps an appropriate change is
>
>< %!assert (nchoosek(100,45), bincoeff(100,45))
> > %!assert (1 - nchoosek(100,45) / bincoeff(100,45), 0, sqrt(eps))
>
>Thoughts?

The test that substitutes that one in the latest changeset I posted is 
%!assert (nchoosek(80,10), bincoeff(80,10))

-- 
Francesco Potortì (ricercatore)        Voice: +39 050 315 3058 (op.2111)
ISTI - Area della ricerca CNR          Fax:   +39 050 315 2040
via G. Moruzzi 1, I-56124 Pisa         Email: Potorti at isti.cnr.it
(entrance 20, 1st floor, room C71)     Web:   http://fly.isti.cnr.it/