Miller-Rabin prime test (Mathematica)

This page shows three Mathematica functions which together implement the Miller-Rabin prime test as described in your textbook, page 140.

The first one, bc, takes a number n as argument and returns the exponent b and odd factor c such that n == 2^bc. The underline has to be present on the left, but never on the right. Mathematica has some very sophisticated constructs, and this is one of them. You could define a function f[{a_,b_}], which would take one argument which is a list of two elements - but you could refer to the components a and b on the right hand side of the definition.

A Module is a construct that lets you define local variables. They are defined as a list in the first argument, and whatever the module computes is the second argument.

Sequential actions are indicated by statement terminators (semicolons). The last thing computed in the module is the list {b,c}, and that will be returned.

The While function has two arguments. The first is the condition, and the second is the action.

If you want integer division, you have to use the Quotient function.

bc[n_] := Module[ {b = 0, c = n},
While[ Mod[c,2] == 0, b++; c = Quotient[c,2] ];
{b, c}
]

The next function tries step 3 of your textbook algorithm. If n is certifiably not prime (cnp) it returns True, otherwise False.

The syntax for If is If[condition, then-action, else-action].

The algorithm is

Initialize the power to a and the product to 1.
While the exponent is nonzero
{
	If the low-order bit of the exponent c is 1, multiply the power into
            the product.
	Square the power and shift off the low-order bit of the exponent.
        If the squaring produced -1 (mod n), return False.
	If it produced 1, then we have a square root of 1 other than 1 or -1,
            so return True.
}
Now square the product b-1 times. If any square is -1 (mod n), return False,
    but if it's 1, return True.
At the end, if we have either 1 or -1 (mod n), then return(False), otherwise
    return True.

And here's the Mathematica code:

cnp[n_, a_, b_, c_] := Module[{ex = c, pr = 1, po = a, twos = b},
While[ex > 0,
   If[ Mod[ex, 2] == 1, pr = Mod[pr * po, n], ];
   po = Mod[po * po, n];
   ex = Quotient[ex, 2];
   If[ po == n-1, Return[False],];
   If[ po == 1, Return[True], ];
];
If[ pr == n - 1, Return[False], ];
If[pr == 1, Return[False], ];
While[twos > 1,
   pr = Mod[pr*pr, n];
   twos--;
   If[ pr == n - 1, Return[False],];
   If[pr == 1, Return[True],];
];
If[ (pr == 1) || pr == n-1,Return[False],Return[True]];

]

Now probablyprime is easy. Specify the number n to test and the number of trials to do.


probablyprime[n_,ntrials_] := Module[
      {bcres = bc[n - 1]},             
For[i = 0, i < ntrials, i++,
    If[cnp[n, Random[Integer,{2, n - 2}],First[bcres],Last[bcres]],Return[False],]
];
Return[True]
]

Now if you input to Mathematica
probablyprime[1111111111111111111, 20]
there is a noticeable delay, and then out comes
True
The number really is prime.

With a little modification we can make the function show us what it's doing:


probablyprime1[n_,ntrials_] := Module[
      {bcres = bc[n - 1], a},             
For[i = 0, i < ntrials, i++,
    a = Random[Integer, {2, n-2}];
    Print["Trying ", a];
    If[cnp[n, a ,First[bcres],Last[bcres]],Return[False],]
];
Return[True]
]

with results like this:


probablyprime1[1111111111111111111, 20]
Trying 815863354527691421
Trying 52387242841935937
Trying 211783923785200936
Trying 834816113346351625
Trying 1101383424989427154
Trying 792994771316060541
Trying 1098213212338499285
Trying 248230583861041812
Trying 186850293519980942
Trying 850654053640527681
Trying 666888673440624253
Trying 197099378035149258
Trying 151919387204106385
Trying 514500642065846973
Trying 504152880762141267
Trying 492512123277809370
Trying 712823375202928580
Trying 399155884802181009
Trying 647765321150763977
Trying 723569789019027366
True