echo on
%This script uses the function logfac(n)
%to calculate the log of n! and to compare
%the exact result to two asymptotic formulas.
%
%First we set up the range in n:
n=1:20;
%And then we plot:
plot(n,logfac(n),'o')

%We compare this graphically to two approximate
%formulas:
hold on
stir=n.*log(n)-n;
plot(n,stir)
asymp=(n+0.5).*log(n)-n+0.5*log(2*pi);
plot(n,asymp,'--')
%And now we do a bit of labelling:
ylabel('log(n!)')
xlabel('n')
legend('exact','stirling','asymptotic')
hold off

pause

%It is actually possible to do the entire
%problem in a much more compact way using
%specialized matlab commands.  We can
%make use of the command 'cumsum' to take
%the cumulative sum of the elements of
%a vector:

plot(n,cumsum(log(n)),'o',n,stir,n,asymp,'--')
ylabel('log(n!)')
xlabel('n')
legend('exact','stirling','asymptotic')

%which eliminates the need for the subprogram
%entirely!  Of course, part of the reason I
%asked you to write a subprogram (and didn't
%tell you about cumprod or cumsum) is that
%you need practice in both writing functions
%and doing loops.

echo off



______________________________

function y=logfac(n)
%This program computes the log of n! where
%n is permitted to be a vector.  It will do
%it the fancy way, so that overflow can be
%avoided even for large values of n.  You did
%not have to get this fancy on your answer.

y=zeros(size(n));  %We initialize the output array.
npt=length(n);
b=10^10;
for i=1:npt
 nt=n(i);
 y(i)=1;
 m=0;
 for j=1:nt
  y(i)=y(i)*j;
  if y(i) > b
   y(i)=y(i)/b;
   m=m+1;
  end;
 end;
 y(i)=log(y(i))+m*10*log(10);
end;