failed nchoosek test in mercurial tip

[Ben Abbott]

On Dec 2, 2008, at 3:34 PM, Thorsten Meyer wrote:

> Hi,
>
> while doing
> 	make check
> with the latest mercurial tip on an intel pc running debian, I got  
> this test failure:
>
>>>>>> 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
>
> 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?
>
> 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

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

 From both Octave 3.0.3 and my build from this morning

octave:1> fprintf('nchoosek=%d\nbincoeff=%d 
\n',nchoosek(100,45),bincoeff(100,45))
nchoosek=61448471214136542200573460480
bincoeff=61448471214133050151643643904

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

Perhaps an appropriate change is

< %!assert (nchoosek(100,45), bincoeff(100,45))
 > %!assert (1 - nchoosek(100,45) / bincoeff(100,45), 0, sqrt(eps))

Thoughts?

Ben