Miller-Rabin test

(Michael Rabin and Gary Miller)
Miller-Rabin primality test is to determine whether a given number is composite or not, so if the result is not composite then the number is called pseudoprime which means it might be a prime and further test to determine primality need to be done. Miller-Rabin test is based on Fermat test.
Algorithm (pseudocode):

1 Choose a random number P to test for primality.
2 Calculate B and M, M is odd, defined as P = 1 + ((2^B) * M).
3 Define random values R, 1 < R < P-1, these random values are witnesses to verify that P is NOT composite, the more of these random value (witnesses) the better the result will be BUT the operation will be more slow, 30 random values are enough.
4 Calculate the first Y = R^M (mod P).
5 If Y = 1 or Y = P-1 then continue with different R.
6 Else set J = 1 to B.
7 calculate Y = Y^2 (mod P).
8 If Y = P-1 then continue with different R.
9 If after the loop finished no Y = P-1 then P is COMPOSITE.
10 If all the random values pass the test then P is probable prime (pseudoprime).

Below is the C++ code, this is only one simple implementation which only handle 32-bits value, remember to put more randomWitness if the "prime" value to be tested is bigger when using big number integer like in OpenSLL or other or mine (yes, I have the big integer class for C/C++ on my web site as well):

BOOL MillerRabinPrimeTestDWORD(DWORD prime)
{

    // randomWitness = witness that the "prime" is NOT composite
    // 1 < randomWitness < prime-1
    DWORD totalWitness = 15;
    DWORD randomWitness[] = { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47};
    DWORD primeMinusOne = prime - 1;
    DWORD oddMultiplier;
    DWORD bitLength;
    DWORD i, j;
    DWORD result;

    // find oddMultiplier, defined as "prime - 1 = (2^B) * M"
    // get bitLength (B) and find the oddMultiplier (M)

    // init value multiplier
    oddMultiplier = primeMinusOne;

    bitLength = 0; // reset
    while(!(oddMultiplier & 1))
    { // it is even

      oddMultiplier >>= 1; // right shift
      bitLength++;

    }

    for(i = 0; i < totalWitness; i++)
    {

    // init value of result = (randomWitness ^ oddMultiplier) mod prime
    result = ExpModDWORD(randomWitness[i], oddMultiplier, prime);

    // is y = 1 ?
    if(result == 1)

      continue; // maybe prime


      // is y = prime-1 ?
      if(result == primeMinusOne)

        continue; // maybe prime

      // loop to get AT LEAST one "result = primeMinusOne"
      for(j = 1; j <= bitLength; j++)
      {

        // square of result
        result = ExpModDWORD(result, 2, prime);

        // is result = primeMinusOne ?
        if(result == primeMinusOne)
          break; // maybe prime
      }

      if(result != primeMinusOne)
        return(TRUE); // COMPOSITE
    }

    // it may be pseudoprime/prime
    return(FALSE);
}

DWORD ExpModDWORD(DWORD value, DWORD exp, DWORD mod)
{

    // use ULONGLONG 64bits to prevent overflow when calculate 32bit x 32bit
    ULONGLONG ullResult = 1;
    ULONGLONG ullValue = value;

    while(exp)
    {
      if(exp & 1)
      { // odd
        // result = (result * value) % mod
        ullResult *= ullValue;
        ullResult %= mod;
      }

      // value = (value * value) % mod
      ullValue *= ullValue;
      ullValue %= mod;

      // exp = exp / 2
      exp >>= 1;

    }

    return((DWORD) ullResult);
}

The odds of a composite passing decreases faster with this test than with previous ones. Three-quarters of the possible values of a are guaranteed to be witnesses. This means that a composite number will slip through t tests no more than (1/4 * T) of the time, where T is the number of iterations. Actually, these numbers are very pessimistic. For most random numbers, something like 99.9 percent of the possible a values are witnesses.
Author Site Map Disclaimer
HMaxF Ultimate Recursive Lossless Compression Research
2001 - 2003 (c) All Rights Reserved.