Miller-Rabin primality test

General FreeBASIC programming questions.
frisian
Posts: 236
Joined: Oct 08, 2009 17:25

Miller-Rabin primality test

Postby frisian » Nov 15, 2012 16:57

The program is a translation of the PureBasic program found on Rosetta Code page.
When you get wrong results you must set the numbers of loops higher, say 10 instead of 5.

http://rosettacode.org/wiki/Miller-Rabin_primality_test#Run_BASIC

Code: Select all

#Define prob_prime -1
#Define composite 0

Function pow_mod(b As ULongInt, power As UInteger, modulus As UInteger) As UInteger

  Dim As ULongInt x = 1

  While power > 0
    If (power And 1) = 1 Then
      x = (x * b) Mod modulus
    EndIf
    b = (b * b) Mod modulus
    power = power Shr 1
  Wend
  Return x

End Function

Function prob_prime_test(n As ULongInt, k As Integer) As Byte

  If n > 4294967295 Then                ' limit, pow_mod can't handle bigger numbers
    Print "number is to big, program will end"
    Sleep
    End
  EndIf
  ' 2 is a prime, if n is smaller then 2 or n is even then n = composite
  If n = 2 Then Return prob_prime
  If (n < 2) Or ((n And 1) = 0) Then Return composite

  Dim As ULongInt x, n_one = n - 1, d = n_one
  Dim As UInteger s, r, a

  While (d And 1) = 0
    d = d \ 2
    s = s + 1
  Wend

  While k > 0
    k = k - 1
    a = Int(Rnd * (n - 2)) + 2          ' 2 <= a < n
    x = pow_mod(a, d, n)
    If (x = 1) Or (x = n_one) Then Continue While
    For r As Integer = 1 To s-1
      x = (x*x) Mod n                   ' x = pow_mod(x, 2, n)
      If x = 1 Then Return composite
      If x = n_one Then Continue While
    Next
    If x <> n_one Then Return composite
  Wend
  Return prob_prime

End Function

' -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
' MAIN
' -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

Randomize timer

Dim As UInteger total, limit = 100000000

For y As UInteger = 1 To limit
  If prob_prime_test(y,5) = prob_prime Then
    total = total + 1
  EndIf
Next

Print
Print total;" primes found below ";limit+1


Sleep
End
integer
Posts: 375
Joined: Feb 01, 2007 16:54
Location: usa

Re: Miller-Rabin primality test

Postby integer » Nov 22, 2012 22:00

From a earlier post:
counting_pine wrote:
AlejandroPadrino wrote:Thank you, Counting_pine. I was seeing some prime tests at Wikipedia, but all they results are for "probably prime numbers".


Miller-Rabin has deterministic variants for certain limits - i.e. you can run a set of specific tests, and it will guarantee primality for a specific set of numbers, as verified by others. (e.g. it's easy to check that a witness of 2 is reliable for any number up to 2047.)

For example, the page for Miller-Rabin shows three tests that you can run which guarantee primality for any number up to 4.7 billion, i.e. any 32-bit integer. Additionally, http://miller-rabin.appspot.com shows a set of 7 tests that reportedly works for any number up to 2^64.


Using those 7 base numbers with a slight modification of frisian's code and richard's Biginteger functions the following will give a
deterministic result about the prime status of any specific numbers < 2^64.

There are some lines of code (part of a larger test) that can be removed, but it's Thanksgiving, so it stays in.

Code: Select all

'' MillerRabinDeterministic.bas

#include once "Big_Integer_2.bas"   '' The overloaded version is super!

#define AlternateTest 0     
'       AlternateTest 0     '' default value; use the 7 base number test n<2^64 DETERMINISTIC
'       AlternateTest 1     '' use 1st 12 primes as base numbers to test n<2^64 DETERMINISTIC
'       AlternateTest 2     '' use 1st 127 primes as base numbers to test n<2^256 Probable
'       AlternateTest 3     '' use 200 Random Numbers as base numbers to test n<2^2048 Probable

DECLARE FUNCTION MillerRabin64(byref anumber as bigint ) as byte        'as close to boolean as possible here
DECLARE FUNCTION modpow(byref basenum as bigint, byref exponentnum as bigint, byref modulusnum as bigint ) as bigint
DECLARE SUB GeneratePrimeNumberList()


    '' using these numbers the MR test is deterministic for numbers less than 2^64
    dim shared as bigint BaseNumber(1 to 7) = { 2, 325, 9375, 28178, 450775, 9780504, 1795265022}   
   
    '' use first 12 primes as an alternate test (also deterministic for n < 2^64 )
    dim shared as bigint PrimeNumber()
    redim shared as bigint PrimeNumber( 1 to 128 )  '={2,3,5,7,11,13,...}

'------------------------------------------
    if AlternateTest > 2 then stop : end if     ' this test excluded here
'------------------------------------------

    if AlternateTest then
        print "generating primenumber list"
        GeneratePrimeNumberList()
        print "primes list generated."
    end if
   
    dim as BigInt n
    dim as byte t
    dim as integer count = 0
    dim as integer primecount = 0, compositecount = 0
    dim as double t1, t0
   
    ' start testing somewhere
    n =  &hffffffffffffffffull      ' Largest Unsigned Long Integer is ODD

    print "  starting with : ";n
    t0 = timer
    do                              ' a few test numbers
        t = MillerRabin64(n)
        if t = 1 then primecount += 1 : end if      ' I'm sure it's one prime number
        if t = 0 then compositecount += 1 : end if  ' nada prime
        n -= bigint_2               ' decrement to the next odd number
        count += 1
    loop until ( count >= 100 )
    t1 = timer
    print "    ending with : ";n
    print "odd #'s only"
    print "    PRIMES found: ";primecount
    print "composites found: ";compositecount
    print " seconds elapsed: ";t1-t0
    sleep
    end
   
   
' results for 1000 oddnumbers less than maximum ULONGINT
' approximately  8 seconds using 7 base numbers
' approximately 10 seconds using 12 primes as base numbers
' approximately 45 seconds using 127 primes as base numbers


' frisian's MR test for Big Integers
'_______
FUNCTION MillerRabin64(byref anumber as bigint ) as byte
   
    dim as bigint n = anumber       ' protect the source data
    dim as bigint nm1               ' N Minus 1
    dim as bigint d                 ' as in d*2^s
    dim as bigint answer            ' temp/intermediate variable
    dim as byte ans                 ' as close to boolean
    dim as integer s                ' could use a short

    'exclude negatives, 0, 1 & even numbers
    if n = bigint_2 then return(1):end if     ' if it's a 2 then it's PRIME
    if n < bigint_2 then return(0):end if     ' ignore NEGATIVE numbers & 1
    if remainder(n,bigint_2) = bigint_0 then return(0):end if '  & even #'s
    '' n is an oddnumber greater than 2
   
    'writing n-1 as d*2^s (with d being odd)   
    Nm1 = N - bigint_1              ' = n-1
    S = 0
    D = Nm1                         ' initial condition: d = n - 1 {even}
    'while (D is even):D\=2:S+=1:wend   
    dim as bigint quotient, residue   
    do
        divmod(d, 2, quotient, residue )
        if residue = bigint_0 then  ' d is even
            d = quotient            ' d = d\2
            s += 1                  ' exponent count; 2^s ;  n-1 = d * 2^s
        else
            exit do                 ' d is odd, job done.
        end if
    loop
       
       
    ' testing the seven bases: deterministic for n < 2^64
    ' as alternate use first 12 primes: 2,...,37} : deterministic for n<2^64
    ' or 127 primes  if testing in range 2^64 - 2^128
    ' or 200 random numbers if testing greater than 2^128   
   
    dim as integer TestLimit = 7   
/'
    if alternateTest = 0 then TestLimit = 7   '' the default value
'/   
    if alternateTest = 1 then TestLimit = 12  '' 1ST 12 PRIMES
    if alternateTest = 2 then TestLimit = 127 '' 128 = negative BYTE#
    'if alternateTest = 3 then TestLimit = 200 '' 200 random numbers
   
    dim as bigint b, x   
    for I as integer = 1 to TestLimit
       
        if alternatetest = 0 then b = BaseNumber(I) : end if
       
        if alternatetest = 1 or alternatetest = 2 then     ''which test to use
            b = primenumber(i)
        'elseif alternatetest = 3 then           
        '    dim as double TheNumber
        '    TheNumber = BigInt_2_Double ( n )
        '    b = Double_2_Bigint ( TheNumber * rnd)
        '    if b < 2 then b = 2 : end if
        '    if b > (n - bigint_2) then b = n - bigint_2 : end if
        end if
       
        x = modpow(B, D, N )        ' x = b^d mod n
        if ( x = bigint_1 or x = Nm1) then
            continue for            ' n passed the initial basenum(i) test
        end if                      '            so skip the squaring part
       
        'run out the possible residues
        for R as integer = 1 to S-1
            x = remainder(x*x, N)
            if x = bigint_1 then    ' if +1 occurs before n-1 then NOT prime
                return(0)
            end if
            if x = Nm1 then         ' n passed this base test; get another base number
                exit for
            end if
        next R
        if x <> nm1 then            ' n-1 NOT FOUND, thus n is COMPOSITE!
            return(0)
        end if
       
    next I                          ' get next base number to test for primality
   
    '' passed all base tests.  The number IS definitely PRIME if n < 2^64
    return(1)                       ' other wise it is higly probable as prime.
   
end function



'' calculation of  base^exponent % modulus;  b^e mod n
'_______
FUNCTION modpow(_
        byref basenumber as bigint, _
        byref exponentnumber as bigint, _
        byref modulusnumber as bigint _
        ) as bigint               
   
    'protect the source values
    dim as bigint b = basenumber
    dim as bigint exponent = exponentnumber
    dim as bigint number = modulusnumber   
   
    'local variables
    dim as bigint result = bigint_1
    dim as bigint quotient, residue
   
    while exponent > bigint_0                           ' b^0 = 1 & DONE
        divmod( exponent, Bigint_2, quotient, residue ) ' exponent \= 2
        if residue = bigint_1 then                      ' exponent was odd
            result = remainder( result * b, number )    ' so multiply result
        end if
        b = remainder ( b*b, number )                   ' square the bas
        exponent = quotient                             ' quotient = exponent \ 2
    wend                                                ' exit when zero
    return( result )
   
end function


'' Generate a small list of prime numbers (ONLY FOR ALTERNATE = 1,2 )
'_______
SUB GeneratePrimeNumberList()
   
    dim as integer count = 0
    dim as bigint candidate = bigint_1
    dim as integer i, j, Limit              'prime number indices   
    PrimeNumber(1) = Bigint_2               ' = 2
    PrimeNumber(2) = Bigint_2 + Bigint_1    ' = 3   
    j = 2                                   ' # of primes in the list
   
    if AlternateTest = 1 then limit = 12  : end if
    if AlternateTest = 2 then limit = 127 : end if
   
    candidate = 3
    do       
        candidate += Bigint_2               ' next ODD number
        for i = 1 to j                      ' result = candidate mod SOMEPRIME
           if remainder( candidate, PrimeNumber(i) ) = 0 then exit for ' not PRIME
        next i                              ' at exit i = j+1
        if i > j then                       ' save the number, it is a prime
            j += 1
            PrimeNumber(j) = candidate
            'print j,PrimeNumber(j)
        end if
    loop while ( j <= Limit )
   
end sub
frisian
Posts: 236
Joined: Oct 08, 2009 17:25

Re: Miller-Rabin primality test

Postby frisian » Nov 23, 2012 17:32

Hi integer,

Nice program you made, using Richard's Big_Integer program.
I mostly GMP if I need to use big number but for my Miller-Rabin routine I decided to use only FB, GMP has a build in Miller-Rabin routine, using that would make it to easy for me.

But I have to point out some error with your program, it will not work correctly if the number to test is (N) is smaller then the witness (B).

I changed N to 207 and insert in the line if t = 1 then primecount += 1 : end if ' I'm sure it's one prime number before the end if a : print N, statement.

This gave the following output.

Code: Select all

  starting with : +207
+199          +197          +191          +181          +179
+173          +167          +163          +157          +151
+149          +139          +137          +131          +127
+113          +109          +107          +103          +101
+97           +89           +83           +79           +71
+67           +61           +59           +53           +47
+43           +41           +37           +31           +29
+23           +17           +11           
    ending with : +7
odd #'s only
    PRIMES found:  38
composites found:  62
 seconds elapsed:  0.0354757998152877

It's clearly to see that at the low numbers the primes 13 and 19 are missing. (two others are missing also)

Simply insert a test to see if N <= B, if this is true then N is a prime.
integer
Posts: 375
Joined: Feb 01, 2007 16:54
Location: usa

Re: Miller-Rabin primality test

Postby integer » Nov 25, 2012 18:27

@frisian
Good catch.

There were no statements, from the web site that, that imposed additional constraints or restrictions.
With those 7 bases, you must also, check to see if the number being tested is:
2, 3, 5, 13, 19, 73, 193, 407521, 299210837 <--- those 9 values are factors of some of the 7 bases,
which give (fill in your own adjective) _________ results.
frisian
Posts: 236
Joined: Oct 08, 2009 17:25

Re: Miller-Rabin primality test

Postby frisian » Dec 01, 2012 18:59

A improved version of the Miller-Rabin program.
It can handle numbers up to 9223372036854775808 (2^63).

Code: Select all

#Define prob_prime -1
#Define composite 0

Function mul_mod(a As ULongInt, b As ULongInt, modulus As ULongInt) As ULongInt
  ' return a * b mod modulus

  Dim As ULongInt x , y = a ' a mod modulus, but a is already smaller then modulus

  While b > 0
    If (b And 1) = 1 Then
      x = (x + y) Mod modulus
    EndIf
    y = (y Shl 1) Mod modulus
    b = b Shr 1
  Wend
  Return x

End Function

Function pow_mod(b As ULongInt, power As ULongInt, modulus As ULongInt) As ULongInt
  ' return b ^ power mod modulus

  Dim As ULongInt x = 1

  While power > 0
    If (power And 1) = 1 Then
      ' x = (x * b) Mod modulus
      x = mul_mod( x,b,modulus)
    EndIf
    ' b = (b * b) Mod modulus
    b = mul_mod(b,b,modulus)
    power = power Shr 1
  Wend
  Return x

End Function

Function prob_prime_test(n As ULongInt, k As Integer) As Byte

  If n > 9223372036854775808ull Then ' limit 2^63, pow_mod/mul_mod can't handle bigger numbers
    Print "number is to big, program will end"
    Sleep
    End
  EndIf
  ' 2 is a prime, if n is smaller then 2 or n is even then n = composite
  If n = 2 Then Return prob_prime
  If (n < 2) Or ((n And 1) = 0) Then Return composite

  Dim As ULongInt a, x, n_one = n - 1, d = n_one
  Dim As UInteger s

  While (d And 1) = 0
    d = d Shr 1
    s = s + 1
  Wend

  While k > 0
    k = k - 1
    a = Int(Rnd * (n - 2)) + 2          ' 2 <= a < n
    x = pow_mod(a, d , n)
    If (x = 1) Or (x = n_one) Then Continue While
    For r As Integer = 1 To s-1
      x = pow_mod(x, 2, n)
      If x = 1 Then Return composite
      If x = n_one Then Continue While
    Next
    If x <> n_one Then Return composite
  Wend
  Return prob_prime

End Function
' -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
' MAIN
' -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

Randomize Timer

Dim As ULongInt y, total, limit = 9223372036854775807

For y = limit To limit-1000 Step -1
  If prob_prime_test(y,10) = prob_prime Then
    total = total + 1
    Print y,
  EndIf
Next

Print
Print
Print total;" primes between ";limit+1;" and ";y

Print
Print "Hit any key to stop"
Sleep
End


/'
23 primes between 9223372036854775808 and 9223372036854774806
found by using a trial division program
9223372036854775783         9223372036854775643         9223372036854775549
9223372036854775507         9223372036854775433         9223372036854775421
9223372036854775417         9223372036854775399         9223372036854775351
9223372036854775337         9223372036854775291         9223372036854775279
9223372036854775259         9223372036854775181         9223372036854775159
9223372036854775139         9223372036854775097         9223372036854775073
9223372036854775057         9223372036854774959         9223372036854774937
9223372036854774917         9223372036854774893
'/

Return to “General”

Who is online

Users browsing this forum: No registered users and 3 guests