Special Functions math library: binomial(n,d) - File Exchange

No BSD License

Highlights from
Special Functions math library

bern(n) Bern Bernoulli number
betad(z) BETAD Dirichlet Beta function
binomial(n,d) BINOMIAL calculate the binomial coefficient
deta(z,k) DETA Calculates Dirichlet functions of the form
erfz(zz) ERFZ Error function for complex inputs
eta(z) ETA Dirichlet Eta function
euler(n) Euler Euler number
eulergamma Euler-Mascheroni constant = -Psi(1) = 0.5772156649015328606...
fact(n) FACT Vectorized Factorial function
factd(n) FACTD Double Factorial function = n!!
gamma(z) GAMMA Gamma function valid in the entire complex plane.
gammaln(z) GAMMALOG Natural Log of the Gamma function valid in the entire complex plane.
genocchi(z) Genocchi number
harm(z) Harm Harmonic sum function valid in the entire (complex) plane.
lambda(z) LAMBDA Dirichlet Lambda function
poch(z,n)
psi(z) Psi Psi (or Digamma) function valid in the entire complex plane.
psin(n,z) Psin Arbitrary order Polygamma function valid in the entire complex plane.
totient(n) TOTIENT calculates the totient function (also
zeta(z) ZETA Riemann Zeta function
View all files

from Special Functions math library by Paul Godfrey
Collection of Special Functions programs.

binomial(n,d)

function [f] = binomial(n,d)
%BINOMIAL  calculate the binomial coefficient
%          n and d may be complex and any equal size
%
%        n!
% f = --------  
%     d!(n-d)!
%
%usage: b = binomial(n,d)
%
%tested on version 5.3.1
%
%see also: Gamma, Prod
%see also: mhelp binomial
%
%not correct for negative integer arguments!

%Paul Godfrey
%pgodfrey@conexant.com
%8-12-00

nsize=size(n); n=n(:);
dsize=size(d); d=d(:);

if min(nsize==dsize)==0
   error('Input argument size mismatch!')
end

if (0<n & 0<=d & round(n)==n & round(d)==d & n>=d & max(size(n))==1 & max(size(d))==1)
   f = prod((d+1):n)/prod(1:(n-d));
else
   z=gamma([n d n-d]+1);
   f=z(:,1)./(z(:,2).*z(:,3));
end

f=reshape(f,nsize);

return

%a demo of this routine is
dx=1/16;
dy=dx;
x=-4:dx:4;
y=-4:dy:4;
[X,Y]=meshgrid(x,y);
f=binomial(X,Y);
p=find(abs(f)>4);
%f(p)=sign(real(f(p)))*4;
f(p)=NaN;
mesh(x,y,real(f))
view([23 54]);
rotate3d

Highlights from Special Functions math library

Highlights from
Special Functions math library