Can someone explain this Miller-Rabin Primality test pseudo-code in simple terms?

Question

Here it is...

Input: n > 3, an odd integer to be tested for primality;
Input: k, a parameter that determines the accuracy of the test
Output: composite if n is composite, otherwise probably prime
Write n − 1 as (2^s)·d with d odd by factoring powers of 2 from n − 1
WitnessLoop: repeat k times:
   pick a random integer a in the range [2, n − 2]
   x ← a^d mod n
   if x = 1 or x = n − 1 then do next WitnessLoop
   repeat s − 1 times:
      x ← x^2 mod n
      if x = 1 then return composite
      if x = n − 1 then do next WitnessLoop
   return composite
return probably prime

I got this from the Wikipedia article on the Miller-Rabin primality test. But I've not been able to comprehend it...... I'm not looking to understand the math behind it but only to implement it in a program. This algorithm seems to me, to be kind of confusing. A better, more simpler pseudo-code or implementation of it in vb.net, would be helpful.

EDIT code written so far:

Function Miller_Rabin(ByVal n As Integer) As Boolean
    If n <= 3 Then : Return True
    ElseIf n Mod 2 = 0 Then : Return False
    Else
        Dim k, s, a, d, x As Integer
        k = 3
        d = n - 1

        While d Mod 2 = 0
            d = d / 2
            s += 1
        End While

        For c = 1 To k
            a = Random(2, n - 1)
            x = a ^ d Mod n
            If x = 1 Or x = n - 1 Then GoTo skip
            For r = 1 To s - 1
                x = x ^ 2 Mod n
                If x = 1 Then
                    Return False
                    Exit Function
                Else
                    If x = n - 1 Then
                        GoTo skip
                    Else
                        Return False
                        Exit Function
                    End If
                End If
            Next
skip:   Next
        Return True
    End If
End Function

Function Random(ByVal x As Integer, ByVal n As Integer) As Integer
    Dim a As Integer = Now.Millisecond * Now.Second
skip:
    a = (a ^ 2 + 1) Mod (n + 1)
    If a < x Then
        GoTo skip
    Else
        Return a
    End If
End Function

Answer 1

Here is simple pseudocode, as requested:

function isStrongPseudoprime(n, a)
    d := n - 1; s := 0
    while d % 2 == 0
        d := d / 2
        s := s + 1
    t := powerMod(a, d, n)
    if t == 1 return ProbablyPrime
    while s > 0
        if t == n - 1
            return ProbablyPrime
        t := (t * t) % n
        s := s - 1
    return Composite

function isPrime(n)
    for i from 1 to k
        a := randInt(2, n-1)
        if isStrongPseudoprime(n, a) == Composite
            return Composite
    return ProbablyPrime

function powerMod(b, e, m)
    x := 1
    while e > 0
        if e % 2 == 1
            x := (b * x) % m
        b := (b * b) % m
        e := e // 2 # integer division
    return x

The isStrongPseudoprime function tests if a is a witness to the compositeness of n; note that if isStrongPseudoprime returns Composite the number is definitely composite, but the opposite of that is ProbablyPrime because there is a chance that the number is still composite. The isPrime function tests k witnesses; by setting the value of k you can determine the likelihood of error as 1 chance in 4^k. Most people use a value of k somewhere between 10 and 25. The powerMod function performs exponentiation by squaring, and is provided on the chance that your language doesn't provide it for you.

If you want to know more about the mathematics behind this test, I modestly recommend this essay at my blog, which also includes implementations in five languages, though none of them is VBA.

EDIT: Although he didn't say so, what the original poster actually wants to do is to find the sum of the primes less than two million and thus solve Project Euler 10. Looping through the numbers from 2 to n is a very inefficient way to sum the primes less than n; instead, the recommended method is to use a sieve. Pseudocode again:

function sumPrimes(n)
    sum := 0
    sieve := makeArray(2..n, True)
    for p from 2 to n step 1
        if sieve[p]
            sum := sum + p
            for i from p * p to n step p
                sieve[i] := False
    return sum

The algorithm used here is the Sieve of Eratosthenes, invented over two thousand years ago by a Greek mathematician. Again, an explanation and code is in the essay at my blog.

Answer 2

Key Ideas and Concepts (p stands for prime here):

1. Fermat's Little Theorem. ( a^(p-1) = 1 ( mod p ))
2. If p is prime and x^2 = 1 ( mod p ), then x = +1 or -1 ( mod p ).

We can prove this as follows:

x^2 = 1 ( mod p )
x^2 - 1 = 0 ( mod p )
(x-1)(x+1) = 0 ( mod p )

Now if p does not divide both (x-1) and (x+1) and it divides their product, then it cannot be a prime, which is a contradiction. Hence, p will either divide (x-1) or it will divide (x+1), so x = +1 or -1 ( mod p ).

Let us assume that p - 1 = 2^d * s where s is odd and d >= 0. If p is prime, then either as = 1 ( mod p ) as in this case, repeated squaring from as will always yield 1, so (a^(p-1))%p will be 1; or a^(s*(2^r)) = -1 ( mod p ) for some r such that 0 <= r < d, as repeated squaring from it will always yield 1 and finally a^(p-1) = 1 ( mod p ). If none of these hold true, a^(p-1) will not be 1 for any prime number a ( otherwise there will be a contradiction with fact #2 ).

Algorithm :

1. Let p be the given number which we have to test for primality.
2. First we rewrite p-1 as (2^d)*s. (where s is odd and d >= 0).
3. Now we pick some a in range [1,n-1] and then check whether as = 1 ( mod p ) or a^(s*(2^r)) = -1 ( mod p ).
4. If both of them fail, then p is definitely composite. Otherwise p is probably prime. We can choose another a and repeat the same test.
5. We can stop after some fixed number of iterations and claim that either p is definitely composite, or it is probably prime.

Small code : Miller-Rabin primality test, iteration signifies the accuracy of the test

bool Miller(long long p,int iteration)
{
    if(p<2)
        return false;

    if(p!=2 && p%2==0){
                return false;

        long long s=p-1;
        while(s%2==0)
        {
                s/=2;
        }
        for(int i=0;i<iteration;i++)
        {
                long long a=rand()%(p-1)+1;
            long long temp=s;
                long long mod=modulo(a,temp,p);
                while(temp!=p-1 && mod!=1 && mod!=p-1)
            {
                    mod=mulmod(mod,mod,p);
                    temp *= 2;
                }
                if(mod!=p-1 && temp%2==0)
             {
                    return false;
                }
         }
         return true;
}

Few points about perfomance:

It can be shown that for any composite number p, at least (3/4) of the numbers less than p will witness p to be composite when chosen as 'a' in the above test. Which means that if we do 1 iteration, probability that a composite number is returned as prime is (1/4). With k iterations the probability of test failing is (1/4)k or 4(-k). This test is comparatively slower compared to Fermat's test but it doesn't break down for any specific composite numbers and 18-20 iterations is a quite good choice for most applications.

PS: This function calculates (a*b)%c taking into account that a*b might overflow WHICH I HAVE USED ABOVE IN MILLER RABBIN TEST.

   long long mulmod(long long a,long long b,long long c)
   {
       long long x = 0,y=a%c;
       while(b > 0)
       {
          if(b%2 == 1)
          {
              x = (x+y)%c;
          }
          y = (y*2)%c;
          b /= 2;
       }
       return x%c;
    }

Answer 3

The VB implementation uses a hexadecimal conversion function to handle large numbers before the modular exponentiation. Example provided in comments:

' USAGE:
' Example: strResult = mpModExp("3c", "03", "face")
' computes (0x3c)^3 mod 0xface = 0x5b56
' or, in decimal, 60^3 mod 64206 = 23382
' Parameters may be hex strings of any length subject to limitations
' of VB and your computer. May take a long time!