From octave-sources-request at bevo dot che dot wisc dot edu Fri Aug 17 05:38:33 2001
Subject: binominal coefficients ( nchoosek.m )
From: Rolf Fabian <fabian at tu-cottbus dot de>
To: "'Paul Kienzle'" <pkienzle at kienzle dot powernet dot co dot uk>
Cc: "'octave-sources UWISC'" <octave-sources at bevo dot che dot wisc dot edu>
Date: Fri, 17 Aug 2001 12:34:32 +0200

Hi Paul,

enclosed please find my version of a binominal coefficient
function 'nchoosek(n|v,k)', which is based on your 'nchoosek.m'
octave-script, released with  matcompat-2001-02-25 packkage.

Besides removal of a minor bug which may occur for
complex vector v input ( your orig. version
invokes a complex conjugated transposition, whereas
ONLY transposition is needed ), my version introduces
considerable extension to the acceptable range of
input variables.

Your method of calculation uses ratios of products like

----------------------------------------------------------
(*)	c = prod( k+1:n ) / prod( 2:n-k )
----------------------------------------------------------

e.g.  at  k=3, this leads to an overflow of the nominator already
for relatively small n=171 . The following sample results are
calculated with your original function-script:

:) nchoosek_org(170,3)
ans =  8.0444e+05
:) nchoosek_org(171,3)
ans = Inf
:) nchoosek_org(172,3)
ans = Inf
:) nchoosek_org(173,3)
ans = Inf
:) nchoosek_org(174,3)
ans = NaN
:) nchoosek_org(200,8)
ans = NaN

On the other hand, my approach for expression (*) in order
to prevent overflow for such relatively small arguments is
the usage of differences of sums of logs instead of a ratio
of large numbers:

----------------------------------------------------------
c = round( exp( sum( log(k+1:n) ) - sum( log( 2:n-k ) ) ))
----------------------------------------------------------

The following resultss are calculated with the
attached version of the 'nchoosek'-function, which
implements this method:

:)nchoosek(170,3)
ans = 804440
:)nchoosek(171,3)
ans = 818805
:)nchoosek(172,3)
ans = 833340
:)nchoosek(173,3)
ans = 848046
:)nchoosek(174,3)
ans = 862924
:)nchoosek(200,8)		% format long 
ans = 55098996177226

This approach allows successfull calculation
even for relatively large n  e.g. 
:)nchoosek(1e3,100)
ans =  6.38505119264157e+139
	-or-
:)nchoosek(1e3,500)
ans =  2.70288240944523e+299
	-or-
:)nchoosek(100000,50)
ans =  3.24791116453536e+185

Hence, such a numerical approach provides
a large performance enhancement to the
original function found in the current
matcompat-package.

Best wishes
Rolf

------------------------ nchoosek.m ---------------------------
--------------------------- snip ------------------------------
## Copyright (C) 2001 Rolf Fabian, Paul Kienzle
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 2 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program; if not, write to the Free Software
## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

##USAGE	c = nchoosek ( n,k )
##	return binominal coefficient  c = n! / (k! * (n-k)!)
##
##	A = nchoosek ( v,k )
##	generate all combinations of the elements of vector v
##	taken k at a time, one row per combination. The resulting
##	A has size [nchoosek(n,k), k], where n is the length of v
##
##AUTHORS Rolf Fabian  <fabian at tu-cottbus dot de>
##	  Paul Kienzle <pkienzle at kienzle dot powernet dot co dot uk>

function A = nchoosek(v,k)

  n = length(v);

  if n == 1
     A = round(exp( sum(log(k+1:v)) - sum(log(2:v-k)) ));

  elseif k == 0
    A = [];

  elseif k == 1
    A = v(:);

  elseif k == n
     A = v(:).';

  else

    m =  round(exp( sum(log(k:n-1)) - sum(log(2:n-k)) ));
    A = [ v(1)*ones( m,1 ), nchoosek(v(2:n), k-1) ; \
    	  nchoosek(v(2:n), k) ];
  endif

endfunction