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')
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