This program should be saved in two parts.  The top
part is the main script and may be saved under
any name you wish, provided that it has a .m
suffix.  The second part should be saved under
the name nfactors.m


___________________________


clear
echo on
%This script file uses the function nfactors.m
%to plot up the density of prime factors for
%numbers up to 10000.  Note that the file
%nfactors.m must be in the same directory as
%this script so that matlab can find it.

xmax=10000;
x=[1:xmax];
%This command generates an array of integers
%ranging from 1 to 10000.
pause

y=nfactors(x);
%This command calls the function nfactors to
%generate the number of prime factors for
%each element in x.  The answer is placed in
%the array y.
pause

figure(1)
plot(x,y)
xlabel('N')
ylabel('number of prime factors')
%This command gives a line plot between x and
%y.  Note that you really can't tell a whole
%lot about the density from this plot.
pause

figure(2)
semilogx(x,y,'o')
xlabel('N')
ylabel('number of prime factors')
%This command produces a semilog plot of x vs y.
%Note that from this type of plot we can see
%much more of the structure of the prime factor
%density.

pause

%We can use this to learn more about prime numbers
%as well using the "find" command:

i=find(y==1);
primes=x(i);
%which are all the prime numbers in the array x.

%We can use the histogram command to look at the
%distribution of prime numbers as well:
figure(3)
nbins=20
hist(primes,nbins)
xlabel('Integer Range')
ylabel('Number of Primes')
hold on
plot(xmax/nbins*[.5:nbins-.5],xmax/nbins./log(xmax*[.5:nbins-.5]/nbins),'r')
legend('number density','dx/log(x) asymptote')
hold off


%The study of the distribution of prime numbers has
%a very long history.  In fact, the Riemann Hypothesis
%(related to this distribution) is still unsolved, and
%if you could prove it you would collect the $1M 
%Millenium Prize from the Clay Mathematical Institute...


echo off


__________________________


function n=nfactors(m)
%This function calculates the number of prime
%factors n of the integers in the array m.  This
%program demonstrates the use of while loops.

np=length(m);
%This command gets the number of integers in the
%array m that we will operate on.

n=ones(size(m));
%This command initializes the output array.

for i=1:np
%We work on the elements one at a time.

mt=m(i);
k=2;
if (mt>1)&(fix(mt)==mt)
%The command & is the logical and.  Fix rounds
%to the integer closest to zero.

  while k<=mt^.5
    while (fix(mt/k)==mt/k)&(k< mt)
          n(i)=n(i)+1;
          mt=mt/k;
        end;
        k=k+1;
  end;
end;

end;

%The program returns the integer array n containing
%the number of prime factors.