nchoosek checks, warning, corner cases and tests

[Francesco Potortì]

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])

nchoosek changeset
Description: "nchoosek changeset"

-- 
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/