function p = PrimeTheorem(N)
close all
% Calculate and remember all PRIME NUMBERS <= N
% Ex. (assuming N=15): p = [2,3,5,7,11,13]
p = primes(N);
np = length(p); % number of primes <= N
% Count the NUMBER OF PRIMES <= n for all values of n up to N 
% Ex. (assuming N=15): piN = [1,2,2,3,3,4,4,4,4,5,5,6,6,6]
n=2:N;
piN=zeros(1,N-1);
for i=n
    piN(i-1)=sum(p<=i);
end

% Plot differences of consecutive pairs of primes.
% Note how larger differences generally don't occur until later in the sequence.
plot(diff(p),'.')
title(sprintf('(Consecutive) differences for the np=%d primes up to N=%d.',np,N))

% Plot exact density of primes and Gauss's 1/log(n) approximation
figure
plot(n,piN./n,'g',n,1./log(n),'r')
legend('exact density \pi(n)/n','Gauss''s approximate density 1/log(n)')
title(sprintf('Exact and approximate density for the primes up to N=%d.',N))

% Plot number of primes up to N along Gauss's estimate (based on his density), 
% i.e., pi(N) \approx 1/log(2) + 1/log(3) + ... + 1/log(N)
figure
plot(n,piN,'g',n,cumsum(1./log(n)),'r')
legend('\pi(n)','1/log(2) + 1/log(3) + ... + 1/log(n)','Location','NorthWest')
title(sprintf('Number of primes up to N=%d, and Gauss''s estimate.',N))

% Plot primes together with their n*log(n) approximation 
figure
n=1:np;
plot(n,p,'g',n,n.*log(n),'r')
legend('primes p_n','approximation n*log(n)','Location','NorthWest')
title(sprintf('The np=%d primes p_n up to N=%d and their n*log(n) approximation.',np,N))
