[binomresult,n,m]=binom(n,m);

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

if isempty(bilnmx), bilnmx=0; end;
if isempty(binomresult), binomresult=0; end;
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(i), i=0; end;
if isempty(k), k=0; end;
%***BEGIN PROLOGUE  BINOM
%***PURPOSE  Compute the binomial coefficients.
%***LIBRARY   SLATEC (FNLIB)
%***CATEGORY  C1
%***TYPE      SINGLE PRECISION (BINOM-S, DBINOM-D)
%***KEYWORDS  BINOMIAL COEFFICIENTS, FNLIB, SPECIAL FUNCTIONS
%***AUTHOR  Fullerton, W., (LANL)
%***DESCRIPTION
%
% BINOM(N,M) calculates the binomial coefficient (N!)/((M!)*(N-M)!).
%
%***REFERENCES  (NONE)
%***ROUTINES CALLED  ALNREL, R1MACH, R9LGMC, XERMSG
%***REVISION HISTORY  (YYMMDD)
%   770701  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)
%   900326  Removed duplicate information from DESCRIPTION section.
%           (WRB)
%***end PROLOGUE  BINOM
if isempty(first), first=false; end;
if firstCall,   sq2pil=[0.91893853320467274e0];  end;
if firstCall,   first=[true];  end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT  BINOM
if( first )
bilnmx = log(r1mach(2));
fintmx = 0.9./r1mach(3);
end;
first = false;
%
if( n<0 || m<0 )
xermsg('SLATEC','BINOM','N OR M LT ZERO',1,2);
end;
if( n<m )
xermsg('SLATEC','BINOM','N LT M',2,2);
end;
%
k = fix(min(m,n-m));
if( k<=20 )
if( k.*log(max(n,1))<=bilnmx )
%
binomresult = 1.;
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;
binomresult = binomresult.*real(n-i+1)./i;
end; i = fix(k+1);
%
if( binomresult<fintmx )
binomresult = fix(binomresult+0.5);
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','BINOM','RESULT OVERFLOWS BECAUSE N AND/OR M TOO BIG',3,2);
end;
%
xn = n + 1;
xk = k + 1;
xnk = n - k + 1;
%
corr = r9lgmc(xn) - r9lgmc(xk) - r9lgmc(xnk);
binomresult = xk.*log(xnk./xk) - xn.*alnrel(-(xk-1.)./xn)- 0.5.*log(xn.*xnk./xk) + 1.0 - sq2pil + corr;
%
if( binomresult>bilnmx )
xermsg('SLATEC','BINOM','RESULT OVERFLOWS BECAUSE N AND/OR M TOO BIG',3,2);
end;
%
binomresult = exp(binomresult);
if( binomresult<fintmx )
binomresult = fix(binomresult+0.5);
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 BINT4