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
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.
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With