Binomial coefficient, Octave

Language: Layout:

function ret = timeit(f)
	st = tic;
	f();
	ret = toc(st);
end
function nk=nchoosek2(n, k)
	if n-k > k
		nk = prod((k+1:n) .* prod((1:n-k).^ (-1/(n-k))));
	else
		nk = prod((n-k+1:n) .* prod((1:k).^ (-1/k)) ) ;
	end
end
warning ('off');
  function P=test_binom()
      k1=150; k2=150; D=200; P=0;
      for i=0:D-1
          for j=0:i
              if (i-j>k2||j>k1) 
                  continue;
              end
              P=P+nchoosek2(k1,j)*nchoosek2(k2,i-j)/nchoosek2((k1+k2),i);          
          end 
      end
  end
function ret = choose (n,k, D, K1, K2)
	persistent binTable = zeros(max([D+1,K1+K2+1]) , D+1);
	if binTable(n+1,k+1) == 0
		binTable(n+1,k+1) = nchoosek2(n,k);
	end
	ret = binTable(n+1,k+1);
end
function P = tst()
    P=0;k1=150; k2=150; D=200; P=0;
    choose(1,0,D,k1,k2);
    for i = 0:D-1
        for j = j=max(i - k2 , 0):min (i,k1)
            P=P+choose(k1,j)*choose(k2,i-j)/choose((k1+k2),i);
        end
    end
end

function C = nCk(n,k)
% use smaller k if available
k(k>n/2) = n-k(k>n/2);
k = k(:);
kmat = ones(numel(k),1)*(1:max(n-k));
kmat = kmat.*bsxfun(@le,kmat,(n-k));
pw = bsxfun(@power,kmat,-1./(n-k));
pw(kmat==0) = 1;
kleft = ones(numel(k),1)*(min(k):n);
kleft = kleft.*bsxfun(@gt,kleft,k);
t = bsxfun(@times,kleft,prod(pw,2));
t (kleft==0) = 1;
C = prod(t,2);
end
function P = binomial_coefficient(k1,k2,D)
warning ('off','MATLAB:nchoosek:LargeCoefficient');
i_ind = nonzeros(triu(ones(D,1)*(1:D)))-1;
j_ind = nonzeros(tril(ones(D,1)*(1:D+1)).')-1;
valid = ~(i_ind-j_ind>=k2 | j_ind>=k1);
i_ind = i_ind(valid);
j_ind = j_ind(valid);
beta = @(ii,jj) nCk(k1,jj).*nCk(k2,ii-jj)./nCk((k1+k2),ii);
b = beta(i_ind,j_ind);
P = sum(b(:));
end

function P = binomial_coefficient_norepeat(k1,k2,D)
	if ~exist('repelem')
		repelem = @(v,n) repelems(v, [1:numel(v); n]);
	end
	warning ('off','MATLAB:nchoosek:LargeCoefficient');
	i_ind = nonzeros(triu(ones(D,1)*(1:D)))-1;
	j_ind = nonzeros(tril(ones(D,1)*(1:D+1)).')-1;
	valid = ~(i_ind-j_ind>=k2 | j_ind>=k1);
	i_ind = i_ind(valid);
	j_ind = j_ind(valid);
	ni= numel(i_ind);
	sb_r = repelem([k1 k2 k1+k2], repmat(ni,1,3))';
	sb_c = [j_ind; i_ind-j_ind; i_ind];
	ranges = 1:numel(sb_c);
	idx =sub2ind([k1+k2 , max(sb_c)+1], sb_r, sb_c+1);
	[idx_srtd idx_srtd_i] = sort(idx);
	[uq_idx uq_idx_i] = unique(idx_srtd);
	hst = [uq_idx_i(1); diff(uq_idx_i)];
	[~, bck_idx] = sort(idx_srtd_i);
	idx_uq_idx = ranges(idx_srtd_i(uq_idx_i));
	nck_args = [sb_r(idx_uq_idx) sb_c(idx_uq_idx)];
	[arg_srtd arg_srtd_i] = sortrows(nck_args);
	[uq_arg uq_arg_i] = unique(arg_srtd(:,1));
	cnt = [0 ;cumsum([uq_arg_i(1); diff(uq_arg_i)])];
	tbl = zeros(size(nck_args,1),1);
	for i=2:numel(cnt)
		tbl(cnt(i-1)+1:cnt(i),:) = nCk(uq_arg(i-1),arg_srtd(cnt(i-1)+1:cnt(i),2));
	end
	[~, bck_arg] = sort(arg_srtd_i);
	tbl = tbl(bck_arg,:);
	tblrep = repelem(tbl,hst');
	tblrep = tblrep(bck_idx);
	beta = @(arg_srtd_i,jj) tblrep(1:ni) .* tblrep(ni +1 : 2*ni) ./ tblrep(2* ni +1 : 3*ni);
	b = beta(i_ind,j_ind);
	P = sum(b(:));
end

test_binom()
tst()
k1=150; k2=150; D=200;
binomial_coefficient(k1,k2,D)
binomial_coefficient_norepeat(k1,k2,D)
timeit(@test_binom)
timeit(@tst)

timeit(@()binomial_coefficient(k1,k2,D))
timeit(@()binomial_coefficient_norepeat(k1,k2,D))

[ + ] Show input
Absolute running time: 9.55 sec, cpu time: 9.38 sec, memory peak: 268 Mb, absolute service time: 9,56 sec