function  MillerRabinTestEnhanced(n, k)
% y = MillerRabinTestEnhanced(n, k)
% Inputs: an odd integer n > 3, and a positive integer k:  Program will apply the 
% Miller-Rabin Enhanced Primality Test (Algorithm 8.4).  
% The second input variable k is optional
% (default value is 1).    The program will indicate that n was found to be
% composite, if at least one of the k tests has found n to be composite, 
% and that n is probably prime, in case all of the tests were inconclusive
% (meaning, informally, they all found n to be probably prime).   In cases 
% that n is found composite, the witness a, along with the corresponding 
% exponent parameters [j m], that resulted in the composite conclusion will also
%be outputted, or a nontrivial prime factor, in case one is detected.
%CAUTION:  This program does not have symbolic functionality, the analogous
%program that has is is named as MillerRabinTestEnhancedSym

if nargin < 2, k = 1; end
%first we determine the largest power f of 2 that divides n - 1:
f = 0;
while floor((n-1)/2^(f+1)) == (n-1)/2^(f+1)
    f = f + 1;
end
m = (n-1)/2^f;

for Test = 1:k
    flag = 0; %will be reset to 1 if Test finds n probably prime
    a = ceil((n - 3)*rand) + 1;%randomly generate a in range 1 < a < n - 1
    A = FastExp(a,m,n);
    if A == 1 | A == n - 1
        %n is declared probably prime; for this particular test,
        %go on to new test.
        flag = 1;
    else
        for j = 1:f-2
            newA = mod(A^2,n);
            if newA == 1
                %n is proved composite; output nontrivial factor, and exit program
                fprintf('Inputted integer was proved composite using the Miller-Rabin test.\r')
                fprintf('Here is a corresponding nontrivial factor: \r')
                gcd(A,n)
                return
            elseif newA == n-1
                %n is declared probably prime; for this particular test,
                %go on to new test.
                flag = 1;
                break  %out of j for loop
            end
            oldA = A; A = newA;
        end
        if flag == 0 %n has not yet been declared probably prime (or composite)
           A = mod(A^2,n); %this is a^(n-1)/2
           if A ~= n-1 & A ~= 1
              fprintf('Inputted integer was proved composite using the Miller-Rabin test.\r')
              fprintf('Here is a corresponding witness: \r')
              a
              fprintf('The corresponding exponent paramters for 2^j*m are [j m]=:\r')
              [k-1 m]
              return
              %in case A == 1, can get a nontrivial factor as above:
           elseif A == 1
                fac1 = gcd(n, oldA);
                fprintf('A nontrivial factor of the inputted integer has been found to be:\r')
                fac1
                return
           else
              flag = 1; %don't really need this
              %n is declared probably prime for this particular test,
              %go on to new test.
           end
        end
    end
end    
%If the program has not yet exited, this means that all k tests resulted in
%a probable prime conclusion, so this will be the overall conclusion.
fprintf('Inputted integer n is declared probably prime after %d applications\r', k)
fprintf('of the Miller-Rabin test.\r')
