[dbinomresult,n,m]=dbinom(n,m);

function [dbinomresult,n,m]=dbinom(n,m);
dbinomresult=[];
persistent bilnmx corr fintmx first firstCall i k sq2pil xk xn xnk ; if isempty(firstCall),firstCall=1;end; 

;
if isempty(i), i=0; end;
if isempty(k), k=0; end;
%***BEGIN PROLOGUE  DBINOM
%***PURPOSE  Compute the binomial coefficients.
%***LIBRARY   SLATEC (FNLIB)
%***CATEGORY  C1
%***TYPE      doubleprecision (BINOM-S, DBINOM-D)
%***KEYWORDS  BINOMIAL COEFFICIENTS, FNLIB, SPECIAL FUNCTIONS
%***AUTHOR  Fullerton, W., (LANL)
%***DESCRIPTION
%
% DBINOM(N,M) calculates the doubleprecision binomial coefficient
% for integer arguments N and M.  The result is (N!)/((M!,N-M)!).
%
%***REFERENCES  (NONE)
%***ROUTINES CALLED  D1MACH, D9LGMC, DLNREL, XERMSG
%***REVISION HISTORY  (YYMMDD)
%   770601  DATE WRITTEN
%   890531  Changed all specific intrinsics to generic.  (WRB)
%   890531  REVISION DATE from Version 3.2
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
%***end PROLOGUE  DBINOM
if isempty(corr), corr=0; end;
if isempty(fintmx), fintmx=0; end;
if isempty(sq2pil), sq2pil=0; end;
if isempty(xk), xk=0; end;
if isempty(xn), xn=0; end;
if isempty(xnk), xnk=0; end;
if isempty(bilnmx), bilnmx=0; end;
if isempty(first), first=false; end;
if firstCall,   sq2pil=[0.91893853320467274178032973640562d0];  end;
if firstCall,   first=[true];  end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT  DBINOM
if( first )
bilnmx = log(d1mach(2)) - 0.0001d0;
fintmx = 0.9d0./d1mach(3);
end;
first = false;
%
if( n<0 || m<0 )
xermsg('SLATEC','DBINOM','N OR M LT ZERO',1,2);
end;
if( n<m )
xermsg('SLATEC','DBINOM','N LT M',2,2);
end;
%
k = fix(min(m,n-m));
if( k<=20 )
if( k.*log(max(n,1))<=bilnmx )
%
dbinomresult = 1.0d0;
if( k==0 )
csnil=dbstack(1); csnil=csnil(1).name(1)~='@';
if csnil&&~isempty(inputname(1)), assignin('caller','FUntemp',n); evalin('caller',[inputname(1),'=FUntemp;']); end
if csnil&&~isempty(inputname(2)), assignin('caller','FUntemp',m); evalin('caller',[inputname(2),'=FUntemp;']); end
return;
end;
for i = 1 : k;
xn = n - i + 1;
xk = i;
dbinomresult = dbinomresult.*(xn./xk);
end; i = fix(k+1);
%
if( dbinomresult<fintmx )
dbinomresult = fix(dbinomresult+0.5d0);
end;
csnil=dbstack(1); csnil=csnil(1).name(1)~='@';
if csnil&&~isempty(inputname(1)), assignin('caller','FUntemp',n); evalin('caller',[inputname(1),'=FUntemp;']); end
if csnil&&~isempty(inputname(2)), assignin('caller','FUntemp',m); evalin('caller',[inputname(2),'=FUntemp;']); end
return;
end;
end;
%
% IF K.LT.9, APPROX IS NOT VALID AND ANSWER IS CLOSE TO THE OVERFLOW LIM
if( k<9 )
xermsg('SLATEC','DBINOM','RESULT OVERFLOWS BECAUSE N AND/OR M TOO BIG',3,2);
end;
%
xn = n + 1;
xk = k + 1;
xnk = n - k + 1;
%
corr = d9lgmc(xn) - d9lgmc(xk) - d9lgmc(xnk);
dbinomresult = xk.*log(xnk./xk) - xn.*dlnrel(-(xk-1.0d0)./xn)- 0.5d0.*log(xn.*xnk./xk) + 1.0d0 - sq2pil + corr;
%
if( dbinomresult>bilnmx )
xermsg('SLATEC','DBINOM','RESULT OVERFLOWS BECAUSE N AND/OR M TOO BIG',3,2);
end;
%
dbinomresult = exp(dbinomresult);
if( dbinomresult<fintmx )
dbinomresult = fix(dbinomresult+0.5d0);
end;
%
csnil=dbstack(1); csnil=csnil(1).name(1)~='@';
if csnil&&~isempty(inputname(1)), assignin('caller','FUntemp',n); evalin('caller',[inputname(1),'=FUntemp;']); end
if csnil&&~isempty(inputname(2)), assignin('caller','FUntemp',m); evalin('caller',[inputname(2),'=FUntemp;']); end
end
%DECK DBINT4