## Miller-Rabin primality test

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

### Miller-Rabin primality test

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 0Function 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 xEnd FunctionFunction 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_primeEnd Function' -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-' MAIN' -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-Randomize timerDim As UInteger total, limit = 100000000For y As UInteger = 1 To limit  If prob_prime_test(y,5) = prob_prime Then    total = total + 1  EndIfNextPrintPrint total;" primes found below ";limit+1SleepEnd`
integer
Posts: 378
Joined: Feb 01, 2007 16:54
Location: usa

### Re: Miller-Rabin primality test

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 ProbableDECLARE FUNCTION MillerRabin64(byref anumber as bigint ) as byte        'as close to boolean as possible hereDECLARE FUNCTION modpow(byref basenum as bigint, byref exponentnum as bigint, byref modulusnum as bigint ) as bigintDECLARE 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: 249
Joined: Oct 08, 2009 17:25

### Re: Miller-Rabin primality test

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 : +7odd #'s only    PRIMES found:  38composites 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: 378
Joined: Feb 01, 2007 16:54
Location: usa

### Re: Miller-Rabin primality test

@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,
frisian
Posts: 249
Joined: Oct 08, 2009 17:25

### Re: Miller-Rabin primality test

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 0Function 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 xEnd FunctionFunction 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 xEnd FunctionFunction 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_primeEnd Function' -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-' MAIN' -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-Randomize TimerDim As ULongInt y, total, limit = 9223372036854775807For y = limit To limit-1000 Step -1  If prob_prime_test(y,10) = prob_prime Then     total = total + 1    Print y,  EndIfNextPrintPrintPrint total;" primes between ";limit+1;" and ";yPrintPrint "Hit any key to stop"SleepEnd/'23 primes between 9223372036854775808 and 9223372036854774806found by using a trial division program9223372036854775783         9223372036854775643         92233720368547755499223372036854775507         9223372036854775433         92233720368547754219223372036854775417         9223372036854775399         92233720368547753519223372036854775337         9223372036854775291         92233720368547752799223372036854775259         9223372036854775181         92233720368547751599223372036854775139         9223372036854775097         92233720368547750739223372036854775057         9223372036854774959         92233720368547749379223372036854774917         9223372036854774893'/`