Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

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

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
like image 302
isados Avatar asked Oct 25 '25 14:10

isados


1 Answers

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.

like image 63
user448810 Avatar answered Oct 28 '25 03:10

user448810