Re: nchoosek checks, warning, corner cases and tests

[Jaroslav Hajek]

On Sat, Dec 13, 2008 at 3:17 AM, Francesco Potortì <address@hidden> wrote:
> Following the discussion about my previous changeset that disrupted some
> of Jaroslav's improvements, here is an alternative changeset that
> corrects that problem and improves the warning test for loss of
> precision.
>
> Please disregard the changeset that I had sent on December, 10th.
>
> Discussion:
> - help string now mentions some differences between this and bincoeff
> - many sanity checks are done on the arguments
> - handle k==0
> - correct k==n case, when a column vector was output
> - add error and warning tests
> - when testing against bincoeff, use a case where no precision is lost
>
> The changeset is both in the mail body and as an attachment.
>
> # HG changeset patch
> # User Francesco Potortì <address@hidden>
> # Date 1229133593 -3600
> # Node ID 1235e389bbd99059b370ddbfbc4e1c6e24268f34
> # Parent  87cca636a6c61e0cc8c3ab88a2652782203aff26
> nchoosek checks, warning, corner cases and tests
>
> diff -r 87cca636a6c6 -r 1235e389bbd9 scripts/ChangeLog
> --- a/scripts/ChangeLog Fri Dec 12 23:25:54 2008 +0100
> +++ b/scripts/ChangeLog Sat Dec 13 02:59:53 2008 +0100
> @@ -1,3 +1,8 @@ 2008-12-11  Jaroslav Hajek  <address@hidden
> +2008-12-13  Francesco Potortì  <address@hidden>
> +
> +       * specfun/nchoosek.m: Check for input arguments, signal loss of
> +       precision, correctly handle k==0 and k==n cases, add proper tests.
> +
>  2008-12-11  Jaroslav Hajek  <address@hidden>
>
>        * optimization/fsolve.m: Optionally allow pivoted qr factorization.
> diff -r 87cca636a6c6 -r 1235e389bbd9 scripts/specfun/nchoosek.m
> --- a/scripts/specfun/nchoosek.m        Fri Dec 12 23:25:54 2008 +0100
> +++ b/scripts/specfun/nchoosek.m        Sat Dec 13 02:59:53 2008 +0100
> @@ -49,31 +49,44 @@
>  ## of @var{n}, taken @var{k} at a time, one row per combination. The
>  ## resulting @var{c} has size @code{[nchoosek (length (@var{n}),
>  ## @var{k}), @var{k}]}.
> +##
> +## @code{nchoosek} works only for nonnegative integer arguments; use
> +## @code{bincoeff} for non-integer scalar arguments and for using vector
> +## arguments to compute many coefficients at once.
>  ##
>  ## @seealso{bincoeff}
>  ## @end deftypefn
>
>  ## Author: Rolf Fabian  <address@hidden>
>  ## Author: Paul Kienzle <address@hidden>
> -
> -## FIXME -- This function is identical to bincoeff for scalar
> -## values, and so should probably be combined with bincoeff.
> +## Author: Jaroslav Hajek
>
>  function A = nchoosek (v, k)
>
> -  if (nargin != 2)
> +  if (nargin != 2
> +      || !isnumeric(k) || !isnumeric(v)
> +      || !isscalar(k) || (!isscalar(v) && !isvector(v)))
>     print_usage ();
> +  endif
> +  if ((isscalar(v) && v < k) || k < 0
> +      || k != round(k) || any (v < 0 || v != round(v)))
> +    error ("nchoosek: args are nonnegative integers with V not less than K");
>   endif
>
>   n = length (v);
>
>   if (n == 1)
> -    if (k > v/2)
> -      k = v - k;
> +    k = min (k, v-k);          # improve precision at next step
> +    A = round (prod ((v-k+1:v)./(1:k)));
> +    if (A*2*k*eps >= 0.5)
> +      warning ("nchoosek", "nchoosek: possible loss of precision");
>     endif
> -    A = round (prod ((v-k+1:v)./(1:k)));
> +  elseif (k == 0)
> +    A = [];
>   elseif (k == 1)
>     A = v(:);
> +  elseif (k == n)
> +    A = v(:).';
>   elseif (k > n)
>     A = zeros (0, k, class (v));
>   else
> @@ -110,6 +123,8 @@ function A = nchoosek (v, k)
>   endif
>  endfunction
>
> -# nchoosek seems to be slightly more accurate (but only allowing scalar 
> inputs)
> -%!assert (nchoosek(100,45), bincoeff(100,45), -1e2*eps)
> +%!warning (nchoosek(100,45));
> +%!error (nchoosek(100,45.5));
> +%!error (nchoosek(100,145));
> +%!assert (nchoosek(80,10), bincoeff(80,10))
>  %!assert 
> (nchoosek(1:5,3),[1:3;1,2,4;1,2,5;1,3,4;1,3,5;1,4,5;2:4;2,3,5;2,4,5;3:5])
>
> --
> 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: address@hidden
> (entrance 20, 1st floor, room C71)     Web:   http://fly.isti.cnr.it/
>
>
>
>
>
>

I appied it.

thanks

-- 
RNDr. Jaroslav Hajek
computing expert
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz