% Scripts and functions for MAT 331,  Sept 21, 2017

% function to find nth prime number 

  function our = nth_prime(n)
    m = round( 4 * n * log(n));
    x=find(isprime([1:m]));
    out = x(n);
  return;

%---------------------------------
% function to find nth prime number 

  function our = nth_prime(n)
    m = round( 4 * n * log(n));
    x=find(isprime([1:m]));
    out = sum(x(1:n));
  return;

%---------------------------------
% script to  count number of primes less than n

   x=[1:n];
   sum((isprime(x));

%------------------------------------
% script to sum primes less than n

   x=[1:n];
   sum(find(isprime(x)));

%----------------------------------
% script to count palindromic primes 

counter=0;
N=100000;
clear save_primes; % use this to save the primes found
for k=1:N
  if isprime(k)
     a=num2str(k);
     b=a([length(a):-1:1]);
     %  c=str2num(b)  % use these to find primes whose 'reverse'
     %  if isprime(c) % is also prime
     if  a==b
       counter=counter+1;
       a;
       save_primes(counter)=k;
     end
  end
end % for k
counter

%--------------------------------
% script to count twin primes 

   L=100000;
   x=isprime(1:L);
   y=x(3:end);
   z=x(1:end-2);
   sum(y.*z)

%-----------------------------------
%  script to count digits in prime numbers

   counter=zeros(10,1);
   N=1000000;
   x=find(isprime([1:N])); % create list of primes <= N
   for k=1:length(x)
      a=num2str(x(k));     %convert to string
      for j=1:length(a)
         n=str2num(a(j))+1;     % add to counter array, starts with 1
         counter(n)=counter(n)+1;
      end
   end % for k
   counter
   figure; hold on; grid on;
   title(['Histogram of digits in primes less than ',num2str(N)]);
   bar([0:9],counter)

%-------------------------------------
% script to count number of primes = 1,3 mod 4

  N=10000000;
  count1=zeros(N,1);
  count3=zeros(N,1);
  for k=3:N
      count1(k)=count1(k-1);
      count3(k)=count3(k-1);
        if isprime(k)
           if mod(k,4)==1
              count1(k)=count1(k)+1;
           else
              count3(k)=count3(k)+1;
           end 
        end 
  end 
  figure; hold on; grid on;
  title('Number of primes equal to 1 or 3 modulo 4');
  plot(count1);
  plot(count3);
  legend('1 mod 4','3 mod 4');

  figure; hold on; grid on;
  title('Difference of number of primes equal to 3 or 1 mod 4');
  plot(count3-count1);

  find(count1>count3)