[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: nchoosek warning, checks and tests
From: |
Jaroslav Hajek |
Subject: |
Re: nchoosek warning, checks and tests |
Date: |
Wed, 10 Dec 2008 15:50:42 +0100 |
On Wed, Dec 10, 2008 at 3:37 PM, Francesco Potortì <address@hidden> wrote:
> 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");
I think this should say "args should be...", not "args are...",
because they aren't, if we got to this point.
> 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)));
Why did you replace
round (prod ((v-k+1:v)./(1:k)))
by
round (exp (sum (log (v-k+1:v)) - sum (log (2:k))))
?
Using your numbers, n = 100, k = 45, I get (format long)
6.14484712141362e+28 from the prod version,
6.14484712141383e+28 from the exp-sum-log version.
Comparing that to the exact value 61448471214136179596720592960 (comp. by PARI),
not only is the first expression faster, it appears also more
accurate. So, what is the reason?
> + 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/
>
--
RNDr. Jaroslav Hajek
computing expert
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz