[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
nchoosek warning, checks and tests
From: |
Francesco Potortì |
Subject: |
nchoosek warning, checks and tests |
Date: |
Wed, 10 Dec 2008 15:37:31 +0100 |
Here is a new version of my previous patch, this time against Jaroslav's
non recursive implementation.
Discussion:
- help string now mentions some differences between this and bincoeff
- many sanity checks are done on the arguments
- log of sum is used instead of product so that overflow is further away
- 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
Jaroslav, should I send you this as an attachment?
# HG changeset patch
# User Francesco Potortì <address@hidden>
# Date 1228919345 -3600
# Node ID 089d2f5463453acca07a0888fa9b6b52cf5cec86
# Parent c8a785d0e867a2f2c735a8f7b828c00e05204e59
nchoosek warning, checks and tests
diff -r c8a785d0e867 -r 089d2f546345 scripts/ChangeLog
--- a/scripts/ChangeLog Wed Dec 10 13:20:24 2008 +0100
+++ b/scripts/ChangeLog Wed Dec 10 15:29:05 2008 +0100
@@ -1,3 +1,8 @@ 2008-12-10 Jaroslav Hajek <address@hidden
+2008-12-10 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-10 Jaroslav Hajek <address@hidden>
* script/linear-algebra/expm.m: New source.
diff -r c8a785d0e867 -r 089d2f546345 scripts/specfun/nchoosek.m
--- a/scripts/specfun/nchoosek.m Wed Dec 10 13:20:24 2008 +0100
+++ b/scripts/specfun/nchoosek.m Wed Dec 10 15:29:05 2008 +0100
@@ -1,5 +1,5 @@
## Copyright (C) 2001, 2006, 2007 Rolf Fabian and Paul Kienzle
-## Copyright (C) 2008 Jaroslav Hajek
+## Copyright (C) 2008 Jaroslav Hajek and Francesco Potortì
##
## This file is part of Octave.
##
@@ -49,31 +49,45 @@
## 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
+## 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 ((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 (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
- 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 +124,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/
- nchoosek warning, checks and tests,
Francesco Potortì <=