nchoosek improvements and cleanup

[Francesco Potortì]

While cleaning up the docs and checking the args, I discovered that the
used algorithm was suboptimal in most common cases, that is, when the
second arg is less than half the first arg.  In these cases, the new
algorithm is faster, uses less memory and is numerically more accurate.

Forgot to mention in the Changelog: I also added a check for the cases
when there is loss of precision and a warning.  I was inspired for this
by the Matlab help page, but my test is possibly more accurate, as theirs
just checks that the result does not exceed 15 digits.

# HG changeset patch
# User Francesco Potortì <address@hidden>
# Date 1228142575 -3600
# Node ID 021be2026afac53045210fb2669331fbe6414079
# Parent  ed030a109bb8bfbeb4def31ac84eebbb975ffa30
nchoosek improvements and cleanup

diff -r ed030a109bb8 -r 021be2026afa scripts/ChangeLog
--- a/scripts/ChangeLog Wed Nov 26 15:43:21 2008 +0100
+++ b/scripts/ChangeLog Mon Dec 01 15:42:55 2008 +0100
@@ -1,3 +1,10 @@ 2008-11-26  Francesco Potortì  <address@hidden
+2008-12-01  Francesco Potortì  <address@hidden>
+
+       * specfun/nchoosek.m: Add checks for numeric, integer arguments
+       with correct bounds and relative tests.
+       Remove the FIXME notice and add some difference info wrt bincoeff.
+       Improve precision, execution time and memory used in common cases.
+
 2008-11-26  Francesco Potortì  <address@hidden>
 
        * specfun/nchoosek.m: Set max_recursion_depth and use a subfunction.
diff -r ed030a109bb8 -r 021be2026afa scripts/specfun/nchoosek.m
--- a/scripts/specfun/nchoosek.m        Wed Nov 26 15:43:21 2008 +0100
+++ b/scripts/specfun/nchoosek.m        Mon Dec 01 15:42:55 2008 +0100
@@ -1,4 +1,5 @@
 ## Copyright (C) 2001, 2006, 2007 Rolf Fabian and Paul Kienzle
+## Copyright (C) 2008 Francesco Potortì
 ##
 ## This file is part of Octave.
 ##
@@ -49,25 +50,37 @@
 ## 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: Francesco Potortì <address@hidden>
 
 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 (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)
-     A = round (exp (sum (log (k+1:v)) - sum (log (2:v-k))));
+    k = min (k, v-k);          # improve precision at next step
+    A = round (exp (sum (log (v-k+1:v)) - sum (log (2:k))));
+    if (A > 1/(exp (2*k*eps) - 1))
+      warning ("nchoosek", "nchoosek: possible loss of precision");
+    endif
   elseif (k == 0)
     A = [];
   elseif (k == 1)
@@ -96,5 +109,8 @@ function A = nck (v, k)
   endif
 endfunction
 
-%!assert (nchoosek(100,45), bincoeff(100,45))
+%!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/