function TF = isPrimePocklington(N,pockiter)
% isPrimePocklington: Pocklington certificate of primality
% usage: TF = isPrimePocklington(N,pockiter)
%
% Produces a Pocklington certificate of primality, which applies only when
% the complete factorization of N-1 is available.
%
%  N - large integer to be tested for primality
%
%  pockiter - an integer at least 1, which indicates the number of
%      random witnesses to be chosen for the certificate.
%      pockiter must be no larger than 1228. This is because the pool of
%      potential witnesses comes from setdiff(primes(10000),2).
% 
%      default: pockiter = 10
%
% Arguments: (output)
%  TF - a result with 3 possible values
%      TF == true  --> N is certainly prime
%      TF == NaN   --> N is probably prime, but not certainly so
%      TF == false --> N is certainly composite

% How many iterations will we try? 10 is the default.
if (nargin < 2) || isempty(pockiter)
  pockiter = 10;
end

% The list of potential witnesses will come from the primes not exceeding
% 10000. Primes are chosen here, because the best witnesses will always be
% square-free. 2 is also excluded, as 2 is the one witness that seems to be
% more commonly a Fermat liar in my experience.
potentialwitnesses = primes(10000);
potentialwitnesses(1) = [];
wlist = potentialwitnesses(randperm(1228,pockiter));

% We must be able to easily factor N-1. We only need the distinct
% factors of N-1.
q = unique(factor(N-1));

TF = NaN;
for iter = 1:pockiter
  w = wlist(iter);

  % perform a Fermat test
  if N == w
    % since all witnesses were prime, then N must be prime
    TF == true;
    return
  elseif gcd(w,N) ~= 1
    % N must share a divisor with w, so is composite
    TF = false;
    return
  end

  if powermod(w,N-1,N) ~= 1
    % N fails a Fermat test for primality. This insures N is composite
    TF = false;
    return
  end

  % Test the pocklington criterion, for EACH distinct factor q of N-1, that
  g = gcd(powermod(w,(N-1)./q-1,N),N);
  if all(g) == 1
    % if we drop into here, then N is certainly prime.
    TF = true;
    return
  end
  % if we fail the above test, then it is not certain that N is composite.
  % In fact, N is very likely prime, as it must have satisfied the Fermat
  % test. So repeat the above sequence with another witness from the set.

end