function FermatTestSym(n, k)
% y = FermatTestSym(n, k)
% Input: an integer n > 3, and an optional second input k, a positive integer.
% Default value for k:  1
% The program will apply the Fermat primality test (Algorithm 8.2) k times.  
% Output: y will be 1 if at least one of the k tests has found n to be 
% composite, and 0 in case all of the tests were inconclusive (meaning, 
% informally, they all found n to be probably prime).
% Program accepts symbolic inputs and calls on the FastExpSym program.
y = 0;  %At the end, y will be 1 if at least one of the k tests has found n to be 
% composite, and 0 in case all of the tests were inconclusive (meaning, 
% informally, they all found n to be probably prime).
if nargin < 2, k = 1; end
for test = 1:k
    a = ceil((n-3)*rand)+1; %randomly select a base a:  1<a<n-1
    pwr = FastExpSym(a,n-1,n);
    if pwr ~=1  %proof that n is composite
        y = 1; 
        fprintf('Inputted number was proved composite using the Fermat text.\r)
        fprintf('Here is a corresponding witness: \r')
        a
        return
    end
end
if y ==0, 
    fprintf('Inputted integer n is declared probably prime after %d applications\r', k)
    fprintf('of the Fermat test.\r')
end

function c = FastExpSym(a,b,n)
% c = FastExpSym(a,b,n)
% Modification of my ModInv function (using my gcdSym function) that
% has symbolic functionality.
% Input:  a, b, < n positive integers, n = the modulus
% Output: c = a^b (mod n) 
% Program uses the fast exponentiation method to calculate c = a^b (mod n).  
% Beginning with factor = a^1: 
% on each pass, it is determined whether this is a desired factor.  If so, then 
% it is multiplied into a variable containing a cumulative product of desired 
% factors.  Then the factor is squared, and the process repeated.

exp = b;
c = ones(size(a));
factor = a;

while(exp ~= 0)
    if (mod(exp,2) == 1)
        c = mod(c .* factor,n);
    end
    
    exp = floor(exp/2);
    factor = mod(factor.*factor,n);
end
