Special Functions math library: fact(n) - 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.

fact(n)

function [f] = fact(n)
%FACT  Vectorized Factorial function
%
%usage: f = fact(n)
%
%tested under version 5.3.1
%
%     This function computes the factorial of 
%     the elements of N.
%     N can be any size but must contain
%     Real, Non-Negative, Integers.
%
%     This routine is much more robust than
%     the built in FACTORIAL function.
%
%see also: Gamma, Prod, Factorial, Binomial

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

[row,col]=size(n);
n=n(:);

f=NaN*n;

pp=find(imag(n)==0 & real(n)>=0 & round(n)==n);
%find integer values
nn=n(pp);
ff=zeros(length(pp),1);

s=1;
p=[];
if ~isempty(nn)
   p=find(nn==0);
end
if ~isempty(p)
   ff(p)=s;
end

%upper limit here depends upon realmax
if ~isempty(nn)
for k=1:170
    s=s*k;
%empty=scalar warning
    p=find(nn==k);
    if ~isempty(p)
        ff(p)=s;
    end
end
end

p=find(nn>170);
if ~isempty(p)
    ff(p)=Inf;
end

f(pp)=ff;
f=reshape(f,row,col);

return

Highlights from Special Functions math library

Highlights from
Special Functions math library